Epigenomic study identifies a novel mesenchyme homeobox2-GLI1 transcription axis involved in cancer drug resistance, overall survival and therapy prognosis in lung cancer patients

INTRODUCTION

Lung cancer remains the leading cause of death due to malignant disease in both the U.S.A. and worldwide, with 1.6 million cases annually [1, 2, 3, 4]. Non-small-cell lung carcinomas (NSCLCs) have been associated with exposure to carcinogenesis-promoting environmental risk factors, including tobacco smoke [5], as well several genetic and epigenetic aberrations in the lung [6]; however, notably, the number of lung malignancy-related deaths that are not tobacco-associated is continuously increasing [2, 6].

Several molecular mechanisms and cell-signaling pathways that may or may not be induced by exposure to environmental risk factors contribute to lung tumor biology, as reflected in the “Hallmarks of Cancer,” including genetic [8], transcriptional [9], and epigenetic aberrations in the cancer epigenome that are involved in [10–13] progression, survival and/or therapeutic responses in lung cancer patients [14–16].

Several studies have described epigenome-wide aberration profiles in lung cancer patients; some have identified DNA hypermethylation profiles at gene promoter sequences for transcription factors, while others have characterized the homeobox-related genes (HOX) MSX1, OTX1, OSR1, IRX2, PAX6 [17], SIX, LHX, PAX, DLX and ENGRAILED [18], as well HOXA cluster genes such as HOXA7 and HOXA9 [19]. Notably, some HOX gene transcription factors have been proposed to be potential biomarkers for early diagnosis and/or to monitor treatment outcomes in lung cancer patients [20].

In NSCLC, certain mesenchyme-HOX (MEOX) family genes, such as MEOX2, have previously been associated with histopathological progression, poor clinical prognosis and oncological chemoresistance [21]. However, the mechanisms associated with HOX-related genes such as MEOX2 in the context of lung cancer drug resistance, overall survival, and clinical prognosis has yet to be fully elucidated.

In this study, we investigated the lung cancer epigenome to profile target gene promoters that are occupied and likely regulated by the transcription factor HOX-related gene MEOX2 accompanied by transcriptionally active RNA Pol II and epigenetic activation histone markers in human solid lung carcinomas. Bioinformatics analysis allowed us to identify a molecular signature consisting of 78 gene promoters (FDR≤0.2), particularly the Hedgehog-GLI1 gene promoter (with the stringent statistical values FDR≤0.1 and FDR≤0.2), accompanied by the enrichment of active RNA Pol II. In addition, quantitative validation and functional analyses confirmed that expression of the MEOX2-GLI1 transcriptional axis, accompanied by active RNA Pol II and the epigenetic activation histones H3K27Ac and H3K4me3, were involved in resistance to the cancer drug cisplatinum and lung cancer cell migration and proliferation. We also demonstrated that the MEOX2-GLI1 axis is clinically associated with poorer overall survival in lung cancer patients with both Epidermal Growth Factor Receptor (EGFR)-non-mutated and EGFR-mutated status. Thus, for the first time, we have described a novel mechanism by which the MEOX2-GLI1 axis is involved in lung cancer malignancy and platinum-based therapy resistance, as well as an association with clinical progression and poor overall survival, in both EGFR-non-mutated and EGFR-mutated lung cancer patients receiving EGFR-tyrosine kinase inhibitor (TKI)-based therapies.

RESULTS

MEOX2 immunoprecipitation on promoter tiling arrays reveals new MEOX2 gene promoter targets in lung cancer patients

The primary goal of this study was to identify gene promoter targets of the HOX-related gene transcription factor MEOX2 in human solid lung carcinomas, to identify and propose new mechanisms, e.g., transcriptional axes, involved in lung oncological treatment resistance, and to characterize the association between clinical progression and overall survival in lung cancer patients. To address these questions, we implemented an epigenomics approach involving immunoprecipitation assays targeting MEOX2 and the RNA Pol II active enzyme using fragmented chromatin derived from 13 solid lung adenocarcinomas obtained from NSCLC patients ranging in age from 62 to 74 years, with clinical outcomes identified as lower (1 to 16 months) or higher (upper to 70 months) overall survival (Table I). The majority of cases demonstrated partial responses to first/second-line oncological treatment regimens based on cisplatinum/paclitaxel-navelbine and were selected to isolate immunoprecipitated DNA (IP-DNA) in MEOX2 and RNA Pol II immunoprecipitation assays that was used for subsequent competitive hybridization on NimbleGen promoter tiling sequence arrays, following a pipeline strategy descripted in Supplementary Figure 1A. Through bioinformatics analysis, we observed the enriched MEOX2 occupation of gene promoter sequences (Log2 above 532 nm fluorescence) versus RNA Pol II occupation (Log2 below 635 nm fluorescence) in the lung cancer patients designated P-13, P-6 and P-3 (Figure 1A). Differential MEOX2/RNA Pol II occupancy throughout the epigenome at gene promoter sequences was likely accompanied by the enrichment of active versus repressive epigenetic histone markers, as determined when we compared the epigenome of solid lung adenocarcinomas (identified in the present work) with the previously reported lung cancer epigenome from the lung cell line A549 obtained from the ENCODE project database (Figure 1B).

Additionally, bioinformatics analysis of 1,000-bp segments using the fluorescence index Log2-Ratio (MEOX2_532/RNA Pol II_635) identified a total of 6,340, 3,707 and 2,512 fluorescence peaks in the lung cancer patients P-13, P-6 and P-3, respectively, with statistical significance at FDR≤0.1 (Supplementary Figure 1B), representing a total of 5,419, 3,486 and 2,048 gene promoters in each patient, respectively (Figure 2A). Furthermore, 8,440, 5,577 and 4,317 fluorescence peaks were statistically significant at FDR≤0.2 (Supplementary Figure 1B), corresponding to a total of 6,787, 5,536 and 3,347 gene promoters in patients P-13, P-6 and P-3, respectively (Figure 2A).

Importantly, these findings suggest a statistically significant molecular signature consisting of 13 (FDR≤0.1) and/or 78 (FDR≤0.2) MEOX2 gene target promoters (Figure 2A; listed in Supplementary Table I) common to patients P-13, P-6 and P-3 despite their differing clinical survival data (Table I), likely controlling and/or modulating gene expression for cell signaling pathways.

MEOX2 gene promoter targets are associated with cellular and embryonic differentiation signaling pathways in lung cancer patients

Using the MEOX2 gene promoter target data (Log2-Ratio MEOX2/RNA Pol II, FDR≤0.2), we predicted biological, molecular and cellular signaling pathways involved in embryonic development, protein-protein binding, ion binding, nucleic acid and nucleotide synthesis processes. Additionally, we investigated several cell-signaling pathways, such as TGF-b, WNT, NOTCH, MAPK, Jak-STAT, and mTOR, and Sonic Hedgehog (SHH) (Figure 2B), for which association with the Hedgehog-GLI1 gene promoter was statistically significant using both stringent FDR≤0.1 and FDR≤0.2 values (Figure 2A and Supplementary Table I). Interactome prediction analysis revealed a high probability of protein interactions among MEOX2 and protein nuclear factors such as the transcription factors PAX1, PAX3, RELA, and ZEB2, as well PCNA and CDKN1A-2A, which are involved in DNA replication and cell cycle control, respectively (Figure 2C).

MEOX2 and RNA Pol II are positioned throughout GLI-1 gene promoter sequences

Bioinformatics analyses of GLI-1 gene promoter sequences between -3,200 bp and +800 bp were performed by evaluating the fluorescence absolute intensities obtained from individually imprinted probes and determining the fluorescence ratio of the MEOX2_532 ”upper” and RNA Pol II_635 “below” channels (first track, Figure 3A).

Second-track (Figure 3A) fluorescence peaks were based on the fluorescence probe ratios (MEOX2/RNA Pol II Ratio) for the patients P-13 and P-6, which demonstrated fluorescence intensity peaks at FDR≤0.05 (red segments), while patient P-3 exhibited peaks (orange segments) at FDR≤0.1, primarily located from -3,200 to +1 bp in the GLI-1 gene promoter (Figure 3A).

Similarly, as shown in Figure 3A (third track), we identified MEOX2 fluorescence peaks in the GLI-1 gene promoter demonstrating a range of significance levels at FDR≤0.05, FDR≤0.1 and FDR≤0.2 (red, orange and gray color segments, respectively). Interestingly, in the fourth track, the significant presence of RNA Pol II on GLI-1 gene promoter sequences suggested that the GLI-1 gene is transcriptionally active, given the presence of the MEOX2 transcription factor, as well as epigenetically altered, demonstrated by a histone activation profile involving H3K4me3 and H3K27Ac and active enzyme RNA Pol II, which were previously identified in lung adenocarcinoma cells A549 obtained from the ENCODE project database analyzed here (Figure 3B).

The GLI-1 gene promoter is occupied by MEOX2 and active RNA Pol II and is epigenetically associated with the activation histones H3K27Ac and H3K4me3

The bioinformatics results suggested the physical occupation and positive regulation of GLI-1 gene promoter sequences by the transcriptional factor MEOX2. To confirm this hypothesis, we quantitatively validated the MEOX2 occupancy of GLI-1 gene promoter sequences in 13 lung samples from patients with adenocarcinoma (Table I) and in the lung adenocarcinoma cell lines A427 and A549.

We designed an oligonucleotide set (Supplementary Table II) that covered 7 of 12 “TAATTA” MEOX putative transcription binding sites (spanning a sequence promoter region from -2,192 to -109 bp) identified via in silico prediction analysis using the TF-Search online program Computational Biology Research Consortium (see lower section of Figure 3A), and additionally we predicted 29 MEOX2 putative domains between -3,200 and -200 bp on the GLI-1 gene promoter using the bioinformatics program Jaspar (data not shown).

Quantitative qPCR analysis confirmed MEOX2 (Figure 4A) and RNA Pol II (Figure 4B) occupancy on GLI-1 gene promoter sequences in the region spanning -2,192 to -109 bp, and we detected higher levels of MEOX2 and RNA Pol II, particularly at MEOX2-putative 6 and 3 transcription binding sites (Figure 4A–4B).

MEOX2 and RNA Pol II occupancy on the GLI-1 gene promoter was differentially enriched, likely correlating with different clinical parameters such as overall survival time and/or cancer drug response (Table I). We identified significant (p≤0.05) quantitative enrichment in domain 3 compared to domain 6 of the GLI-1 gene promoter, including in the lung adenocarcinoma cell lines A427 and A549 (Figure 4C–4D). These findings indicate that MEOX2 occupancy is compatible with positive transcriptional modulation accompanied by enrichment of the epigenetic activation histones H3K27Ac and H3K4me3 (Figure 4E–4F), specifically the significant enrichment of H3K27Ac (Figure 4G) and H3K4me3 (Figure 4H) in domain 3 compared with 6 domain (p≤0.05), including in the A427 and A549 cell lines.

In contrast, decreased enrichment of the repressive histones H3K27me3 and H3K9me3 in both lung adenocarcinoma patients and the lung adenocarcinoma cell lines A427 and A549 was also detected on GLI-1 gene promoter (data not shown); instead an increased active histone profiling may also be confirmed for GLI-1 gene and promoter sequences in A549 lung cancer cells based on the bioinformatics analyses obtained from the ENCODE project database (Figure 3B). These results strongly suggest that positive transcriptional and epigenetic regulation of the MEOX2-GLI1 axis is involved in malignant lung tumor biology and cancer drug resistance versus responses, as previously suggested for MEOX2 in human lung cancer [21]; however, the mechanisms involved remain uncharacterized, and it is possible that a previously undescribed MEOX2 protein interaction on the GLI-1 gene promoter is responsible for lung cancer drug resistance and/or lung oncological therapy responses.

Inducible mRNA and protein GLI-1 expression occur in a MEOX2- and cisplatinum dose-dependent manner in lung cancer

We observed a MEOX2-dependent reduction in GLI-1 mRNA expression both in the absence and the presence of cisplatinum and confirmed this observation using a specific anti-MEOX2 siRNA cocktail (Figure 5A). Specifically, reduced MEOX2 protein expression following the application of anti-MEOX2 siRNAs correlated with reduced GLI-1 protein co-expression levels (Supplementary Figure 2A).

Furthermore, an inducible, cisplatinum dose-dependent MEOX2 protein expression pattern was identified in lung adenocarcinoma cells (Supplementary Figure 2B), corresponding to cisplatinum-induced GLI-1 protein expression in the lung adenocarcinoma cell lines A427 and A549 that was markedly reduced by anti-MEOX2 siRNAs (Figure 5B).

To demonstrate the probable contribution of the MEOX2-GLI1 axis to chemoresistance, we performed cellular cytotoxic assays to evaluate cisplatinum dose responses using the lung adenocarcinoma cells lines A427 and A549. First, inhibitory concentrations (IC:50) were calculated by constructing cisplatinum dose-dependent curves; the IC:50 was approximately 30 μM for both cell lines (Figure 6A). Next, genetic silencing was performed using siRNA-MEOX2, and we observed a reduction in the IC:50s for both A427 (10.44 uM) and A549 (16.75 uM) cells. Similar effects were seen with siRNA-GLI1 or a mixture of siRNA-MEOX2 plus siRNA-GLI1 (Figure 6B), but no significant reduction in cellular viability occurred with either anti-MEOX2 or anti-GLI-1 siRNAs in an EGFR-mutated lung cancer cell line H1975 with non-detectable MEOX2 mRNA expression (Supplementary Figure 2C, Figure 6C). We concluded that MEOX2 transcription factor occupancy on the GLI-1 gene promoter acts as a positive modulator of MEOX2-GLI1 axis expression via a cisplatinum-dependent resistance mechanism associated with lung tumor cellular viability, cellular proliferation and migration.

The inducible MEOX2-GLI-1 axis is involved in cellular migration and cellular proliferation in lung cancer cells

Cellular migration and proliferation were systematically analyzed in functional assays, and significant (p≤0.001) MEOX2-dependent GLI-1 protein expression was present in A549 cells, which contrasted with non-detectable MEOX2 levels in the H1975 lung adenocarcinoma cell line (Figure 7A and Supplementary Figure 2C). Moreover, reductions in lung tumor cellular migration were observed in the A427 (Supplementary Figure 3A) and A549 (Supplementary Figure 3B) cell lines in an in vitro scratch-migration horizontal assay, and MEOX2 and/or GLI-1 gene silencing resulted in a significant reduction (p≤0.01 and p≤0.001) of the migration index in A549 cells (Supplementary Figure 3B), compared with a lower migration index in A427 cells (Supplementary Figure 3A). We also validated the role of the MEOX2-GLI1 axis by performing transwell-migration vertical assays with both the A549 and NH2347 lung adenocarcinoma cell lines and observed significant differences (p≤0.001) compared with H1975 lung cancer cells (Figure 7B), which express non-detectable MEOX2 levels (Figure 7A and Supplementary Figure 2C).

Finally, the MEOX2-GLI1 axis was validated by investigating the proliferation capacity of the A549 and NH2347 lung cancer cell lines in clonogenic assays, and we identified statistically significant differences (p≤0.01, and p≤0.001) compared with H1975 lung cancer cells (Figure 7C). All of these data are likely associated with clinical stage progression, malignant capacity (invasion/metastasis), and/or cancer drug resistance as well as overall survival in lung cancer patients.

In Situ positive co-expression of MEOX2 and GLI-1 supports a MEOX2-GLI1 axis in two cohorts of lung cancer patients

Clinical validation analyses were performed in two independent lung cancer patient cohorts to identify positive nuclear co-expression of the MEOX2 and GLI-1 proteins in situ in solid lung carcinomas derived from 20 lung cancer patients in the “INER cohort” or 90 lung cancer patients in the “INCAN cohort” in the context of clinical outcomes (Table II), supporting our hypothesis regarding the existence of a novel co-regulatory axis involving association of the MEOX2 protein with GLI-1 gene promoter sequences in human NSCLC.

In this context, double immunostaining analyses of cellular nuclear area expressed as an average percentage (MEOX2-DAB “brown” versus GLI1-AP “red”) were carried out for the INER patient cohort, and proportional positive co-expression levels for the MEOX2 and GLI-1 proteins (i.e., high MEOX2 expression levels correlate with high GLI-1 expression levels) were observed with p≤0.005 and p≤0.0001 significance levels, respectively (Figure 8A). Subsequently, this positive co-expression among both low and high protein expression lung cancer patient groups was quantitatively detected in the INCAN patient cohort by undertaking an Intensity Index analysis based on Algorithm Positive Pixel Bin Counts for MEOX2 and GLI-1 (Figure 8B). In addition, robust quantitative analyses were performed based on the Algorithm Intensity Index plus the average nuclear area percentage, confirming a positive correlation in terms of MEOX2-GLI1 axis co-expression levels for each lung cancer patient at both the p≤0.05 and p≤0.0001 significance levels (Figure 8C). Furthermore, based on clinical outcomes for the INCAN lung cancer patient cohort, which predominantly received platinum-based therapy (Table II), a multivariate analysis emphasized MEOX2 overexpression as a significant hazard risk, with a significance level of p=0.046 (Table III); these results support the existence of a role for the MEOX2-GLI1 axis in clinical malignancy and/or cancer drug-based therapy resistance versus response, reflected in overall survival.

Overexpression of the MEOX2-GLI1 axis is associated with poorer overall survival and cancer drug therapy responses

Based on the analysis described above, a Kaplan-Maier survival curve analysis was constructed for the INCAN cohort, and better overall survival (statistically significant at p≤0.005) was observed for all lung carcinoma patients with lower MEOX2 and GLI-1 protein co-expression levels (statistically significant at p=0.009 and p=0.002, respectively) (Figure 9A). Additional Kaplan-Maier survival curve analyses of lung adenocarcinoma patients representing 84% (76/90) of the total INCAN lung cancer patient cohort indicated poorer overall survival (statistically significant at p≤0.005) in lung adenocarcinoma patients with higher MEOX2 and GLI-1 protein co-expression levels (p=0.079 and p=0.002, respectively) (Figure 9B).

Furthermore, survival curve analysis of all INCAN lung carcinoma patients with EGFR-non-mutated status identified a better overall survival index (statistically significant at p≤0.005) with lower MEOX2 and GLI-1 protein expression levels (p=0.122 and p=0.002, respectively) (Figure 9C). However, survival analysis of only lung adenocarcinoma patients with EGFR-non-mutated status identified poorer overall survival (statistically significant at p≤0.005) with higher MEOX2 and GLI-1 protein co-expression levels (p=0.045 and p=0.006, respectively) (Figure 9D).

Finally, overall survival analysis of all INCAN lung cancer patients with EGFR-mutated status identified better overall survival prognosis with lower MEOX2 and GLI1 protein co-expression levels (p=0.024 and p=0.011, respectively) (Supplementary Figure 4A-4B), which is likely associated with combined-cisplatinum plus TKI-based treatments for a specific group (33/90) of EGFR-mutated lung cancer patients (Table II). In summary, our results provide evidence for a novel transcriptional-epigenetic MEOX2-GLI1 axis that is involved in cellular viability, cancer drug resistance, cellular proliferation and migration capacity, and overall clinical survival and therapy prognosis in lung cancer.

DISCUSSION

Genomic and epigenomic studies have described several genetic and epigenetic aberrations, some of which are identified as Homeobox-related gene transcription factors such as the HOXA and HOXB cluster genes, in human lung carcinomas [17, 18] that are involved in lung cancer biology [19] and are associated with early lung cancer diagnosis [20] and/or malignant stage upon clinical diagnosis [17]. Recently, MEOX2 has been implicated in fetal lung development and histopathological lung cancer progression [22], and MEOX2 and TWIST1 were previously associated with chemoresistance and lung cancer prognosis [21]. Recent published reports suggest MEOX2 overexpression accompanied by epigenetic events constitutes a probable new cancer drug resistance mechanism that must be considered in the context of lung cancer clinical management [23]. However, the identities of MEOX2 gene promoter targets, as well as related epigenetic profiles and molecular mechanisms involved in cancer drug resistance and clinical lung cancer prognosis, have remained unclear.

Here, we describe for the first time a MEOX2 gene promoter target profile in the context of transcriptionally active RNA Pol II enzyme in lung cancer. Our epigenome-wide data were directly obtained from human solid lung carcinomas and identified a range of 3,347 to 6,787 gene promoters (FDR≤0.2) representing several cell signaling pathways, including TGFb, WNT, Notch, Jak-STAT, MAPK, and mTOR; the involvement of some of these pathways in cancer progression, tumor cellular invasion and cancer metastasis has previously been described [24].

A statistically significant molecular signature consisting of 78 gene promoters (FDR≤0.2) or 13 gene promoters (FDR≤0.1) was shared by different lung carcinoma patients; these promoters included Sonic Hedgehog-GLI-1 gene promoter sequences, which have not been previously been reported to be under MEOX2 transcriptional influence. Importantly, the GLI-1 gene was recently proposed to be involved in ABC transporter-independent oncological resistance, in which the TWIST1 and Snail gene transcription factors are linked to the Sonic Hedgehog pathway and induce tumor-initiating cancer stem-like cells, in both non-solid tumor-derived leukemia cells and solid-epithelial tumor-derived cervical cancer cells [25]. Additionally, we suggest the MEOX2 transcription factor promotes platinum-based treatment resistance at the Hedgehog-GLI1 gene promoter level, in accordance with a previous report describing the involvement of Hedgehog-GLI1 in the ABCG2 transporter-dependent oncological resistance exhibited by a lung cancer stem-like side cell population [26] and is involved maintaining the self-renewal of cancer stem-like cells as recently reported for arsenic trioxide-treated lung carcinoma cells [27]. Here, we suggest not only the positive modulation of Hedgehog-GLI1 gene expression but also positive induction of the mTOR pathway (Figure 2), in line with very recent reports describing an association between the mTOR and 4E-BP1 pathways, which are under TWIST1 transcriptional control, and poorer overall survival in NSCLC patients [28], as well as the involvement of TWIST1-Snail, linked to the Hedgehog pathway, in a chemoresistant cancer stem cell-like phenotype [25].

The autocrine effects of the SHH pathway have previously been linked to cellular proliferation in NSCLC cell lines (A549, H322M, and HOP-62) and solid lung carcinomas derived from patients; specifically, GLI-1 genetic silencing combined with an SHH antagonist reduces lung tumor cellular proliferation [29]. Additionally, GLI-1 gene transcription factors have been functionally related to cancer stem cell self-renewal capacity and tumorigenicity in solid glioma tumors [30]. Furthermore, SHH-GLI1 signaling inhibition is involved in sensitization to EGFR-TKI (gefitinib/erlotinib)-based therapy, thereby reducing cellular viability, as well as the self-renewal function of cancer stem-like cells derived from the EGFR-mutated lung cancer cell lines H1650 and H1975 [31]. Here, we identified significant participation of the MEOX2-GLI1 axis in both cellular cytotoxic resistance capacity and cellular migration proliferation, which is likely explained in part by previous studies showing the involvement of the SHH-GLI1 signaling pathway in cellular proliferation, cellular transmigration and epithelial-mesenchymal transition in the NSCLC cell lines A549 and H520 [32]. In this regard, recent reports have demonstrated a correlation between MEOX2 overexpression and greater cellular migration capacity, increased senescence and cell cycle progression, promoting the development of vessel cells [33].

Several previous reports implicated the GLI-1 gene in docetaxel, methotrexate and etoposide chemoresistance through transcriptional modulation of the ABC transporter family genes in esophageal adenocarcinoma, prostate carcinoma and metastatic squamous cell carcinoma cell lines [34]. Other studies described the contribution of the canonical cell signaling SHH-GLI1 axis to acquired chemoresistance to doxorubicin, vincristine, and etoposide in both non-solid tumor-derived K-562 erythroid leukemia cells and solid tumor-derived epidermoid carcinoma KB cells [35]. These effects have also been observed in NSCLC cell lines under gefitinib/erlotinib exposure [31] and in mouse- and human-derived SCLC cells treated with the cancer drugs etoposide, carboplatin, cisplatinum and/or the SHH-antagonist NVP-LDE225 [36]. Additionally, a significant reduction in cisplatinum chemoresistance in experimental in vitro models is achieved via inhibition of the SHH-SMO-GLI1 axis signaling with an SMO antagonist (GDC-0449) in both human NSCLC and SCLC cell lines [37].

However, few studies have described chemoresistance at the transcriptional level in terms of gene promoter sequences. The ABCG2 transporter gene promoter is modulated by the GLI-1 gene transcription factor through gene sequences in other domains located from -416 to -408, accompanied by increased levels of the histone marker H3K18Ac, following cellular stimulation of the SHH-GLI1 axis [38]. Such epigenetic markers promote GLI-1 occupation of ABCG2 gene promoter sequences and induce the expression of ABCG2, which is involved in chemoresistance in the presence of doxorubicin and methotrexate in non-solid tumors and B-cell non-Hodgkin lymphoma cells [38]. These and other reports indicate the functionality of the SHH-GLI1-ABCG2 axis is highly correlated with chemoresistance, as shown for the lung cancer stem cell-like population derived from H460 parental lung cancer cells [39].

Our results obtained from solid lung carcinomas and the A427 and A549 NSCLC cell lines constitute the first evidence demonstrating significant epigenome occupancy of the homeobox-related gene MEOX2. Additionally, GLI-1 gene promoter sequences are transcriptionally occupied by active RNA Pol II and epigenetically accompanied by the active histones H3K27Ac and H3K4me3 with statistically significant enrichment at two hotspot regions ranging from -1,830 to -1673 and -822 to -665, supporting the existence of an active MEOX2-GLI1 axis, which was also indicated by the lung cancer epigenome analysis derived from A549 lung adenocarcinoma cells (Figure 1B and Figure 3B) previously published and deposited by Sabo PJ and colleagues in the ENCODE Project database [40].

Here, our findings are in accordance with intrinsic chemoresistance and/or also likely associated with acquired chemoresistance, which have been linked to and/or involved in epigenetic mechanisms based on aberrations in the histone repressive profiles H3K27me3 on MEOX2 gene promoter [21], but also histone activation profiles detected for H3K27Ac and H3K4me3 on GLI-1 gene promoter sequences (Figure 3B and Figure 4E-4H). But in other cases is modulated by long non-coding RNAs (lncRNAs), as was recently been reported for lncRNA HOTAIR, which promotes cisplatinum drug resistance by controlling p21 protein expression to increasing cell cycle control in A549 lung cancer cells [41]. Additionally, HOTAIR overexpression has been detected in solid tumor NSCLC patients with poorer prognosis as well as functionally associated with higher lung tumor cellular migration and cellular invasion in vitro lung in A549 and SPC-A1 cancer cells, with lung tumor metastasis in vivo in SPC-A1 lung cancer cells [42], and with poorer prognosis in lung cancer patients [37, 38]. Other possible explanations include a role for HOTAIR-coordinated assembly “binding” to the Polycomb Repressive Complex 2 (PRC2) in a SUZ12/EZH2 and/or LSD1/CoREST/REST complex-dependent manner, which are reportedly involved in increased repressive histone H3K27 methylation and H3K4 demethylation, respectively [43]. Additionally, lncRNA HOTAIR is epigenetically involved in repression of the p15, p21 and p27 genes at the promoter sequence level in lung cancer [44]. Down-regulated expression of lncRNA MEG3 is involved in cisplatinum resistance through p53 and Bcl-xl controlled gene expression in lung carcinoma cells [45], although additional epigenetic evidences must be obtained in lung cancer patients. Additionally, lncRNA BC087858 overexpression is clearly associated with acquired resistance to EGFR-TKI based “gefitinib” therapy [46], while lncRNA GAS5 down-regulation is involved in EGFR-TKI resistance in lung cancer patients and NSCLC cell lines [47]; in both cases, the underlying epigenetic mechanisms have yet to be fully clarified.

Our functional epigenomics strategy also suggests novel MEOX2-GLI1 transcriptional axis signaling is involved in resistance to the cancer drug cisplatinum in a dose-dependent manner, as previously proposed for lung cancer [21]. MEOX2 overexpression resulted in cisplatinum dose-dependent inducible GLI-1 gene expression at the gene promoter level associated with chemoresistance in lung cancer cells, which likely correlates with previous reports in which inhibition of the SHH-GLI1 axis is involved in cancer drug resistance by a lung cancer stem-like cell side population [37] through a direct transcriptional target action in the ABCG2 gene promoter sequence [38]. Additionally, this mechanism is involved in lung adenocarcinoma cell survival and cancer drug resistance conducted via the SHH-GLI1 axis but negatively regulated by Hedgehog-Interacting Protein (HHIP) in EGFR-non-mutated A549 or EGFR-mutated H1975 and HCC827 NSCLC cell lines under TKI-based treatment [48]; however, here this mechanism is positively regulated by MEOX2 at the transcriptional level.

In conclusion, here we have reported for the first time the functional and clinical correlations of the overexpressed MEOX2-GLI1 transcriptional axis, both as a novel molecular candidate axis involved in lung cancer drug resistance and as a predictor of therapy prognosis. Certain reports investigating NSCLC patients have described an overexpressed SHH-GLI2 axis associated with poorer progression-free survival [49] as well as an overexpressed GLI-2 axis involved in cellular proliferation and apoptosis resistance in lung squamous cell carcinomas [50]. However, as determined in SMO inhibition assays, these axes promote the down-regulation of GLI-1 gene expression and facilitate significant EGFR-TKI-based treatment susceptibility, thereby reducing tumor cell proliferation in both EGFR-wild type as well as EGFR-mutated NSCLC cell lines [51], suggesting an interdependent oncogenic addiction between the EGFR and SHH-GLI1 axis pathways that is involved in lung cancer cell survival and cancer drug resistance [52]. We obtained similar results for both EGFR-wild type and EGFR-mutated lung cancer patients (Figure 9C–9D and Supplementary Figure 4A). However a non-canonical GLI-1 activation likely occurs in a Rack1-dependent manner in NSCLC patients as a proposed necessary step for lung cancer tumorigenicity [53]. Additionally, this process is likely involved in EGRF-TKI-based therapy resistance in both EGFR-non mutated (A549) and EGFR-mutated (H1975) lung cancer cells [54].

Despite developing a significant correlation index based on the overexpressed MEOX2-GLI1 axis, our results depict a novel transcriptional and likely epigenetic mechanism involved in lung tumor biology and cancer drug resistance mechanisms to platinum-based and/or target EGFR-TKI-based therapy responses, affecting lung tumor cell viability, cellular migration/proliferation, clinical prognosis and overall survival in lung cancer patients.

MATERIALS AND METHODS

Sample selection and cell lines

A total of 123 tumor samples from 2 cohorts of NSCLC patients were included in the present study. The training cohort consisted of 13 fresh-frozen (FF) lung tumors collected from the Thoracic Surgical Service and 20 formalin-fixed paraffin-embedded (FFPE) lung tumors from the Pathology Service at the National Institute of Respiratory Diseases (INER) in Mexico City. An independent validation cohort consisted of 90 randomly selected patients with advanced NSCLC enrolled in the Thoracic Oncology Unit at the National Cancer Institute (INCAN) in Mexico City from July 2010 to December 2011. We included samples from patients with clinical stages corresponding to the TNM Staging System (6th edition) who had an Eastern Cooperative Oncology Group Scale of Performance Status (ECOG PS) score of 0–2 and were candidates for platinum-based chemotherapy and/or TKI-based target therapy.

In addition, the NSCLC cell lines A427, A549, NH2347 and H1975 from the American Type Culture Collection (ATCC, Manassas VA, USA) were used to perform experimental assays to investigate cisplatinum drug resistance and epigenetic analyses. The Institutional Review Board reviewed and approved protocols for INER (B17-07, B09-08) and INCAN (015/016/ICI, 964/15). All patients had a histologically confirmed lung cancer diagnosis, underwent surgical procedures, and provided informed consent before tumor resection for this study and future genetic/epigenetic biomarker studies.

Clinical data collection, treatment regimen and patient follow-up

Clinical outcome data for all patients, including their clinical history, histopathological diagnosis, performance status and smoking history, were obtained from medical records. Smoking history was defined as ever smokers (combining current or former smokers). Current smokers were defined as those patients smoking who have smoked at least 100 cigarettes in life and were still smoking at the very moment of the interview. Former smokers were defined as those patients who were not smoking at the moment of the interview but have smoked at least 100 cigarettes in their life. All patients were treated according to international guidelines for lung cancer treatment [55]. All received platinum-based chemotherapy as a first-line treatment in combination with paclitaxel or vinorelbine. After progression with cytotoxic chemotherapy, the patients received second- or third-line treatment with an epidermal growth factor receptor (EGFR)-targeting tyrosine kinase inhibitor (TKI), either gefitinib or erlotinib, based on their mutation status. EGFR non-mutated and mutated patients (exons 18–21) were detected using a Therascreen RGQ EGFR Real-time PCR Kit (QIAGEN, Scorpions ARMS method, Dusseldorf, Germany), which was previously used to genotype NSCLC patients [56].

Chromatin immunoprecipitation (ChIP)

FF lung tumor (LT; 3 mm³) samples consisting of almost 70% tumor tissue were pulverized under liquid nitrogen, fixed with 1% formaldehyde for 10 min, and immediately neutralized with 0.125 M glycine for 5 min. The resulting cells were washed with 1X PBS and treated with lysis buffer (50 mM HEPES-KOH pH 7.5, 140 mM NaCl, 1 mM EDTA pH 8, 1% Triton X-100, 0.1% sodium deoxycholate, and 0.1% SDS) containing 10 ml with protease inhibitors (complete Mini, Roche, Indianapolis, IN, USA).

The chromatin from the FF-LT samples was sonicated for 10 pulses of 20 seconds w/u 60 watts. Then, 10 μg of chromatin was immunoprecipitated (ChIP) using a commercial EZ-Magna ChIP™ G kit (Millipore, Temecula, CA, USA) and 2.5 μg of anti-MEOX2 (Santa Cruz Biotechnology, Dallas, TX, USA) as well as 1 μg of activated anti-H3K27Ac, anti-H3K4me3, anti-H3K27me3, anti-H3K9me3 or anti-RNA Pol II antibodies (Abcam, Cambridge Science Park, Cambridge, U.K.). A 1-μg aliquot of anti-mouse IgG was used as a negative control for ChIP (Millipore, Temecula, CA, USA).

IP-DNA amplification, labeling and hybridization using promoter sequence tiling arrays

Sequence promoter tiling arrays (RefSeq, NimbleGen Assembly HG18, version 3x720K. Roche NimbleGen, Germany) representing 30,893 transcripts, 22,542 promoter sequences, and enhancer regions were employed to identify binding sites for the MEOX2 transcription factor. Microarrays containing 50-75-bp probes with 100 bp between each probe, representing -3,200 bp upstream and +800 bp downstream of the transcription start site (TSS), were used.

Immunoprecipitated DNA (IP-DNA; 20 ng) was linearly amplified using a Whole Genome Amplification (WGA2) GenomePlex Kit (IP-DNA-WGA2; 500 ng) and Reamplification Kit (WGA3) (Sigma, St. Louis, MO, U.S.A.). In both cases, WGA2 and WGA3 were amplified under the following conditions: Initial denaturation step at 95°C for 3 min, 14 cycles of denaturation/annealing/extension at 94°C for 15 seconds, and a final step at 65°C for 5 min. Re-amplified IP-DNA was purified using the phenol-chloroform method.

Next, a commercial NimbleGen Dual-Color DNA Labeling Kit (Roche, Germany) was used to label 1 μg of IP-DNA obtained with anti-RNA Pol II and anti-MEOX2 via incubation with a solution of nonamers (Cy3 and Cy5, respectively), 10 mM dNTP mix and 50 U/μl exo-Klenow fragment (3’-> 5’ exo-) for 3 hours at 37°C in the dark. Finally, 10 μl of buffer (0.5 M EDTA) was added to stop the fluorescence labeling reaction.

Total labeled IP-DNA was purified and washed with both absolute isopropanol and 70% ethanol and quantified at 260 nm using a spectrophotometer (Bioteck, model EPOCH, USA). Then, 34 μg of previously labeled IP-DNA obtained using both the MEOX2 and RNA Pol II antibodies was mixed and hybridized for 48 hours at 42°C. Finally, the microarrays were washed using the NimbleGen Wash Buffer kit (Roche Germany) and subjected to dry centrifugation (Roche NimbleGen, Germany) for 2 min at 1,500 rpm, and images were captured.

ChIP-on-chip bioinformatic analysis

Microarray images were obtained at 532 nm and 635 nm using a scanner (GenePix model 4000a, USA). In addition, bioinformatics analysis was performed using NimbleScan software version 2.6 (NimbleGen, Germany).

Both images (one from each channel) were subjected to microarray alignment with values close to 0.1290, adjusting for brightness and contrast, and stored using the .tiff format. Next, ChIP-on-chip analysis generated the following bioinformatics reports containing information from both fluorescence channels: Pairs reports (.pair), Probe reports (.probe) and GFF reports (General Feature Format: .gff).

Additionally, a complete genome bioinformatics sub-analysis was performed using the CGH method to show the presence or absence of the fluorescence signals obtained from IP-DNA segments attached to MEOX2 and/or RNA Pol II, indicating possible transcriptionally active MEOX2 gene promoter sequence targets. Later, NimbleGen SignalMap (software version 1.9) was used to map the positions of the signal segments resulting from the hybridization data; for this analysis, a Log₂-Ratio scale was established, and the index (ratio) was obtained using the formula 532/635 (MEOX2_532/RNA Pol II_635). Peaks were detected based on fluorescence intensity and determined by the alignment of 4 probes near one another (min probe > cutoff “cutoff peak”), contained in a sliding window of 500 bp.

While the cutoff value was evaluated based on a hypothetical 95 maximum percentage defined as the mean+6 [standard deviation], data from the index (ratio) were randomly assessed 20 times to establish the probability of a positive false discovery rate (FDR); the values FDR≤0.05, FDR≤0.1 and FDR≤0.2 indicate fluorescence peaks. (all epigenome data associated have been deposited in the online Dryad Digital Repository System, with the provisional doi:10.5061/dryad.rm7dd).

Signaling prediction pathways (from MEOX2 gene promoter targets)

Hierarchical data analyses based on FDR probability and the in silico prediction of gene expression levels identified and selected genes that were most likely to be under the transcriptional control of MEOX2 in lung cancer patients. These analyses to predict metabolic pathways or signaling pathways were undertaken using the online KEGG (Kyoto Encyclopedia of Genes and Genomes) database, and in silico prediction analysis was performed using the online webpage http://www.webgestalt.org with the following statistical parameters: a hyper-geometric method and Bonferroni correction tests with a significance level of 0.05, and a minimum number of 2 genes per category.

qPCR assays to quantitatively validate ChIP analyses

The IP-DNAs were analyzed by performing qPCR absolute quantification using a LightCycler 480 System (Roche, Mannheim, Germany), SYBR Green Master Mix (KAPA Science, Foster City, CA, USA) and oligonucleotide sets for the proposed gene promoter target sites (Supplementary Table II).

A total of 20 ng of IP-DNA was used per reaction, followed by the amplification of standard curves through serial dilutions of native DNA (diploid control genomic DNA) derived from the peripheral blood mononuclear cells of healthy donors (100 ng, 10 ng, 1.0 ng, 0.1 ng and 0.01 ng).

The following amplification conditions were used: initial denaturation at 95°C for 10 min and 40 cycles of denaturation at 95°C for 15 seconds, annealing at 55°C for 30 seconds and extension at 72°C for 30 seconds. Oligonucleotides (Supplementary Table II) were designed using the OLIGO and NT VECTOR programs and synthesized by SIGMA-ALDRICH (USA).

Genetic silencing (siRNA transfection assays)

We used the anti-MEOX2 (sc-106233), and GLI-1 (sc-37911) siRNAs, each consisting of a pool (cocktail) of 3 target-specific siRNAs (19-25 nt) designed to knock down corresponding human gene expression. Two non-human related siRNAs obtained from Santa Cruz Biotechnology (Dallas, TX, USA) were used as negative controls (sc-37007 and sc-44230). Briefly, 3×10⁵ cells in RPMI-1640 (Biowest, USA) and 3% antibiotic-free FBS (Biowest, USA) were incubated in 6-well plates for 12 hours and then transfected for a total of 48 hours with Lipofectamine 2000® (Invitrogen, Foster City, CA, USA) and siRNAs at 50 nM in the presence or absence of a cisplatinum drug in a total volume of 2,500 μl.

RT-qPCR assays for mRNA expression analyses

The TRIzol (Invitrogen, California, USA) method was used for total mRNA extraction and purification, and cDNA was synthesized from 2 μg of total RNA using the ImProm-II Reverse Transcription System (Promega, Madison Wisconsin, USA). mRNA expression detection was performed using a Probe Master kit (Roche, Germany) and hydrolysis probes (Universal Probe Library “UPL”, Roche, Germany) on a LightCycler® 480 Real-Time PCR System (Roche, Germany) following the manufacturer’s instructions. mRNA expression levels were normalized to expression of the endogenous gene GAPDH, and expression analyses were performed using the ΔΔT method.

Oligonucleotide sequences and UPL probes (Supplementary Table III) were applied under the following amplification conditions: initial denaturation at 95°C for 5 min and 40 cycles of denaturation at 95°C for 10 seconds, annealing at 60°C for 17 seconds and extension at 72°C for 1 seconds. Oligonucleotides were designed using the ProbeFinder program version 2.51 (Roche, Germany) and synthesized by SIGMA-ALDRICH (USA).

Western blot assays

Total proteins were extracted from cultured cells and purified using the TRIZOL method (Invitrogen), suspended in 5% SDS with protease inhibitors (cOmplete Mini, Roche, Indianapolis, IN, USA) and quantified using a commercial protein assay kit (Bio-Rad DC). Next, 30 μg of total protein was resolved by vertical electrophoresis in 12% acrylamide gels and then transferred to nitrocellulose membranes using Trans-Blot-Turbo equipment/Transfer System Pierce G2 (Thermo Fisher Scientific, Foster City, CA, USA). The membranes were blocked at room temperature for 2 hours with 5% low-fat milk in 1X PBS/0.1% Tween-20 and incubated overnight with antibodies against GLI1 (Rockland, 100-401-223, Pottstown, PA USA; and Abcam AB134906, Cambridge Science Park, Cambridge, U.K.) (1:3,000), MEOX2 Sc-81971 (1:1,500), and GAPDH Sc-25778 (1:3,000) from Santa Cruz Biotechnology (Dallas, TX, USA). The membranes were washed with 1X PBS supplemented with 0.1% Tween-20 at room temperature and then incubated for 2 hours with a secondary antibody (anti-mouse HRP/anti-rabbit HRP; 1:10,000) at room temperature, followed by washing with 1X PBS supplemented with 0.1% Tween-20 at room temperature. The membranes were incubated with Super Signal West Femto Maximum Sensitivity Substrate (Thermo Fisher Scientific, Fermentas Life Science, Foster City, CA, USA), and a C-Digit Scanner was used for image acquisition (LI-COR, Lincoln, NE, USA). The change ratio was calculated using Image Studio Software (version 4.0.21).

Cellular viability (Cytotoxicity) assays

Cellular viability analyses were performed in 96-well plates using 3-(4,5-dimethylthiazol-2-yl)-5-(3-carboxymethoxyphenyl)-2-(4-sulphophenyl)-2H-tetrazolium (MTS) (Promega, Madison, WI, USA) in the presence or absence of cisplatinum. Briefly, 96-well plates were seeded with 2,000 to 3,000 cells in 100 μl of RPMI-1640 culture medium for 48 hours and subsequently incubated with 20 μl of MTS (5 mg/ml) for 3 hours to be read later at 490 nm; cell viability was calculated using the following formula: cell viability = (O.D. of the sample/Control O.D.) × 100. Similarly, drug resistance curves were constructed to evaluate dose responses using serial dilutions of cisplatinum.

Cellular scratch (Migration) assays

Wounding (approximately 800 microns was carried out on cells 8 hours after monolayer cultures were established using 2x10⁵ cells per well (24-well plates)). The cells were deprived of Fetal Bovine serum (FBS) for the remainder of the 48-hours experiment; then, cultures were washed 3 times with 1X PBS. Photographs (40x) were taken at 0, 24, 36 and 48 hours post-wounding. The wound area was quantified using ImageJ, employing the “Duplicate” “Process”, “Binary” and “Make binary” applications to transform pixels into binary code. Finally, adjustments to “Scale Set” and “Analyze Particles” facilitated quantitative representation of the cell migration index expressed as a percentage (%) with respect to the wound (scratch) at time 0 hours.

Cellular transwell (Migration) assays

Six-well plates were seeded in triplicate at a density of 3×10⁵ cells per well. Following overnight adherence, cells were transfected with corresponding siRNAs for 6 hours. Subsequently, 25×10³ transfected cells were transferred into upper transwell chambers in 24-well plates, while the bottom chambers contained RPMI-1640 medium supplemented with 10% FBS. The total culture time was 48 hours at 37°C and 5% CO₂. Then, cells that adhered below the transwell inserts were carefully fixed using 70% methanol, stained with 0.1% crystal violet, washed with running water, and mounted on slides using entellan. Final images were taken with light field microscope at a total magnification of 100X.

Clonogenic assays

Six-well plates were seeded at a density of 3x10⁵ cells per well. After overnight adherence, cells were transfected with corresponding siRNAs. The cells were then transferred in triplicate into 6-well plates at a density of 1x10⁴ per well and into 12-well plates at a density of 5x10³ per well. Cells were cultured for 6 to 21 days at 37 °C and 5% CO₂, followed by fixation with methanol and acetic acid (3:1 v/v), washing in running water, and staining with 0.5% crystal violet. Colony growth was quantified using the ImageJ Colony Number plugin tool [57].

Immunohistochemical assays using histopathological samples

Immunohistochemical analyses were performed on FFPE-LT tissue samples, which were used to construct tissue microarrays (TMAs). Lung cancer specimens were blindly selected by a pathologist who identified triplicate lung tumor regions characterized by high cellularity for each lung tumor. The arrays were constructed using the Chemicon Advanced Tissue Arrayer (ATA100) with a 0.5-mm needle. We used 1-2 μg of specific anti-MEOX2 (Santa Cruz Biotechnology, Dallas, Texas, USA) and anti-GLI-1 antibodies (Rockland, Pottstown, PA USA). UltraView Universal DAB and UltraView Universal Alkaline Phosphatase Red detection kits were used in conjunction with a Ventana System BenchMark GX (Ventana Medical System Inc., Tucson, AZ, USA). The slides were dehydrated, cleared in xylol and mounted in synthetic resin. A Leica DM750 transmitted light microscope (Wetzlar, Germany) was used to obtain photomicrographs at 200X and 400X magnification.

The image analysis software program ImageJ and the ImmunoRatio Plugin (Version 1.42q), developed by Jorma Isola and Vilppu Tuominen at the Institute of Biomedical Technology, University of Tampere (http://153.1.200.58:8080/immunoratio/), were used to quantitatively analyze the percentage corresponding to positive nuclear area (average) on histopathological slides (3 fields for each slide sample) at 200X and 400X. In addition, Aperio ePathology ImageScope Software version 11.0.2.725 (Leica Biosystems, USA) and the Intensity Index based on the Algorithm Positive Pixel Count 2004-08-11 (Version 8.100) were used to categorize positive pixel counts into three Intensity Bins, expressed as the Number of Strong Positive Reactions Index (described as Nsr=Nsp/Nwp+Np+Nsp). All data were statistically analyzed to obtain significant differences among the NSCLC patient groups.

Statistical analyses

T-tests, one-way ANOVA, Tukey’s, Dunn’s and Mann-Whitney´s multiple comparison tests were performed to analyze the differences between patient groups; p values ≤ 0.05 were considered statistically significant. Overall survival analyses were calculated using the Kaplan-Meier method, including comparisons among subgroups with the Breslow-Wilcoxon and Log Rank “Mantel-Cox” tests; p values ≤ 0.01, 0.005 and 0.0025 were considered statistically significant. Adjustments for potential confounders were calculated using a multivariate Cox regression model, and the relative risk (RR) ratio was calculated with a corresponding 95% C.I. as a measure of association. All statistical analyses were carried out using the software programs GraphPad (Mac-OS X, Version 5.0) and IBM-SPSS Statistics (Windows, version 20).

ACKNOWLEDGMENTS

We would like to thank postdoctoral fellow Herrera Solorio Abril Marcela as well as Peralta Alvarez CA, Muñiz Hernández Sae, Velázquez C, Salgado Bernabe E, and Pineda Villegas P for their technical assistance. We also thank Aguilar Lemarroy A and Jave Suárez L for access to the genome hybridization and scanning platform at CIBO-IMSS, Ramírez Tirado L Alejandra for her professional assistance with the lung cancer patient databases and Zatarain Barron Zyanya Lucia for her manuscript edition helpful. This work was funded by the National Institute of Respiratory Diseases (INER) and the INER-Research Unit through the collaboration agreement register: UNAM 43355-3065-17-XI-15; Patricia Piña Sánchez thanks as the recipient of the IMSS-FOUNDATION Research Support Grant. This work fulfills one of the requirements for Leonel Armas López to obtain a Ph.D. degree in Posgrado en Ciencias Biológicas (UNAM): he was the recipient of a doctoral fellowship from CONACyT (Fellowship 226096, CVU 294950).

CONFLICTS OF INTEREST

The authors disclose no potential conflicts of interest.

GRANT SUPPORT

This work was funded by projects PAPIIT-DGAPA: IB202512, RR282512 and IN226317 from the National University Autonomous of Mexico (UNAM), as well as FOSISS 2016-1: Project 00272655 from the National Council for Science and Technology (CONACYT). “The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.”

REFERENCES

1. Jemal A, Bray F, Center MM, Ferlay J, Ward E, Forman D. Global cancer statistics. CA Cancer J Clin. 2011; 61:69–90. doi: 10.3322/caac.20107.

2. Siegel R, Naishadham D, Jemal A. Cancer statistics, 2013. CA Cancer J Clin. 2013; 63:11–30. doi: 10.3322/caac.21166.

3. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2017. CA Cancer J Clin. 2017; 67:7–30. doi: 10.3322/caac.21387.

4. Torre LA, Siegel RL, Jemal A. Lung Cancer Statistics. Lung Cancer and Personalized Medicine [Internet] Cham: Springer; 2016. pp. 1–19. [cited 2017 Aug 16]. doi: 10.1007/978-3-319-24223-1_1.

5. Sun S, Schiller JH, Gazdar AF. Lung cancer in never smokers--a different disease. Nat Rev Cancer. 2007; 7:778–90. doi: 10.1038/nrc2190.

6. Besaratinia A, Pfeifer GP. Second-hand smoke and human lung cancer. Lancet Oncol. 2008; 9:657–66. doi: 10.1016/S1470-2045(08)70172-4.

7. Sun S, Schiller JH, Spinola M, Minna JD. New molecularly targeted therapies for lung cancer. J Clin Invest. 2007; 117:2740–50. doi: 10.1172/JCI31809.

8. Mitchell KR, Warshawsky D. Frequent Ha-ras mutations in murine skin and liver tumors induced by 7H-dibenzo[c,g]carbazole. Mol Carcinog. 1999; 25:107–12.

9. Delgado J, Martinez LM, Sánchez TT, Ramirez A, Iturria C, González-Avila G. Lung cancer pathogenesis associated with wood smoke exposure. Chest. 2005; 128:124–31. doi: 10.1378/chest.128.1.124.

10. Bartke T, Borgel J, DiMaggio PA. Proteomics in epigenetics: new perspectives for cancer research. Brief Funct Genomics. 2013; 12:205–18. doi: 10.1093/bfgp/elt002.

11. Baylin SB. Resistance, epigenetics and the cancer ecosystem. Nat Med. 2011; 17:288.

12. Brena RM, Costello JF. Genome-epigenome interactions in cancer. Hum Mol Genet. 2007; 16:R96–105. doi: 10.1093/hmg/ddm073.

13. Roessler S, Long EL, Budhu A, Chen Y, Zhao X, Ji J, Walker R, Jia HL, Ye QH, Qin LX, Tang ZY, He P, Hunter KW, et al. Integrative genomic identification of genes on 8p associated with hepatocellular carcinoma progression and patient survival. Gastroenterology. 2012; 142:957–966.e12. doi: 10.1053/j.gastro.2011.12.039.

14. Levina V, Marrangoni AM, DeMarco R, Gorelik E, Lokshin AE. Drug-selected human lung cancer stem cells: cytokine network, tumorigenic and metastatic properties. PLoS One. 2008; 3:e3077. doi: 10.1371/journal.pone.0003077.

15. Chen J, Emara N, Solomides C, Parekh H, Simpkins H. Resistance to platinum-based chemotherapy in lung cancer cell lines. Cancer Chemother Pharmacol. 2010; 66:103–11. doi: 10.1007/s00280-010-1268-2.

16. Kratz JR, Yagui-Beltrán A, Jablons DM. Cancer stem cells in lung tumorigenesis. Ann Thorac Surg. 2010; 89:S2090–5. doi: 10.1016/j.athoracsur.2010.03.038.

17. Rauch TA, Zhong X, Wu X, Wang M, Kernstine KH, Wang Z, Riggs AD, Pfeifer GP. High-resolution mapping of DNA hypermethylation and hypomethylation in lung cancer. Proc Natl Acad Sci USA. 2008; 105:252–7. doi: 10.1073/pnas.0710735105.

18. Pfeifer GP, Rauch TA. DNA methylation patterns in lung carcinomas. Semin Cancer Biol. 2009; 19:181–7. doi: 10.1016/j.semcancer.2009.02.008.

19. Rauch T, Wang Z, Zhang X, Zhong X, Wu X, Lau SK, Kernstine KH, Riggs AD, Pfeifer GP. Homeobox gene methylation in lung cancer studied by genome-wide analysis with a microarray-based methylated CpG island recovery assay. Proc Natl Acad Sci USA. 2007; 104:5527–32. doi: 10.1073/pnas.0701059104.

20. Nawaz I, Qiu X, Wu H, Li Y, Fan Y, Hu LF, Zhou Q, Ernberg I. Development of a multiplex methylation specific PCR suitable for (early) detection of non-small cell lung cancer. Epigenetics. 2014; 9:1138–48. doi: 10.4161/epi.29499.

21. Ávila-Moreno F, Armas-López L, Álvarez-Moran AM, López-Bujanda Z, Ortiz-Quintero B, Hidalgo-Miranda A, Urrea-Ramírez F, Rivera-Rosales RM, Vázquez-Manríquez E, Peña-Mirabal E, Morales-Gómez J, Vázquez-Minero JC, Téllez-Becerra JL, et al. Overexpression of MEOX2 and twist1 is associated with H3K27me3 levels and determines lung cancer chemoresistance and prognosis. PLoS One. 2014; 9:e114104. doi: 10.1371/journal.pone.0114104.

22. Cortese R, Hartmann O, Berlin K, Eckhardt F. Correlative gene expression and DNA methylation profiling in lung development nominate new biomarkers in lung cancer. Int J Biochem Cell Biol. 2008; 40:1494–508. doi: 10.1016/j.biocel.2007.11.018.

23. Ahmad A. Epigenetics in Personalized Management of Lung Cancer. Adv Exp Med Biol. 2016; 890:111–22. doi: 10.1007/978-3-319-24932-2_6.

24. Steelman LS, Chappell WH, Abrams SL, Kempf RC, Long J, Laidler P, Mijatovic S, Maksimovic-Ivanic D, Stivala F, Mazzarino MC, Donia M, Fagone P, Malaponte G, et al. Roles of the Raf/MEK/ERK and PI3K/PTEN/Akt/mTOR pathways in controlling growth and sensitivity to therapy-implications for cancer and aging. Aging (Albany NY). 2011; 3:192–222. doi: 10.18632/aging.100296.

25. Kong Y, Peng Y, Liu Y, Xin H, Zhan X, Tan W. Twist1 and snail link Hedgehog signaling to tumor-initiating cell-like properties and acquired chemoresistance independently of ABC transporters. Stem Cells. 2015; 33:1063–74. doi: 10.1002/stem.1955.

26. Zhang S, Wang Y, Mao JH, Hsieh D, Kim IJ, Hu LM, Xu Z, Long H, Jablons DM, You L. Inhibition of CK2α down-regulates Hedgehog/Gli signaling leading to a reduction of a stem-like side population in human lung cancer cells. PLoS One. 2012; 7:e38996. doi: 10.1371/journal.pone.0038996.

27. Chang KJ, Yang MH, Zheng JC, Li B, Nie W. Arsenic trioxide inhibits cancer stem-like cells via down-regulation of Gli1 in lung cancer. Am J Transl Res. 2016; 8: 1133–43.

28. Lv T, Wang Q, Cromie M, Liu H, Tang S, Song Y, Gao W. Twist1-mediated 4E-BP1 regulation through mTOR in non-small cell lung cancer. Oncotarget. 2015; 6: 33006–18. doi: 10.18632/oncotarget.5026.

29. Yuan Z, Goetz JA, Singh S, Ogden SK, Petty WJ, Black CC, Memoli VA, Dmitrovsky E, Robbins DJ. Frequent requirement of hedgehog signaling in non-small cell lung carcinoma. Oncogene. 2007; 26: 1046–55. doi: 10.1038/sj.onc.1209860.

30. Clement V, Sanchez P, de Tribolet N, Radovanovic I, Ruiz i Altaba A. HEDGEHOG-GLI1 signaling regulates human glioma growth, cancer stem cell self-renewal, and tumorigenicity. Curr Biol. 2007; 17: 165–72. doi: 10.1016/j.cub.2006.11.033.

31. Bora-Singhal N, Perumal D, Nguyen J, Chellappan S. Gli1-mediated regulation of Sox2 facilitates self-renewal of stem-like cells and confers resistance to EGFR inhibitors in non-small cell lung cancer. Neoplasia. 2015; 17:538–51. doi: 10.1016/j.neo.2015.07.001.

32. Bermudez O, Hennen E, Koch I, Lindner M, Eickelberg O. Gli1 mediates lung cancer cell proliferation and Sonic Hedgehog-dependent mesenchymal cell activation. PLoS One. 2013; 8:e63226. doi: 10.1371/journal.pone.0063226.

33. Gohn CR, Blue EK, Sheehan BM, Varberg KM, Haneline LS. Mesenchyme Homeobox 2 enhances migration of endothelial colony forming cells exposed to intrauterine diabetes mellitus. J Cell Physiol. 2017; 232:1885–92. doi: 10.1002/jcp.25734.

34. Sims-Mourtada J, Izzo JG, Ajani J, Chao KS. Sonic Hedgehog promotes multiple drug resistance by regulation of drug transport. Oncogene. 2007; 26: 5674–9. doi: 10.1038/sj.onc.1210356.

35. Zhan X, Wang J, Liu Y, Peng Y, Tan W. GPCR-like signaling mediated by smoothened contributes to acquired chemoresistance through activating Gli. Mol Cancer. 2014; 13: 4. doi: 10.1186/1476-4598-13-4.

36. Park KS, Martelotto LG, Peifer M, Sos ML, Karnezis AN, Mahjoub MR, Bernard K, Conklin JF, Szczepny A, Yuan J, Guo R, Ospina B, Falzon J, et al. A crucial requirement for Hedgehog signaling in small cell lung cancer. Nat Med. 2011; 17: 1504–8. doi: 10.1038/nm.2473.

37. Tian F, Mysliwietz J, Ellwart J, Gamarra F, Huber RM, Bergner A. Effects of the Hedgehog pathway inhibitor GDC-0449 on lung cancer cell lines are mediated by side populations. Clin Exp Med. 2012; 12: 25–30. doi: 10.1007/s10238-011-0135-8.

38. Singh RR, Kunkalla K, Qu C, Schlette E, Neelapu SS, Samaniego F, Vega F. ABCG2 is a direct transcriptional target of hedgehog signaling and involved in stroma-induced drug tolerance in diffuse large B-cell lymphoma. Oncogene. 2011; 30: 4874–4886.

39. Shi Y, Fu X, Hua Y, Han Y, Lu Y, Wang J. The Side population in human lung cancer cell line NCI-H460 is enriched in stem-like cancer cells. PLoS One. 2012; 7:e33358. doi: 10.1371/journal.pone.0033358.

40. Sabo PJ, Hawrylycz M, Wallace JC, Humbert R, Yu M, Shafer A, Kawamoto J, Hall R, Mack J, Dorschner MO, McArthur M, Stamatoyannopoulos JA. Discovery of functional noncoding elements by digital analysis of chromatin structure. Proc Natl Acad Sci USA. 2004; 101:16837–42. doi: 10.1073/pnas.0407387101.

41. Liu Z, Sun M, Lu K, Liu J, Zhang M, Wu W, De W, Wang Z, Wang R. The long noncoding RNA HOTAIR contributes to cisplatin resistance of human lung adenocarcinoma cells via downregualtion of p21(WAF1/CIP1) expression. PloS One. 2013; 8: e77293. doi: 10.1371/journal.pone.0077293.

42. Liu XH, Liu ZL, Sun M, Liu J, Wang ZX, De W. The long non-coding RNA HOTAIR indicates a poor prognosis and promotes metastasis in non-small cell lung cancer. BMC Cancer. 2013; 13: 464. doi: 10.1186/1471-2407-13-464.

43. Tsai MC, Manor O, Wan Y, Mosammaparast N, Wang JK, Lan F, Shi Y, Segal E, Chang HY. Long noncoding RNA as modular scaffold of histone modification complexes. Science. 2010; 329: 689–93. doi: 10.1126/science.1192002.

44. Cao W, Ribeiro RO, Liu D, Saintigny P, Xia R, Xue Y, Lin R, Mao L, Ren H. EZH2 promotes malignant behaviors via cell cycle dysregulation and its mRNA level associates with prognosis of patient with non-small cell lung cancer. PLoS One. 2012; 7:e52984. doi: 10.1371/journal.pone.0052984.

45. Liu J, Wan L, Lu K, Sun M, Pan X, Zhang P, Lu B, Liu G, Wang Z. The Long Noncoding RNA MEG3 contributes to cisplatin resistance of human lung adenocarcinoma. PLoS One. 2015; 10:e0114586. doi: 10.1371/journal.pone.0114586.

46. Pan H, Jiang T, Cheng N, Wang Q, Ren S, Li X, Zhao C, Zhang L, Cai W, Zhou C. Long non-coding RNA BC087858 induces non-T790M mutation acquired resistance to EGFR-TKIs by activating PI3K/AKT and MEK/ERK pathways and EMT in non-small-cell lung cancer. Oncotarget. 2016; 7:49948–60. doi: 10.18632/oncotarget.10521.

47. Dong S, Qu X, Li W, Zhong X, Li P, Yang S, Chen X, Shao M, Zhang L. The long non-coding RNA, GAS5, enhances gefitinib-induced cell death in innate EGFR tyrosine kinase inhibitor-resistant lung adenocarcinoma cells with wide-type EGFR via downregulation of the IGF-1R expression. J Hematol Oncol. 2015; 8:43. doi: 10.1186/s13045-015-0140-6.

48. Lin EH, Kao YR, Lin CA, Kuo TY, Yang SP, Hsu CF, Chou TY, Ho CC, Wu CW. Hedgehog pathway maintains cell survival under stress conditions, and drives drug resistance in lung adenocarcinoma. Oncotarget. 2016; 7:24179–93. doi: 10.18632/oncotarget.8253.

49. Giroux Leprieur E, Vieira T, Antoine M, Rozensztajn N, Rabbe N, Ruppert AM, Lavole A, Cadranel J, Wislez M. Sonic Hedgehog pathway activation is associated with resistance to platinum-based chemotherapy in advanced non-small-cell lung carcinoma. Clin Lung Cancer. 2016; 17:301–8. doi: 10.1016/j.cllc.2015.12.007.

50. Huang L, Walter V, Hayes DN, Onaitis M. Hedgehog-GLI signaling inhibition suppresses tumor growth in squamous lung cancer. Clin Cancer Res. 2014; 20:1566–75. doi: 10.1158/1078-0432.CCR-13-2195.

51. Della Corte CM, Bellevicine C, Vicidomini G, Vitagliano D, Malapelle U, Accardo M, Fabozzi A, Fiorelli A, Fasano M, Papaccio F, Martinelli E, Troiani T, Troncone G, et al. SMO Gene Amplification and Activation of the Hedgehog Pathway as Novel Mechanisms of Resistance to Anti-Epidermal Growth Factor Receptor Drugs in Human Lung Cancer. Clin Cancer Res. 2015; 21:4686–97. doi: 10.1158/1078-0432.CCR-14-3319.

52. Rovida E, Stecca B. Mitogen-activated protein kinases and Hedgehog-GLI signaling in cancer: A crosstalk providing therapeutic opportunities? Semin Cancer Biol. 2015; 35:154–67. doi: 10.1016/j.semcancer.2015.08.003.

53. Shi S, Deng YZ, Zhao JS, Ji XD, Shi J, Feng YX, Li G, Li JJ, Zhu D, Koeffler HP, Zhao Y, Xie D. RACK1 promotes non-small-cell lung cancer tumorigenicity through activating Sonic Hedgehog signaling pathway. J Biol Chem. 2012; 287: 7845–58. doi: 10.1074/jbc.M111.315416.

54. Bai XY, Zhang XC, Yang SQ, An SJ, Chen ZH, Su J, Xie Z, Gou LY, Wu YL. Blockade of Hedgehog Signaling Synergistically Increases Sensitivity to Epidermal Growth Factor Receptor Tyrosine Kinase Inhibitors in Non-Small-Cell Lung Cancer Cell Lines. PLoS One. 2016; 11:e0149370. doi: 10.1371/journal.pone.0149370.

55. Arrieta O, Guzmán-de Alba E, Alba-López LF, Acosta-Espinoza A, Alatorre-Alexander J, Alexander-Meza JF, Allende-Pérez SR, Alvarado-Aguilar S, Araujo-Navarrete ME, Argote-Greene LM, Aquino-Mendoza CA, Astorga-Ramos AM, Austudillo-de la Vega H, et al. National consensus of diagnosis and treatment of non-small cell lung cancer. [Article in Spanish]. Rev Investig Clin. 2013 (Suppl); 65:S5–84.

56. Arrieta O, Cardona AF, Federico Bramuglia G, Gallo A, Campos-Parra AD, Serrano S, Castro M, Avilés A, Amorin E, Kirchuk R, Cuello M, Borbolla J, Riemersma O, et al, and CLICaP. Genotyping non-small cell lung cancer (NSCLC) in Latin America. J Thorac Oncol. 2011; 6:1955–9. doi: 10.1097/JTO.0b013e31822f655f.

57. Guzmán C, Bagga M, Kaur A, Westermarck J, Abankwa D. ColonyArea: an ImageJ plugin to automatically quantify colony formation in clonogenic assays. PLoS One. 2014; 9:e92444. doi: 10.1371/journal.pone.0092444.