Aberrant DNA methylation in non-small cell lung cancer-associated fibroblasts.

Epigenetic changes through altered DNA methylation have been implicated in critical aspects of tumor progression, and have been extensively studied in a variety of cancer types. In contrast, our current knowledge of the aberrant genomic DNA methylation in tumor-associated fibroblasts (TAFs) or other stromal cells that act as critical coconspirators of tumor progression is very scarce. To address this gap of knowledge, we conducted genome-wide DNA methylation profiling on lung TAFs and paired control fibroblasts (CFs) from non-small cell lung cancer patients using the HumanMethylation450 microarray. We found widespread DNA hypomethylation concomitant with focal gain of DNA methylation in TAFs compared to CFs. The aberrant DNA methylation landscape of TAFs had a global impact on gene expression and a selective impact on the TGF-β pathway. The latter included promoter hypermethylation-associated SMAD3 silencing, which was associated with hyperresponsiveness to exogenous TGF-β1 in terms of contractility and extracellular matrix deposition. In turn, activation of CFs with exogenous TGF-β1 partially mimicked the epigenetic alterations observed in TAFs, suggesting that TGF-β1 may be necessary but not sufficient to elicit such alterations. Moreover, integrated pathway-enrichment analyses of the DNA methylation alterations revealed that a fraction of TAFs may be bone marrow-derived fibrocytes. Finally, survival analyses using DNA methylation and gene expression datasets identified aberrant DNA methylation on the EDARADD promoter sequence as a prognostic factor in non-small cell lung cancer patients. Our findings shed light on the unique origin and molecular alterations underlying the aberrant phenotype of lung TAFs, and identify a stromal biomarker with potential clinical relevance.

Introduction

Lung cancer is the leading cause of cancer-related deaths worldwide. Non-small cell lung cancer (NSCLC) accounts for ~85% of all lung cancers, and includes two major histologic subtypes: adenocarcinoma (ADC) and squamous cell carcinoma (SCC) (1). Although both NSCLC subtypes are epithelial in origin, there is growing awareness that tumor progression in NSCLC and other solid tumors is driven by the aberrant coevolution of carcinoma cells and surrounding stromal cells (2,3). Among the latter, tumor-associated fibroblasts (TAFs) are the most abundant cell type, and have been implicated in all major steps of tumor progression including cancer cell growth, invasion, chemorresistance and stemness (4). Importantly, there is growing evidence that TAFs exhibit enhanced tumor-promoting effects compared to fibroblasts from unaffected tissue (5). Thus, a better understanding of the aberrant molecular differences between normal fibroblasts and TAFs is needed to unveil their tumor-promoting effects.

Most of our knowledge of the tumor-promoting effects of lung TAFs has been obtained from cell culture assays and animal models (5–7). These studies have consistently reported that the aberrant phenotype of TAFs is maintained for some passages in culture in the absence of continuous interaction with carcinoma cells. Similar observations have been reported in other cancer types, strongly supporting that critical phenotypic alterations in TAFs are maintained through epigenetic mechanisms (3,8,9).

DNA methylation is the most well-studied epigenetic alteration in cancer, owing in part to recent developments in genome-wide DNA methylation profiling techniques (10). DNA methylation involves the covalent modification of the cytosine in a cytosine-phosphate-guanine (CpG) island within genomic DNA, which is catalyzed by DNA methyltransferases. Previous studies have shown that global loss of DNA methylation (hypomethylation) and promoter hypermethylation-associated gene inactivation are common epigenetic hallmarks of cancer cells (11,12). In NSCLC, several DNA methylation alterations have been described in association with the neoplastic transformation, and some of them have been pointed as potential biomarkers with clinical relevance for diagnosis, prognosis and response to therapy (10,13). However, former DNA methylation studies in lung cancer examining either whole tumor tissue samples or cancer cell lines have omitted the epigenetic alterations specifically affecting TAFs or other stromal components (13,14). Indeed, our current knowledge of genome-wide epigenetic alterations within the tumor stroma in human cancer is very scarce, and has been only explored in breast and gastric cancers (8,15).

Here, we present a detailed analysis of the DNA methylation patterns of low-passage primary cultures of TAFs from 12 surgical patients diagnosed with early stage NSCLC, and paired control fibroblasts (CFs) isolated from unaffected lung parenchyma. We found widespread DNA hypomethylation concomitant with focal gain of DNA methylation in TAFs compared to CFs. These epigenetic changes had a global impact on gene expression and, remarkably, a selective impact in the promoters of critical transcription factors of the TGF-β pathway, including SMAD3, which was associated with an aberrant response to exogenous TGF-β1. Pathway enrichment analysis of the aberrant genomic methylation in TAFs provided new insights on their partial bone marrow origin. Moreover, we found that aberrant DNA methylation of selected candidates was retrospectively associated with shorter survival on NSCLC patients, thereby uncovering DNA methylation biomarkers with potential clinical value.

Materials and methods

Tissue samples and primary human lung fibroblasts

Lung tissue samples were obtained from a cohort of 20 early stage NSCLC surgical patients (10 ADC, 10 SCC) prior to 2013 at the Hospital Clinic de Barcelona (HCB, Spain) with the approval of the Ethics Committees of the HCB and the Universitat de Barcelona. All patients gave their informed consent. Selected patients were male, chemo-naïve, Caucasian, ≥55 years old and current smokers with confirmed ADC or SCC diagnosis (further clinical characteristics are shown in Supplementary Table 1, available at Carcinogenesis Online). Samples from tumor and paired tumor-free lung parenchyma were collected. A fraction of each sample was paraffin-embedded for histologic analysis, whereas the remaining was used to isolate CFs and TAFs by outgrowth of tissue explants as described elsewhere (6). In brief, tumor samples were chopped in ~1mm2 fragments, gently distributed into a six-well plate, and incubated in regular culture medium containing 10% fetal bovine serum (FBS) and antibiotics for 3–4 weeks, changing the medium every 2–3 days. During this time frame, fibroblasts outgrew from tissue fragments and proliferated over the well surface. Tissue explants were removed and fibroblasts were subcultured by trypsinization before reaching confluence. Fibroblasts were used up to passage 6. Primary cultures were tested negative for mycoplasma. The mesenchymal origin of the fibroblasts was confirmed by their positive and negative immunofluorescence staining with vimentin and pan-cytokeratin antibodies, respectively. DNA methylation profiling of primary fibroblasts was conducted on 12 patients, whereas all other molecular biology and microscopy analyses were carried out on randomly selected patients from our cohort.

Histology

Immunohistochemical analyses of α-SMA, CD34, CD45 and EDARADD were performed on tissues from our cohort (n = 20) using the Bond automated immunohistochemistry system (Leica Microsystems). Nuclei were counterstained with hematoxylin. Primary antibodies included anti-α-SMA (1A4), anti-CD34 (QBEnd-10), anti-CD45 (2B11+, Dako) and anti-EDARADD (ab121581, Abcam). Histological staining were imaged with a bright-field microscope (BX43) coupled to a digital camera (DP72) using a ×40 objective (Olympus). All image processing henceforth was carried out with Image J (16) under the guidance of our pathologist (JR). To assess the presence of fibrocyte-like cells in CD34 and CD45 staining, fibroblasts were identified within stromal mesenchymal cells according to their elongated spindle-shaped nuclei, discarding those lining blood vessels (6). EDARADD staining in fibroblasts was scored in a semiquantitative manner as described elsewhere (6). In brief, EDARADD staining was scored blind by two independent observers according to three categories (negative, weak and strong) (see Supplementary Material, available at Carcinogenesis Online for further details).

Cell culture

Primary fibroblasts were maintained in culture medium as reported elsewhere (6). Unless otherwise indicated, experiments were conducted on tissue culture plastic substrata coated with 0.1mg/ml collagen-I solution (Millipore) overnight at 4°C. Fibroblasts were seeded at 8×103 cells/cm2. For methylation array studies, fibroblasts were cultured for 5 days in serum-free culture medium (SFM) supplemented with 0.1% FBS. In some experiments, fibroblasts were cultured on collagen-coated polyacrylamide gels engineered to exhibit normal (~1 kPa) or tumor-like (~30 kPa) Young′s elastic moduli (E) (6) in the presence or absence of 2.5ng/ml TGF-β1 (R&D Systems) for 5 days. For extracellular matrix (ECM) expression analysis, cells were cultured on collagen-coated rigid substrata in SFM supplemented with 2.5ng/ml TGF-β1 for 5 days.

Infinium 450K DNA methylation profiling

We used the EZ DNA Methylation Kit (Zymo Research) for bisulfite conversion of 500ng genomic DNA of CFs and paired TAFs from 12 randomly selected patients from of our cohort (clinical characteristics shown in Supplementary Table 1, available at Carcinogenesis Online). The corresponding DNA methylation profiles were obtained with the Infinium 450K Methylation Array as described previously (17), which quantifies methylation levels (β-value) of ~450000 CpGs located both at gene promoter and non-promoter regions (10). Raw fluorescence intensity values were normalized in Illumina Genome Studio software (V2010.3) using ‘control normalization’ without background correction. Normalized intensities were used to calculate β-values (GSE68851). All methylation data analysis was carried out henceforth with the R Software for Statistical Computing (v3.1.1). Data points with insufficient fluorescent intensities (P > 0.01) were excluded from the analysis. Likewise, genotyping probes present on the chip as well as DNA methylation probes overlapping with known single-nucleotide polymorphisms were removed. A differential methylation analysis between CFs and TAFs was conducted applying paired comparisons via moderated t-statistics provided by the linear models implemented in limma Bioconductor package (18) to identify statistically significant differential DNA methylation differences. Clustering analysis was applied to visualize the differential DNA methylation patterns between groups. Further details are given in Supplementary Material, available at Carcinogenesis Online.

Gene enrichment analysis

A global gene enrichment analysis (GEA) was applied to the statistically significant differentially methylated CpG sites (referred to as differential CpGs) between CFs and TAFs using R software. All possible genes related to each differential CpG were individually considered, sharing the same β-value. The resulting list of genes was first filtered to avoid redundances, assigning the maximum β-value observed among repeated genes; secondly, by selecting genes with absolute │Δβ│ = │βTAF − βCF│> 0.2. The package clusterProfiler (v2.0.0) (19) was used to compute GEA over the final list of genes. Queried biological pathways were from KEGG (20) and Reactome (21) databases, accessed through the packages KEGG.db and reactome.db, respectively. In both cases, a P value for each gene was calculated based on a hypergeometric distribution test (5% FDR). The package topGO was used to conduct GEA for the Gene Ontology (GO) Biological Processes database (22). Results were prioritized based on P values obtained with Fisher’s exact test statistic and weight method (23). GEA corresponding to the TGF-β pathway was conducted with the KEGG TGF-β signalling (hsa04350) through the package KEGGgraph v1.24.0 (24). Genes annotated in this pathway were matched to the former list of genes with │Δβ│>0.2. The same procedure was applied to the TGF-β receptor signaling pathway from Netpath database (NetPath_7) (25). Complete package references and further details are provided in Supplementary Material, available at Carcinogenesis Online.

Pyrosequencing

Bisulfite-treated DNA was used as a template for PCR. The primers for PCR amplification and sequencing were designed with PyroMark assay design software version 2.0.01.15. Primer sequences (Supplementary Material, available at Carcinogenesis Online) were designed, when possible, to hybridize with CpG-free sites to ensure methylation-independent amplification. Pyrosequencing analyses were performed as described previously (26).

qRT-PCR

For transcriptional analysis of selected genes that exhibited differential DNA methylation between CFs and TAFs, cells were cultured using the same protocol than for DNA methylation profiling. Total RNA was isolated using the RNeasy Mini kit (QIAGEN) and reverse-transcribed into cDNA using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems). qRT-PCR reactions were performed in triplicates on a 7900HT Fast Real-Time PCR system (Applied Biosystems) using 20ng cDNA, SYBR Green PCR Master Mix (Applied Biosystems) and specific primers for SMAD3, EDARADD, CHI3L1 and ACTB (used as housekeeping gene). Primer pairs were designed with PerlPrimer v1.1.14 software and validated by gel electrophoresis to amplify specific single products. For ECM expression analysis, RNA extraction and reverse transcription was conducted using the same procedure. Real-time PCR reactions were performed on 50ng of each cDNA sample using TaqMan Gene Expression Master Mix and TaqMan gene-specific primer pairs and probes for COL1A1, LOX, SPARC, TNC and POLR2A (used as a housekeeping gene). Primers and probes for detection of EDA-FN were customary designed based on the sequences reported elsewhere (27). Relative gene expression with respect to a housekeeping gene was assessed as 2−ΔCt as described previously (28). Details of primers and probes are given in Supplementary Material, available at Carcinogenesis Online.

Immunofluorescence

Fibroblast cultures were fixed with 4% paraformaldehyde, permeabilized and blocked with 0.2% Triton X-100, 1% bovine serum albumin (Sigma) and 6% FBS (Gibco), firstly incubated with an anti-α-SMA mouse antibody (clone 1A4, Sigma), and secondly with a Cy3 goat anti-mouse IgG secondary antibody (Jackson). Nuclei were counterstained with Hoechst 33342 (Molecular Probes). Fluorescence images were acquired with an Eclipse TE2000 microscope (Nikon) at nine randomized locations with an EM-CCD C9100a camera (Hamamatsu) using Metamorph software (Molecular Devices) and a ×20 objective. Each image was background corrected and its total intensity (I) and cell number (N) was measured with Image J. All I versus N data for each culture condition were least-squares fitted to a linear function with MATLAB (Mathworks) to assess the α-SMA intensity per cell as the fitted slope.

Flow cytometry

The percentages of cells positive for the cell surface markers CD34 and CD45 were assessed by flow cytometry. For this purpose, primary fibroblasts (CFs and paired TAFs) were cultured in 10% FBS culture medium for 24h before the experiment, and detached with trypsin-EDTA (Sigma) for 1min at 37°C. Fresh 10% FBS medium was added to quench the trypsin and cells were pelleted by centrifugation and kept with 10% FBS in suspension at room temperature for 2h to enable the restoration of cell surface epitopes. Suspended cells were washed with PBS, incubated for 30min with 2 µg/ml APC-conjugated antibodies against either CD34 or CD45 (Biolegend) in blocking solution (5% FBS in PBS) on ice in dark conditions. Cells were washed with 0.5% bovine serum albumin in PBS solution, resuspended with PBS and analyzed as single cells with a flow cytometer (Gallios, Beckman Coulter) using FlowJo 10.0 software software. Additional details are provided in Supplementary Material, available at Carcinogenesis Online.

Traction force microscopy

Maps of traction forces were assessed in single fibroblasts using constrained Fourier transform traction microscopy as described elsewhere (29,30). Briefly, collagen-coated polyacrylamide gels with embedded fluorescent nanobeads were prepared by mixing 7.5% acrylamide and 0.06% bisacrylamide to achieve E ~ 6 kPa. Cells were cultured in SFM in the absence or presence of 2.5ng/ml TGF-β1 for 5 days. A map of gel displacement and the corresponding traction force map were computed as described previously (29,30). The average traction force per unit area was computed and averaged over ~10 cells per condition.

Re-analysis of public cancer DNA methylation and expression datasets

Available Infinium 450K DNA methylation experimental data of clinical NSCLC specimens with detailed clinical annotation (13) (GSE39279) were correlated with survival parameters. For each selected gene, we collected its differentially methylated CpGs and established a β-value threshold to split the data in two groups: high (β > 0.33) and low (β < 0.33). The corresponding Kaplan–Meier survival curves were obtained with R software, and their statistical significance was assessed with the log-rank test. A similar procedure was applied on gene expression data and survival parameters from stages I to II lung NSCLC available elsewhere (31) (GSE31210).

Statistical analysis

Two-group comparisons of non-methylation array data were performed with the Student’s t-test unless otherwise indicated. Statistical significance was assumed at P < 0.05.

Results

DNA methylation profiling reveals a global hypomethylation in TAFs

A major hallmark of lung TAFs in vivo is the expression of alpha-smooth muscle actin (α-SMA) (6), which is indicative of an activated/myofibroblast-like phenotype (32). Primary TAFs recapitulate this hallmark in culture, according to their larger α-SMA expression compared to CFs in all patients examined (Figure 1A). To shed light on the epigenetic alterations underlying the persistent TAF activation in culture, we conducted a genome-wide DNA methylation profiling on TAFs and paired CFs from 12 randomly selected patients from our cohort using the Infinium 450K Methylation Array. Comparing β-values in CFs and TAFs identified 18520 statistically significant differentially methylated CpG sites outside X chromosome (Δβ ≠ 0; P < 0.002, 5% FDR) (Supplementary Table 2, available at Carcinogenesis Online). Differentially methylated sites in CFs and TAFs were preferentially located in non-promoter (76%) rather than gene promoter (24%) sequences. In contrast, no statistically significant differentially methylated sites were found within TAFs from different histologic subtypes. Therefore, only two groups (CFs and TAFs) were considered in further analyses.

Figure 1.

Primary lung TAFs exhibit global DNA hypomethylation and focal gain of DNA methylation. (A) Representative fluorescence images illustrating α-SMA overexpression in lung TAFs compared to paired CFs obtained with a ×20 objective (top). Scale bars here and thereafter, 30 μm. The bottom plot shows the quantification of fold α-SMA intensity per cell of fibroblasts from four randomized patients. (B) Unsupervised clustering of 1452 CpG sites with marked differential methylation in TAFs and CFs from 12 randomized patients of our cohort and (C) normalized distribution (relative density) of the corresponding β-values. Dashed vertical lines indicate median β-values. *P < 0.05; **P < 0.01; ***P < 0.001 (here and thereafter).

Among the list of 18520 differential CpGs in CFs and TAFs, 1452 exhibited marked DNA methylation differences (taken as absolute Δβ ≥ 0.2 henceforth), which corresponded to 750 distinct genes. An unsupervised clustering analysis of the latter 1452 probes shown in Figure 1B revealed a clear sorting of CFs and TAFs, with a global DNA hypomethylation concomitantly with DNA hypermethylation of a smaller set of genomic regions in TAFs compared to CFs. To analyze quantitatively the global loss of DNA methylation, we assessed the distribution of β-values within the 1452 probes with marked differential methylation. A clear shift towards lower β-values was observed in TAFs compared to CFs (Wilcoxon test, P < 0.001), with a 19.4% reduction of median DNA methylation values in TAFs (dashed vertical lines in Figure 1C).

Global hypomethylation in TAFs is partially elicited by TGF-β1 but not matrix stiffening

To examine the relationship between an activated phenotype and global loss of DNA methylation, we analyzed the impact of two fibroblast activation signals frequently observed in the desmoplastic tumor stroma (9): TGF-β1 and ECM stiffening. For this purpose, we first cultured CFs from two randomly selected patients (P5 and P28) in stiff substrata for 5 days in the presence or absence of TGF-β1 and analyzed their corresponding DNA methylation profiles with the 450K Methylation Array. TGF-β1 increased α-SMA expression in CFs (Figure 2A) and elicited a statistically significant reduction on DNA methylation (2.2%, P < 0.01) (Figure 2B). However, the latter drop could not account for the 18.8% reduction of DNA methylation observed in TAFs from the same patients in stiff substrata in the absence of TGF-β1 (P < 0.001, Figure 2C), indicating that such 18.8% reduction is largely imputable to the transformation process. Moreover, we found that 14.3% of the CpG sites in CFs with the largest variability upon TGF-β1 treatment (Supplementary Table 3A and Supplementary Material, available at Carcinogenesis Online) overlapped the 1452 differentially methylated CpG sites in CFs and TAFs. Secondly, we cultured CFs on polyacrylamide gels engineered to exhibit either normal- (~1 kPa) or tumor-like (~30 kPa) rigidities in the presence of TGF-β1 (Figure 2D). Unlike Figure 2B, DNA methylation differences of CFs cultured in these conditions were very modest, and only showed 0.7% reduction in the median β-value in cells cultured in the stiffest gels (non-statistically significant, Figure 2E). Likewise, only one of the differentially methylated CpG sites with the largest variability induced by matrix stiffening (Supplementary Table 3B and Supplementary Material, available at Carcinogenesis Online) overlapped with the 1452 differentially methylated CpG sites in CFs and TAFs. Altogether, these results indicate that TGF-β1 stimulation of CFs in a stiff microenvironment contributes to the aberrant DNA methylation pattern observed in TAFs but does not fully recapitulate it, whereas matrix stiffening in a TGF-β1-rich microenvironment does not have any significant contribution.

Figure 2.

Activating CFs with TGF-β1 partially mimics the genomic methylation changes in TAFs. (A) Effect of TGF-β1 on α-SMA fluorescence staining in CFs. (B) Effect of TGF-β1 on the DNA methylation distribution of the list of 1452 differential CpG sites in CFs from two randomly selected patients (P5 and P28). Dashed horizontal lines indicate the median of each distribution. (C) DNA methylation distribution of the list of 1452 differential CpG sites in CFs and paired TAFs from the same patients. (D) Outline of the culture assay based on polyacrylamide gels with normal- or tumor-like rigidities in the presence of TGF-β1. (E) Effect of culture conditions shown in (D) on the DNA methylation distribution of the list of 1452 differential CpG sites in CFs from the same patients.

DNA methylation changes affect pathways associated with ECM/focal adhesions and the FCγ receptor

To gain insights on the functional consequences caused by the observed changes on DNA methylation, we conducted a pathway enrichment analysis selecting those genes from our list with marked methylation differences using three complementary databases: KEGG, Reactome and GO Biological Processes. Figure 3A shows a plot with the overrepresented KEGG pathways (P < 0.05, 5% FDR) (Supplementary Table 4A, available at Carcinogenesis Online). The corresponding enrichments obtained with Reactome and GO are shown in Supplementary Figures 1–2 and Supplementary Tables 4B–E, available at Carcinogenesis Online. Among the statistically significant enriched pathways or processes, only two were coincidentally reported in all three databases. Such pathways were closely related to the ECM/focal adhesions and the FCγ receptor (FCγR). Enhanced ECM deposition and cell–ECM interactions have been described previously as myofibroblast hallmarks (4,9). In contrast, FCγRs have been associated with immune cells rather than myofibroblasts (33,34). A straightforward explanation for the FCγR pathway enrichment is that it reflects the presence of fibrocytes, since these are bone marrow-derived cells that express FCγRs as well as leukocyte (CD45) and hematopoietic progenitor (CD34) markers, and can differentiate into fibroblasts/myofibroblasts in response to TGF-β (33,34). To test this possibility, we examined CD34 and CD45 expression in tumor histologic sections from our cohort (n = 20) upon the guidance of our pathologist (JR), and observed several CD34+ and CD45+ stromal mesenchymal cells (Figure 3B and Supplementary Figure 3, available at Carcinogenesis Online). In addition, we assessed CD34+ and CD45+ fibroblasts in culture by flow cytometry. Although the percentages of CD34+ and CD45+ fibroblasts were collectively modest, they were consistently larger in TAFs compared to paired CFs in 3 randomly selected patients (Figure 3C,D and Supplementary Figure 3, available at Carcinogenesis Online). These results strongly support that a fraction of primary lung TAFs may be fibrocytes or fibrocyte-like cells in origin.

Figure 3.

Pathway-enrichment analysis reveals that a fraction of TAFs are fibrocyte or fibrocyte-like cells in origin. (A) Statistically significant overrepresented KEGG pathways within the 750 distinct genes corresponding to the list of 1452 CpG sites with marked differential methylation between TAFs and CFs. Pathway’s circle size is proportional to the number of annotated genes. Genes annotated to each pathway are color-coded according to their fold DNA methylation. None, one or two circles around each gene indicate that differential methylation was found within non-promoter, promoter, or both non-promoter and promoter regions, respectively. (B) Representative images of tumor histologic sections of a randomly selected patient from our cohort stained for fibrocyte markers CD34 (left) and CD45 (right). Black arrows point to CD34+ and CD45+ (non-endothelial) spindle-shaped stromal mesenchymal cells. (C) Histograms of CD34+ (left) and CD45+ (right) cultured CFs and TAFs from the same patient as in (B) assessed by flow cytometry. (D) Fold (TAFs/CFs) percentages of CD34+ and CD45+ fibroblasts from three randomly selected patients assessed as in (C).

DNA methylation alterations of genes and miRNAs involved in fibroblast activation

TGF-β signaling is essential for fibroblast activation, and its corresponding KEGG pathway was found overrepresented in TAFs (Figure 3A). To delimit the scope of altered TGF-β signaling in TAFs, we sought to identify those genes annotated to the TGF-β signaling pathway available in either KEGG or NetPath databases that appeared in our list of differentially methylated genes in TAFs with respect to CFs. Among the 274 genes annotated to the TGF-β pathway in both databases, a marked differential methylation was observed selectively in 24 genes (9%) (Supplementary Table 5, available at Carcinogenesis Online). Of note, only 8 of the latter 24 genes exhibited a marked differential DNA methylation in promoter regions, including 4 important transcription factors: SMAD3 and the runt domain-containing family RUNX1, RUNX2 and RUNX3. Among them, SMAD3 was hypermethylated whereas RUNX genes exhibited a loss of DNA methylation in TAFs. Thus, even though methylation changes were only found in a modest percentage of genes annotated to the TGF-β pathway, these changes affected the critical component SMAD3.

In addition to TGF-β signaling, there is evidence that fibroblast activation may be epigenetically regulated by a growing list of miRNAs (9,35). Among those described previously in the literature, none exhibited a robust differential methylation in promoter regions in our dataset. In contrast, we identified three miRNAs with a marked hypomethylation affecting different CpGs in promoter regions that have not been previously associated with fibroblast activation: miR-296, miR-298 and miR-1249.

Impact of aberrant DNA methylation in transcription

Aiming to examine the impact of DNA methylation alterations on gene transcription, we analyzed the differential gene expression of genes undergoing large DNA methylation changes at promoter sequences, which are sorted in Tables 1 and 2 according to their number of differentially methylated CpG sites and maximum DNA methylation change (full list is given in Supplementary Table 6, available at Carcinogenesis Online). Close examination of the top genes in Tables 1 and 2 identified SMAD3 and EDARADD among those with the most consistent DNA methylation changes in the lists of promoter hypermethylated and hypomethylated genes, respectively. To expand our candidate selection, we took advantage of the differential transcriptional profiling of lung TAFs and paired CFs obtained from a cohort of 15 early stage NSCLC patients analyzed by Navab and coworkers (5). This study identified 46 differentially expressed genes in TAFs and CFs with an absolute fold change >2. Among them, 22 (48%) were coincidental with the probes differentially methylated in our dataset (Supplementary Table 7, available at Carcinogenesis Online), yet only CHI3L1 consistently exhibited a marked DNA methylation change (│Δβ│ > 0.2) in its promoter region. Based on these analyses, SMAD3, EDARADD and CHI3L1 were selected for subsequent validation by pyrosequencing. Our results confirmed the increased DNA methylation in SMAD3 promoter sequence (P < 0.001). Likewise, DNA hypomethylation was validated in both EDARADD and CHI3L1 promoter sequences (P < 0.001) in TAFs compared to paired CFs (Figure 4A). To determine the impact of these epigenetic changes on transcriptional activity, we assessed mRNA levels of the three selected genes by qRT-PCR (Figure 4B). SMAD3 levels were down-regulated in TAFs compared to paired CFs in all patients (P < 0.05), in agreement with the transcriptional repression associated with promoter DNA hypermethylation. Conversely, increased expression levels were observed in both EDARADD (P < 0.05) and CHI3L1 (P < 0.05) in TAFs with respect to paired CFs, in agreement with transcriptional reactivation upon loss of DNA methylation at the promoter.

Table 1. Differential methylation in promoter sequences in lung TAFs. Hypermethylated genes with marked differential methylation (│Δβ │ > 0.2) in promoter sequences in lung TAFs in ≥ 2 CpG sites

	Gene symbol	CpG sites	Highest Δβ = βTAF − βCF
Hypermethylated	SMAD3	4	0.236
	SYNPO	3	0.211
	GPR88	2	0.276
	C7orf54	2	0.272
	SND1	2	0.272
	TMEM212	2	0.248
	LOC404266	2	0.241
	TTC39C	2	0.230
	ZMIZ1	2	0.225
	EYA4	2	0.220

Table 2.

Hypomethylated genes with marked differential methylation (│Δβ │ > 0.2) in promoter sequences in lung TAFs in ≥ 2 CpG sites

	Gene symbol	CpG sites	Highest Δβ = βTAF − βCF
Hypomethylated	IVL	5	−0.261
	RUNX1	4	−0.319
	C22orf9	4	−0.298
	MIR1249	4	−0.298
	NTM	4	−0.245
	CSGALNACT1	3	−0.314
	IPO5	3	−0.303
	EDARADD	3	−0.295
	SLAMF8	3	−0.283
	SLC22A18AS	3	−0.271
	SLC22A18	3	−0.271
	SMCP	3	−0.244
	MIR298	3	−0.230
	MIR296	3	−0.230
	SCT	2	−0.351
	PARP4	2	−0.316
	EPS15	2	−0.297
	WIPF1	2	−0.291
	HRH1	2	−0.280
	FAM49A	2	−0.277
	CHRM5	2	−0.270
	RUNX3	2	−0.263
	TSPAN9	2	−0.262
	CARD14	2	−0.259
	S100A3	2	−0.258
	GNASAS	2	−0.254
	ZC3H12D	2	−0.254
	PLEKHA5	2	−0.249
	LEPR	2	−0.242
	LEPROT	2	−0.242
	CHRNA1	2	−0.231
	ESRRG	2	−0.229
	STRA6	2	−0.224
	HTR1D	2	−0.221
	GSTA3	2	−0.220
	SH3BP4	2	−0.217
	LCE1B	2	−0.207

Figure 4.

Validation of differential methylated promoters and the prognostic value of EDARADD in NSCLC with high stromal mesenchymal expression in-vivo. (A) DNA methylation of selected genes measured by pyrosequencing in TAFs and paired CFs from 12 patients. (B) Fold (TAFs/CFs) relative mRNA expression of selected genes in five randomized patients by qRT-PCR. Horizontal dashed lines here and thereafter are added as a reference. (C) Kaplan–Meier survival curve using publicly available data of a cohort of 204 NSCLC patients (161 ADC, 43 SCC) sorted by their EDARADD DNA methylation status assessed in primary tumoral DNA (13). (D) Kaplan–Meier survival curve using publicly available data of a cohort of 226 NSCLC patients (226 ADC) sorted by their EDARADD expression assessed in primary tumors (31). (E) Representative images of histologic sections from our cohort stained for EDARADD obtained from unaffected lung parenchyma (control, left) or tumor (right). Black arrows within inserts point to EDARADD+ spindle-shaped mesenchymal cells. (F) Scoring of EDARADD protein expression in (non-endothelial) spindle-shaped stromal cells in our cohort of 20 NSCLC patients (10 ADC, 10 SCC).

Clinical value of selected genes

The clinical impact of the three validated genes was assessed in terms of prognosis. For this purpose, we first examined the DNA methylation status of the three genes on a previously published cohort of 204 tumor tissue samples from NSCLC patients (13). Of note, EDARADD promoter hypomethylation was associated with a shorter disease-associated survival (Log-rank test, P < 0.01) in comparison with those patients exhibiting higher levels of promoter methylation (Figure 4C). In agreement with the worst prognosis conferred by EDARADD promoter DNA hypomethylation, high levels of EDARADD expression were consistently associated with shorter survival (Log-rank test, P < 0.01) (Figure 4D) according to a publicly available gene expression database from a cohort of 226 NSCLC patients (31). Unlike EDARADD, no prognostic information was observed in SMAD3 or CHI3L1 methylation status in the same data sets (Supplementary Figure 4, available at Carcinogenesis Online). To further validate the expression of EDARADD in clinical specimens, we examined its protein level in TAFs by immunohistochemistry in our cohort of 20 patients (Figure 4E and Supplementary Figure 5, available at Carcinogenesis Online). In tumor samples, EDARADD expression was heterogeneous and could be observed in both TAFs and cancer cells, whereas it was largely absent in parenchymal cells from unaffected tissue. In agreement with our qRT-PCR data, quantitative image analyses revealed a higher percentage of stromal fibroblasts with strong EDARADD staining compared to unaffected lung parenchyma (used as control) (Figure 4F; P < 0.05).

SMAD3 hypermethylation in TAFs is associated with an aberrant response to TGF-β1 in terms of contractility and ECM expression

To get insights on the functional impact of the selective SMAD3 epigenetic silencing in TAFs, we took advantage of a previous work on Smad3-null mice, which reported two skin fibroblast-associated altered functions in vivo: (i) increased ECM deposition in response to exogenous TGF-β1; and (ii) enhanced wound healing rate possibly due to enhanced contractility upon skin injury (36). To assess whether the latter in vivo phenotype of skin fibroblasts could be extended to lung TAFs in culture, we first examined the expression of a panel of wound-related ECM genes (COL1A1, EDA-FN, LOX, SPARC and TNC) commonly associated with fibroblast activation (32) by qRT-PCR. For this purpose, TAFs and paired CFs from randomly selected patients were cultured with TGF-β1 for 5 days. Fold expression data (TAFs/CFs) showed a marked (>>1-fold) overexpression of all ECM genes except TNC in most patients (Figure 5A). Secondly, we assessed fibroblast contractility by traction force microscopy. Mean traction forces increased in response to TGF-β1 in both TAFs and CFs (Figure 5B); however, such contractility increase was ~40% larger in average in TAFs than in CFs (Figure 5C). To our knowledge, these findings unveil for the first time an aberrant hyperresponsiveness of lung TAFs to TGF-β1 in terms of ECM expression and contractility.

Figure 5.

SMAD3 epigenetic silencing in TAFs correlates with enhanced response to TGF-β1, ECM expression and contractility. (A) Fold (TAFs/CFs) relative mRNA expression of a panel of wound-related ECM genes in response to TGF-β1 in four randomized patients. (B) Illustrative traction maps of a single CF and paired TAF from the same patient under basal conditions or TGF-β1 stimulation. (C) Fold (TAFs/CFs) average traction force in four randomized patients. All data are mean ± SE.

Discussion

In the present study, we interrogated for the first time the DNA methylation landscape changes between CFs and paired TAFs from surgical patients with early stage NSCLC using the high-resolution Infinium 450K DNA Methylation Array. We observed a global DNA hypomethylation and selective promoter DNA hypermethylation in TAFs compared to CFs. Therefore, the DNA methylation alterations in TAFs could be attributed to their neoplastic transformation, and included a ~19.4% reduction in the median DNA methylation levels compared to CFs. Such loss of DNA methylation is in quantitative agreement with that reported in TAFs from two gastric cancer patients (15), and in qualitative agreement with the global DNA hypomethylation reported in TAFs from breast cancer (8) and from a mouse pancreatic cancer model (37). Likewise, we confirmed and expanded the catalogue of DNA methylation alterations in TAFs (15). Collectively, these findings support that global DNA methylation loss concurrent with focal hypermethylation define a general epigenetic hallmark of TAFs from solid tumors, and may underlie their ‘phenotypic memory’ in culture.

It is commonly assumed that lung TAFs are a heterogenous population with similar tumor-promoting effects regardless their histologic subtype of origin. However, some of us recently challenged this assumption by reporting phenotypic differences in TAFs from major histologic subtypes (6). In contrast, we did not find statistically significant diferentially methylated CpGs within subtypes. A possible interpretation of these results is that a larger patient cohort may be required to unveil the fine subtype-specific DNA methylation alterations. Alternatively, it is conceivable that, in addition to DNA methylation, the ‘epigenetic memory’ of the subtype-specific phenotypic differences in TAFs is supported by other epigenetic alterations.

The robust global DNA hypomethylation in TAFs raises the question of its underlying mechanism(s). A potential straightforward molecular mechanism could be an altered DNA methyltransferase activity. However, this mechanism is unlikely according to previous data (15). Alternatively, it is conceivable that chronic fibroblast activation and DNA hypomethylation are mechanistically interconnected (38–40). To begin to explore this possibility in primary lung TAFs, we examined the genomic methylation effects of two potent fibroblast activator cues: TGF-β1 and matrix stiffening (6,9,41). Interestingly, whereas TGF-β1 treatment recapitulated some of the epigenetic changes associated with tumorigenesis, matrix stiffening failed to elicit similar alterations. These results are in agreement with a previous study describing that TGF-β1 might be associated with global DNA hypomethylation in gastric cancer (42). Importantly, our data suggest that TGF-β1 may be necessary but not sufficient to induce the aberrant genomic methylation in TAFs, suggesting that factors other than TGF-β1 reshape their epigenome. In support of the latter interpretation, it has been suggested that chronic inflammation, aberrant differentiation of local or recruited circulating progenitor cells like fibrocytes, or abnormal paracrine signaling by cancer cells may also underlie the global demethylation of TAFs (3).

The elucidation of the pathologic consequences of global hypomethylation remains a matter of intense research. In cancer cells, global loss of DNA methylation has been associated with genomic alterations including chromosome instability, activation of transposable elements and loss of genomic imprinting, which may contribute to cancer development through either increased genomic instability and/or aberrant gene expression (43,44). Currently available data do not support the former events in TAFs (15). In contrast, two lines of evidence support that the aberrant TAF methylome may have a global impact on gene expression. First, within a list of 46 genes with marked differential expression between lung TAFs and CFs reported by Navab and coworkers (5), 22 exhibited also differential DNA methylation in our study (Supplementary Table 7, available at Carcinogenesis Online). Such large overlap is in marked contrast with the single differentially methylated gene from the 46 gene list that was reported in the same study, based on DNA methylation profiling obtained with the Illumina 27K BeadChip Methylation Array (5); however, this discrepancy can be easily accounted for by the larger genomic coverage of the 450K Methylation Array used in our study. Second, a widespread impact of aberrant DNA methylation in gene expression is expected by the identification of three miRNAs with marked DNA methylation changes in promoter regions, according to the growing list of their potential transcriptional targets.

In addition to a global impact on gene expression, DNA methylation changes had a selective impact in the promoter regions of important transcription factors downstream of TGF-β signaling: SMAD3 and the RUNX family of developmental factors. The latter family of factors is frequently deregulated in human cancers, where they may function as either oncogenes or tumor suppressors according to their context. Likewise, RUNX proteins are critical regulators of TGF-β signaling in a context-dependent fashion (45). However, the tumor promoting functions of RUNX proteins in TAFs are currently unknown, and overall their roles in lung cancer remain poorly understood. In contrast, low SMAD3 expression is frequently observed in the stroma of lung tumors of current smokers (46), in agreement with its epigenetic silencing observed in this study. In fibroblasts, SMAD3 is considered a critical transcription factor underlying TGF-β-induced α-SMA expression in lung and other organs (47,48). Yet, lung TAFs overexpressed α-SMA compared to CFs, indicating that α-SMA expression is largely SMAD3-independent in TAFs, and revealing that TAFs may have an aberrant activated phenotype. In support of this interpretation, lung TAFs exhibited enhanced contractility and expression of wound-associated ECM genes compared to CFs, in qualitative agreement with previous observations in skin fibroblasts in Smad3-null mice (36). These novel findings reveal that TAFs are hyperresponsive to exogenous TGF-β1, and strongly suggest that such hyperresponsiveness is intimately associated with the epigenetic silencing of SMAD3. Based on these observations, it is tempting to speculate that TGF-β1 hyperresponsiveness may be a major contributor to the formation and maintenance of the desmoplastic stroma in NSCLC (6).

Pathway enrichment analysis of DNA methylation differences pointed at pathways and processes associated with ECM/focal adhesions and FCγR. The former enrichment was expected, since ECM pathways are directly associated with the myofibroblast phenotype. In agreement with this analysis, enrichment in the KEGG focal adhesion pathway was previously found from differential expression data between lung TAFs and CFs elsewhere (5). Unexpectedly, our analysis identified also an enrichment in pathways related to the FCγR, which has been associated with immune cells rather than myofibroblasts, for its main known function is to provide antigen-presenting activity (34). Nonetheless there is evidence that bone marrow-derived fibrocytes do express FCγRs and can differentiate into fibroblasts/myofibroblasts (33,34), thereby raising the intriguing possibility that lung TAFs are enriched in differentiated fibrocytes. In support of this possibility, we found stromal cells in histologic sections with spindle-shaped nuclei and positivity for the hematopoietic proteins CD34 and CD45, which are well-defined fibrocyte markers. In agreement with these observations, the percentages of CD34+ and CD45+ fibroblasts in culture were more than ~2.5-fold larger in TAFs than CFs. Moreover, it is possible that the actual fibrocyte percentages were underestimated, since fibrocytes that differentiate into fibroblasts/myofibroblasts frequently down-regulate CD34 and/or CD45 (33,34). On the other hand, fibrocyte-like stromal cells have been documented in a growing list of human and murine tumors as well as in chronic inflammatory lung diseases (34,49). Therefore, these observations strongly support that a non-marginal percentage of lung TAFs derive from fibrocytes or fibrocyte-like cells in addition to resident fibroblasts.

More than a dozen hypermethylated genes have been previously identified as biomarker candidates in patients with NSCLC (10,13). Among them, only PCDHGB6 was also found consistently hypermethylated in TAFs, although mostly within non-promoter regions, indicating that most of the former hypermethylation events are restricted to either cancer cells or stromal cells other than fibroblasts. On the other hand, among the validated genes that exhibited a marked differential methylation in promoter regions in lung TAFs, only EDARADD was associated with a worse progression-free survival in NSCLC patients when analyzing public databases of DNA methylation and gene expression. EDARADD is a NF-kappa-B related developmental gene (50) that was found epigenetically upregulated in TAFs in culture and in histologic sections, thereby emerging as a potential stromal biomarker. It also appeared upregulated in cancer cells compared to parenchymal cells, indicating that future studies are warranted to elucidate the tumorigenic roles of EDARADD expression in TAFs and cancer cells. Moreover, since no specific survival data of NSCLC patients concomitant with DNA methylation in TAFs are available, it is conceivable that additional differentially methylated genes other than EDARADD may bear prognostic value in the stroma.

Conclusions

There is an increasing interest in understanding the molecular alterations of TAFs, since they are the major cell type within the tumor-supporting stroma. Our study uncovered unprecedented epigenetic alterations in TAFs from NSCLC patients, including a global impact on gene expression and a selective impact on SMAD3 and other key transcription factors of the TGF-β pathway. Remarkably, these alterations were associated with a hyperresponse to exogenous TGF-β, which may account for the sustained desmoplastic stroma in vivo. Our analysis revealed also that a fraction of TAFs may be recruited from bone-marrow fibrocytes. In addition, we identified EDARADD as a stromal biomarker in terms of prognosis. Altogether, our findings provide novel molecular and cellular insights underlying the tumor-promoting phenotype of lung TAFs.

Supplementary Material

Supplementary Figures 1–5 and Supplementary Tables 1–7 can be found at http://carin.oxfordjournals.org/

Funding

European Community’s Seventh Framework Programme (HEALTH-F2-2010–258677–CURELUNG to M.E.); the Ministerio de Economía y Competitividad co-financed by the European Development Regional Fund (PI13/02368 to J.A., PI10/02992 and PIE13/00022 to M.E., PS09/01377 to N.R., TEC2014-60337-R to A.P.); AECC (10/103 to N.R.-J.A.); SEPAR (2009–853 to N.R.); European Research Council (CoG-616480 to X.T.); Fundación Fuentes Abiertas (“New approaches in the fight against lung cancer” to J.A.); fellowships from the Ministerio de Educación (FPU to M.V.); Fundació Cellex (to M.P.); CONACYT (to R.L.); COLCIENCIAS (to A.V.); Marie-Curie action (CAFFORCE 328664 to A.L.). CIBERES and CIBER-BBN are a initiative of the Spanish ISCIII.