Cellular Model of Colon Cancer Progression Reveals Signatures of mRNAs, miRNA, lncRNAs, and Epigenetic Modifications Associated with Metastasis Free

Metastasis accounts for over 90% of cancer mortality. Therefore, a better understanding of its molecular and cellular bases is of high clinical relevance. Several studies suggested that reversible epithelial-to-mesenchymal transitions (EMT) promote metastasis: Cancer cells first convert from an epithelial- to a mesenchymal-like cell state (EMT) that allows them to leave the primary tumor and migrate to distant sites. When they colonize distant organs adoption of an epithelial-like state via a mesenchymal-to-epithelial transition (MET) allows the growth of metastases (1–4). Similar EMTs also occur during organ morphogenesis. For example, during renal organogenesis the mesenchyme surrounding the ureteric bud, which develops through EMT, gives rise to kidney epithelium via an MET process (5). The required cancer cell plasticity that allows the conversion between the metastable epithelial and mesenchymal states is presumably facilitated by epigenetic mechanisms, which stabilize expression patterns in each state, but still allow conversions into another cellular state (6, 7). The underlying epigenetic switches are triggered by intrinsic or extrinsic factors, but the resulting cellular state is inherited in the absence of the initial signal (8).

In this study, we developed a cellular model of colon cancer progression that recapitulates the alternating EMT/MET processes observed during metastasis. By characterization of the underlying transcriptomic and epigenetic alterations, combined with evaluation of their clinical relevance using patient-derived data, we uncovered numerous regulations and epigenetic modifications that may critically contribute to progression of colon cancer.

DLD1 colon cancer cells were from own stocks and authenticated by STR analysis in 2014 (Eurofins Medigenomix Forensik GmbH, Ebersberg, Germany). DLD1 cells and its derivatives were maintained in McCoy's 5A Medium (Invitrogen) containing 10% FBS (Invitrogen), 100 U/mL penicillin and 0.1 mg/mL streptomycin. SiRNAs [Ambion silencer siRNA: negative control (ID#4611), H19 (ID#n272450)], and miRNA inhibitors [Ambion Anti-miR: negative control (ID#AM17010), anti-miR-100-5p (ID#AM10188), anti-miR-125b-5p (ID#AM10148), and anti-miR-99a (ID#AM10719)] were transfected at a final concentration of 10 nmol/L using Lipofectamine 2000 transfection reagent (Invitrogen). IL6 (Immunotools) was dissolved in water and used at a final concentration of 20 ng/mL. 5-Aza-2′-deoxycytidine (Sigma) was used at a final concentration of 5 μmol/L.

To assess anchorage independent growth, 1,000 cells were seeded in medium containing 0.4% agarose on top of 0.8% agarose/media (GenePure LowMelt Agarose, ISC BioExpress). Covering media were changed twice per week. Colonies were stained with 0.01% Crystal Violet (Alfa Aesar) and examined after 21 days of culture. Colonies larger than 50 μm were enumerated.

Cells (5,000) were seeded in 24-well ultra-low attachment plates in sphere culture medium [DMEM-F12, B27 Supplement (1:50), EGF (20 ng/mL), BSA (0.4%), and Insulin (4 μg/mL)]. After 7 days, the resulting spheres were counted by phase-contrast microscopy, dissociated into single cells using a 0.05% Tryple solution (Invitrogen), and reseeded in sphere medium. This procedure was repeated for five times (generations).

The association of expression with nodal status or metastasis was calculated by a Student t test. qPCR, soft agar colony, colonosphere, migration, invasion, and TOP/FOP Wnt activity data are expressed as mean ± SD. Differences were analyzed by a two tailed Student t test. Calculations were performed using Prism5 (GraphPad Software Inc.) and P values ≤ 0.05 were considered as significant; *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001.

RNA-seq, miRNA, DNA methylation, and H3K4me3, H3K27me3, and H3K79me3 ChIP-seq data generated in this study have been deposited at the NCBI Gene Expression Omnibus repository (GEO) and is accessible through the accession number GSE85688.

Recently, we showed that IL6 treatment induces EMT in the epithelial colon cancer cell line DLD1 (9). To test whether the IL6-induced EMT is transient or permanent, we performed long-term treatment of DLD1 cells with IL6 followed by withdrawal of IL6. The upregulation of the mesenchymal cell state–associated mRNAs VIM, SNAIL, SLUG, and ZEB1 as well as downregulation of the epithelial cell state–associated mRNA CDH1 (encoding for E-cadherin) observed after IL6 exposure for 3 weeks was reversed after IL6 withdrawal (Fig. 1A). Similar effects were observed on the protein level: E-cadherin was reexpressed and F-Actin relocated to cell membranes after IL6 withdrawal (Supplementary Fig. S1A). Therefore, DLD1 cells underwent a transient EMT after IL6 treatment. It has been reported that not all cells switch from an epithelial-like state to a permanent mesenchymal-like state at once in primary tumors, but rather single permanently mesenchymal-like cell clones may give rise to a tumor invasion front (2, 10). Therefore, we hypothesized that during IL6-induced EMT single DLD1 cells might undergo a permanent EMT and retain mesenchymal characteristics after IL6 withdrawal. Because EMT has been associated with anchorage-independent growth and stemness (11, 12), we used the cultivation in soft agar to enrich for anchorage-independent growth or colonosphere culture conditions to enrich tumor stem cells. Treatment with IL6 resulted in a significant increase of soft agar colonies (Fig. 1B) and colonospheres (Fig. 1C). Importantly, cells isolated from single soft agar colonies (designated SA cells) or colonosphere pools (SPH cells) exhibited a mesenchymal morphology (spindle-shaped with scattered growth pattern) in monolayer culture in the absence of IL6 for more than 2 months, suggesting that these cells had undergone a permanent EMT (Fig. 1D). In contrast, vehicle-treated cells isolated from soft agar colonies or colonospheres showed an epithelial-like morphology (cobble-stone shaped and tightly packed) similar to parental DLD1 cells (data not shown). One explanation for the permanently mesenchymal phenotype of SA and SPH cells would be autocrine stimulation by IL6. However, in these cells, expression of IL6 mRNA or protein was not detectable (data not shown). To test the metastatic potential of parental DLD1 cells and the permanently mesenchymal DLD1 cell derivatives the respective cells were injected into the tail veins of immune-deficient NOD/SCID mice. In contrast with the poorly metastatic parental DLD1 cells, which gave rise to lung metastases in only one of 8 mice, SA and SPH cells were highly metastatic, suggesting that the mesenchymal phenotype enhances the metastatic potential (Fig. 1E). We isolated cells from metastases formed by mesenchymal-like SA cells and cultured them under monolayer conditions. Notably, these cells (designated METS) showed an epithelial morphology similar to parental DLD1 cells, suggesting that they had undergone an MET during metastatic colonization (Fig. 1D). Interestingly, although the METS cells displayed an epithelial morphology, they were highly metastatic (Fig. 1E). SA and SPH cells displayed a significantly higher expression of the mesenchymal markers VIM, SNAIL, SLUG, and ZEB1 as well as lower expression of the epithelial marker E-cadherin than DLD1 and METS cells (Fig. 1F and G). Furthermore, SA and SPH derivatives showed a higher invasion and migration when compared with DLD1 and METS cells (Fig. 1H; Supplementary Fig. S1B). In Matrigel, DLD1 and METS cells formed well-defined acinar structures, whereas SA and SPH cells formed disorganized structures with outgrowths reminiscent of invasive filopodia (Supplementary Fig. S1C; ref. 13). Altogether, the parental DLD1 cells and their SA, SPH, and METS derivatives represent a model of colon cancer progression that recapitulates the EMT/MET switching, which has been observed during metastatic progression: DLD1 cells represent the epithelial-like primary tumor cells, SA and SPH cells the invasive, migratory, mesenchymal-like tumor cells that had undergone EMT, and METS cells the epithelial-like tumor cells that form metastases after seeding and MET (Fig. 1D). Notably, SA and SPH cells may also represent circulating tumor cells (CTC), which often display mesenchymal properties (14).

Next, we performed RNA-Seq analyses to determine the mRNA and lncRNA expression profiles of DLD1, SA, SPH, and METS cells. Hierarchical clustering and principal component analyses revealed that SA and SPH cells showed high similarities in their mRNA expression profiles, whereas the profiles of DLD1 and METS cells were substantially different from all other cell derivatives (Fig. 2A and B). These results indicated that METS cells did not fully revert to the state of parental DLD1 cells during metastasis. Altogether, six mRNA signatures were identified: SSM—mRNAs more than 2-fold upregulated in SA, SPH, and METS cells compared with DLD1 cells (N = 33), MES—mRNAs more than 2-fold upregulated in mesenchymal-like SA and SPH cells compared with DLD1 and METS cells (N = 146), DLD—mRNAs more than 2-fold upregulated in DLD1 cells compared with SA, SPH, and METS cells (N = 30), DSS—mRNAs more than 2-fold upregulated in DLD1, SA, SPH cells compared with METS cells (N = 186), METS—mRNAs more than 2-fold upregulated in METS cells compared to DLD1, SA, and SPH cells (N = 80), and DM—mRNAs more than 2-fold upregulated in DLD1 and METS cells compared with SA and SPH cells (N = 34; Fig. 2A; Supplementary Fig. S2A; and Supplementary Dataset S1). Next, we tested whether the mRNA expression profiles of DLD1 cell derivatives correspond to epithelial- or mesenchymal-like cell states of established colon cancer cell lines. To obtain expression profiles of established epithelial- and mesenchymal-like colon cancer cell lines, we used data from the Cancer Cell Line Encyclopedia (CCLE) database and grouped the CCLE colon cancer cell lines (N = 58) according to their epithelial/mesenchymal characteristics: epithelial-like cell lines (high expression of CDH1 and low expression of VIM; N = 35) and mesenchymal-like cell lines (low expression of CDH1 and high expression of VIM; N = 8; Supplementary Fig. S2B). Gene set enrichment analysis (GSEA) revealed highly significant enrichment of mesenchymal-like SA and SPH cell-associated mRNAs (MES mRNA signature) in mesenchymal-like CCLE cell lines (Fig. 2C), whereas other signatures did not show a notable enrichment in epithelial- or mesenchymal-like CCLE cell lines (Fig. 2H and M; Supplementary Fig. S2C, S2F, and S2I). Next, we associated the identified mRNA signatures with clinical and pathological parameters by utilizing data from The Cancer Genome Atlas (TCGA) cohort of patients with colon adenocarcinoma (TCGA COAD, N = 456). The MES, SSM, and METS signatures, consisting of mRNAs with elevated expression in highly-metastatic SA, SPH, and METS cells, showed significant associations with poor overall and relapse-free survival (Fig. 2D, E, I, J, N, and O). Moreover, the expression of these signatures was elevated in primary tumors from patients with distant metastases and positive nodal status (Fig. 2F, G, K, L, P and R). In contrast, DSS, DM, and DLD signatures, consisting of mRNAs with elevated expression in low-metastatic DLD1 cells, did not show significant associations with clinical or pathological parameters (Supplementary Fig. S2D, S2E, S2G, S2H, S2J, and S2K). In two independent cohorts of 566 (GSE39582) and 130 (GSE37892) cases, SSM, MES, and METS signatures were significantly associated with poor relapse-free survival (Supplementary Fig. S2L–S2N). Altogether, the mRNA expression patterns of metastatic SA, SPH, or METS cells mirror the expression patterns of primary tumors from colon cancer patients with distant metastasis.

The miRNA expression profiles of the DLD1 derivatives could be separated into DLD1 and METS as well as SA and SPH groups (Fig. 3A and B). Unexpectedly, only eight miRNAs displayed significant differential expression that could be validated by qPCR analysis (Fig. 3C). The expression of miR-1246, miR-224, and miR-452 was elevated in DLD1 and METS cells (DM miRNAs), miR-675 and miR-135b were preferentially increased in METS cells (METS miRNAs), and miR-125b, miR-100, and miR-99a expression was induced in SA and SPH cells (MES miRNAs). miR-224 and miR-452 are encoded by the same host gene (15). miR-100 and miR-125b-1 are encoded by MIR100HG (hsa11), and miR-99a, and miR-125b-2 are encoded by MIR99AHG (hsa21), with miR-125b-1 and miR-125b-2 having the same mature miRNA sequence (16). In accordance with their common host genes, the expression levels of miR-224 and miR-452 as well as miR-100, miR-99a, and miR-125b show a high degree of correlation among each other within the TCGA COAD cohort (Fig. 3D). The negative correlation between DM and MES miRNA expression observed in the cellular model was also detected in the patient tumors from the TCGA COAD cohort (Fig. 3D). The MES miRNA signature was significantly associated with poor overall and relapse-free survival (Fig. 3E and F). Moreover, all three MES miRNAs were also individually significantly associated with poor survival (Supplementary Fig. S3A–S3C). The DM miRNA signature was associated with good survival, although with borderline significance (Fig. 3G and H), whereas no association of the METS miRNA signature with survival was observed (Fig. 3I and J). Next, we inhibited these miRNAs with specific antagomirs. However, because the mature MES miRNAs (especially miR-100 and miR-99a) have very similar sequences, cross inhibition was observed (Supplementary Fig. S3D and S3E). Inhibition of the MES miRNAs miR-125b, miR-100, and miR-99a in SA and SPH cells suppressed cell migration and invasion to a similar degree (Fig. 3K and L; and Supplementary Fig. S3F). Moreover, the expression of the mesenchymal-state associated gene ZEB1 was suppressed after inhibition of miR-125b, miR-100, and miR-99a (Fig. 3M). Next, we aimed to determine potential target mRNAs of the differentially regulated miRNAs that might represent mediators of their effects on tumor progression. Because miRNAs repress their target mRNAs, potential targets of MES miRNAs are expected among DM mRNAs, whereas potential targets of DM miRNAs might be found among MES mRNAs. Indeed, the negative correlation between the expression levels of the majority of MES mRNAs with the expression of DM miRNAs in the cell model was reflected in patient tumors from the TCGA COAD cohort (Fig. 3N). Moreover, the expression of most of these MES mRNAs was also positively correlated with the expression of MES miRNAs in the TCGA cohort (Fig. 3N). On the contrary, the association between MES miRNA and DM mRNA expressions was not nearly as convincing, since only few of MES miRNA/DM mRNA pairs showed a negative correlation in the TCGA cohort (Fig. 3N). By combining the miRNA/mRNA correlation data from the cellular model and TCGA with miRNA target prediction algorithms, we identified 21, 43, and 35 potential targets of DM miRNAs miR-1246, miR-224, and miR-452 and 3, 4, and 3 potential targets of MES miRNAs miR-100, miR-125b, and miR-99a, respectively (Supplementary Tables S1 and S2). For example, the mRNA encoding the metastasis-associated protein SERPINE1 (17, 18) harbors a conserved miR-224 binding site in the 3′-untranslated region (Supplementary Fig. S3G). Furthermore, expression of SERPINE1 mRNA inversely correlated with expression of miR-224 in the TCGA COAD dataset (Supplementary Fig. S3H). Finally, in SA and SPH cells that express low levels of miR-224 the levels of SERPINE1 mRNA were increased (Supplementary Fig. S3I).

The phenotype of each of the DLD1 derivatives was stable, as they did not change their characteristics during continuous passaging (data not shown). The maintenance of such stable cellular phenotypes might be explained by epigenetic mechanisms that allow metastable states. Analyses of DNA methylation patterns in the cellular model showed that most CpG sites exhibit either a low or high degree of methylation (Supplementary Fig. S4A). Furthermore, global methylation was the highest in SA and SPH cells and lowest in METS cells (Fig. 4A). Principal component analyses and hierarchical clustering once more demonstrated a separation into DLD1/METS and SA/SPH groups (Fig. 4B and C). In addition, we observed a strong negative correlation between expression and methylation in gene regulatory regions (200,178 CpG sites) in all four DLD1 derivatives (Fig. 4D and Supplementary Fig. S4B–S4D). Methylation in gene-regulatory regions was significantly lower in SA, SPH, and METS cells when compared with DLD1 cells. Moreover, methylation was further decreased in METS cells versus SA and SPH cells (Fig. 4E). Also in samples from colon cancer patients within the TCGA database, global methylation in regulatory regions was significantly lower in liver metastasis when compared with the corresponding primary tumor (Fig. 4F). Moreover, methylation in regulatory regions was lower in primary tumors from patients with liver metastasis compared with patients without metastasis (Fig. 4G). Subsequently, we determined whether the changes in gene expression (Δexpression) between cell states are accompanied by alterations in DNA methylation in gene-regulatory regions (Δmethylation) of those genes. Because SA and SPH cells had very similar expression and methylation patterns we calculated the average expression and methylation in those cells and used it as one parameter (MES). For all possible comparisons between cells (MES vs. DLD1, MES vs. METS, and METS vs. DLD1), we observed a negative correlation between Δexpression and Δmethylation with the strongest correlations detected when only genes with Δexpression of more than 2-fold and Δmethylation more than 10% were considered (Fig. 4H). Lists of genes that exhibit highest differential expression and DNA methylation between the DLD1 derivatives are provided in Supplementary Tables S3–S5. A list of all genes with differential expression and changes in DNA methylation between DLD1 derivatives is provided in Supplementary Dataset S2.

Besides DNA methylation, also histone modifications contribute to epigenetic regulation of gene expression. We analyzed three histone modifications: H3K4me3, which is associated with initiation of transcription, H3K79me3, which is linked to transcriptional elongation, and H3K27me3, which is characteristic for transcriptional repression. Indeed, in DLD1 cells, enrichment of H3K4me3 and H3K79me3 at the transcription start site (TSS) was associated with elevated expression, whereas enrichment of H3K27me3 correlated with low expression (Supplementary Fig. S5A). In SA, SPH, and METS cells the associations between histone modifications and expression were similar (data not shown). Differential expression among the different cell derivatives positively correlated with changes in H3K4me3 and H3K79me3 marks at the TSSs of the respective genes (Fig. 5A and B). Lists of all genes with differential expression and changes in H3K4me3/H3K79me3 marks between DLD1 derivatives are provided in Supplementary Datasets S3 and S4. In contrast, changes in expression negatively correlated with alterations in H3K27me3 (Supplementary Fig. S5B). A list of all genes with differential expression and H3K27me3 between DLD1 derivatives is provided in Supplementary Dataset S5). Similar to DNA methylation, the strongest correlations between changes in expression and changes in histone modifications were observed for Δexpression and Δhistone modification of more than 2-fold. Altogether, these results suggest that changes in gene expression that occur during in the cellular model of tumor progression are frequently accompanied by alterations in DNA methylation or histone modifications.

Recently, several long noncoding RNAs (lncRNA) have been implicated in metastasis (19, 20). By RNA-Seq, we detected the expression of 1,427 lncRNAs in DLD1, SA, SPH, and METS cells. Hierarchical clustering and principal component analyses of lncRNA expression separated the DLD1-derivatives into DLD1/METS and SA/SPH sub-groups (Fig. 6A and B). Three main lncRNA signatures were identified: MES, 55 lncRNAs upregulated more than 2-fold in mesenchymal-like SA and SPH cells compared with DLD1 and METS cells; METS, 23 lncRNAs upregulated more than 2-fold in METS cells compared with DLD1, SA, and SPH cells; and DLD, 18 lncRNAs upregulated more than 2-fold in DLD1 cells compared with SA, SPH, and METS cells (Fig. 6A; Supplementary Dataset S6). lncRNAs with elevated expression in SA and SPH cells (MES lncRNAs) were significantly associated with poor survival, metastasis, and positive nodal status in the TCGA COAD cohort (Fig. 6C–F). Among these were several known oncogenic lncRNAs, such as HOTAIR and PVT-1 (Supplementary Dataset S6), which are known to promote metastasis in various cancer types including colon cancer (19, 21, 22). The METS lncRNA signature showed significant association with metastasis and positive nodal status, but not with survival (Supplementary Fig. S6A–S6D), whereas lincRNAs with highest expression in DLD1 cells did not show significant associations with clinical parameters (Supplementary Fig. S6E–S6H).

The most pronounced difference in lncRNA expression between DLD1 derivatives was observed for H19. Although H19 expression was lower in SA and SPH cells when compared with DLD1 cells, H19 was highly expressed in METS cells (Fig. 6G). The elevated expression of H19 in METS cells was accompanied with notably lower DNA methylation of the H19 gene region in METS cells compared with other cells (Fig. 6G). Furthermore, comparison of primary tumors and the corresponding liver metastases from within the TCGA database also demonstrated higher expression and lower methylation of H19 in colon cancer metastases (Fig. 6H). Treatment with the methylation inhibitor 5-Aza-2′-deoxycytidine resulted in the induction of H19 expression DLD1, SA, SPH cells, but not in METS cells (Fig. 6I) implying that DNA methylation in the H19 gene region regulates H19 expression. The enrichment of histone modifications H3K4me3 and H3K79me3 at the H19 TSS was strongly increased in METS cells (Fig. 6J and K). Knockdown of H19 in DLD1 and METS cells resulted in the repression of the epithelial marker CDH1 and induction of mesenchymal markers SLUG and ZEB1 (Fig. 6L and M), suggesting that H19 suppresses EMT. H19 expression was also significantly elevated in SW620 compared to SW480 cells, which represent lymph node metastasis and primary tumor cells from the same patient, respectively (Supplementary Fig. S6I). In 18 pairs of primary colon cancer tumors and matched liver metastases (23), H19 expression was significantly higher in metastases (Supplementary Fig. S6J). Moreover, expression of H19 was significantly higher in the CT26 murine colon cancer cell line derivative that had undergone three cycles of inoculation and isolation from liver metastasis compared with parental CT26 cells (GSE67675; Supplementary Fig. S6K). Furthermore, elevated expression of H19 was associated with poor survival, metastasis, and positive nodal status in TCGA COAD, GSE39582, and GSE37892 cohorts (Supplementary Fig. S6L–S6R). Taken together, these results suggest that H19 expression is repressed via CpG methylation during EMT, but strongly induced during metastatic colonization. Therefore, H19 presumably functions as a suppressor of early and as a promoter of late stages of colon cancer progression.

Next, we used the molecular signatures database (MSigDB) to determine enriched signatures and pathways in the different DLD1 derivatives. Because, SA, SPH, and METS cells were metastatic in mice and their mRNA signatures associated with metastasis in patients, we selected mRNAs that showed at least 1.5-fold elevated expression when compared with parental DLD1 cells. Gene set enrichment analyses of the MSigDB Hallmark gene sets showed the most pronounced enrichment within KRAS signaling, estrogen response, EMT, and hypoxia related gene sets (Supplementary Fig. S7A). To test whether the same gene sets are enriched in primary tumors from patients with metastasis (M1) compared with patients without metastasis (M0), we compared them with 2056 mRNAs with significantly higher expression in tumors from TCGA COAD M1 patients and 1,149 mRNAs that showed a significantly higher expression in M0 patients (Supplementary Fig. S7B). The MSigDB—Hallmark gene sets EMT, KRAS signaling, and hypoxia showed the highest enrichment in M1-associated mRNA (Supplementary Fig. S7C), thereby displaying a notable overlap with the enrichment of SA-, SPH-, and METS-associated mRNAs.

To determine which pathways are altered in the 4 DLD1 derivatives, we performed a global analysis of signaling pathways by using the Cignal 45-pathway reporter array, which allows the detection of 45 signal transduction pathways (Supplementary Fig. S7D). ATF2/3/4, ATF6, MAPK/JNK, Myc, NRF1/2, SP1, and Wnt pathways showed activation in the SA and SPH derivatives when compared with DLD1 and METS cells, whereas the estrogen receptor (ER) pathway showed repression in SA, SPH, and METS cells compared with DLD1 cells (Supplementary Fig. S7D). Higher Wnt activity in SA and SPH cells was verified by TOP/FOP assays (Supplementary Fig. S7E) and immunofluorescence analyses, which showed nuclear β-catenin staining in SA and SPH cells (Supplementary Fig. S7F). Consistent with these results we found that the gene set LEF1_UP.V1_UP showed the most significant enrichment in SA-, SPH-, METS-associated mRNAs (Supplementary Fig. S7G) and was also highly enriched in TCGA COAD M1-associated mRNAs (Supplementary Fig. S7H). This gene set contains the 200 most upregulated genes in DLD1 cells transfected with an LEF-1 expression construct compared with control DLD1 cells (GSE3229). LEF-1 is a Wnt-signaling effector, suggesting that the activation of Wnt-signaling is higher in mesenchymal-like cells. Finally, we analyzed the MSigDB – Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway gene sets and surprisingly, the most significant enrichment of SA, SPH, METS-associated mRNAs was in the axon guidance pathway (Fig. 7A). Interestingly, also the M1-associated mRNAs displayed this association (Fig. 7B). Notably, the majority of the 31 axon guidance genes that are expressed at higher levels in M1 TCGA tumors show lower methylation in these tumors in the TCGA COAD dataset (Supplementary Fig. S8A), suggesting that a decrease in DNA methylation and consequent increase of expression of axon guidance genes is an important event during colon cancer metastatic progression. The KEGG axon guidance expression signature was significantly associated with poor overall and relapse-free survival, as well as with metastasis and positive nodal status in the TCGA COAD cohort (Fig. 7C–F), indicating a potentially important role of this pathway in colon cancer metastasis. Moreover, when SA and SPH cells were grown in Matrigel, protrusions that resemble the filopodia of axon growth cones were observed (Supplementary Fig. S1D). Among the SA-, SPH-, METS- and TCGA COAD M1–associated axon guidance genes we identified several semaphorins and their receptors, the plexins. For example, plexin B3 (PLXNB3) and semaphorin 4A (SEMA4A) were expressed at notably higher levels in highly metastatic SA, SPH, and METS cells when compared with largely nonmetastatic DLD1 cells (Fig. 7G; Supplementary Fig. S9A). High expression of PLXNB3 and SEMA4A was associated with poor survival, metastasis, and positive nodal status in TCGA COAD, GSE39582, and GSE37892 cohorts (Fig. 7H–K; Supplementary Figs. S8B and S8C and S9B–S9F). Furthermore, the enrichment of the histone marks H3K4me3 and H3K79me3, which are associated with transcriptionally active genes, was higher at the TSS of PLXNB3 and SEMA4A in SA, SPH, and METS cells (Fig. 7L and M; and Supplementary Fig. S9G). Moreover, DNA methylation in the regulatory regions of PLXNB3 and SEMA4A genes was lower in SA, SPH, and METS cells (Fig. 7N; Supplementary Fig. S9H). In addition, the methylation of PLXNB3 and SEMA4A gene regulatory regions was lower in tumors from patients with distant metastasis when compared with patients without metastasis in the TCGA COAD cohort (Fig. 7O; Supplementary Fig. S9I). Therefore, epigenetic mechanisms are presumably involved in the upregulation of PLXNB3 and SEMA4A during metastatic progression. Collectively, these results suggest that the activation of axon guidance-associated genes during colon cancer progression may play an important role in metastasis.

In the cellular model of colon cancer progression introduced here, we observed reversible EMT/MET transitions during the metastatic cascade. The mesenchymal-like tumor cells were highly metastatic in mice and the corresponding transcriptomic profile was associated with metastasis in patients, indicating that the transcriptional programs were translated into cellular functions that contribute to metastatic traits. Unexpectedly, the metastatic potential of epithelial-like METS cells remained high after the MET that occurred after colonization of mouse lungs. Therefore, the prometastatic potential of METS cells was retained despite the reversion to an epithelial-like morphology. Interestingly, the transcriptional and epigenetic profiles of METS diverged significantly from that of the parental DLD1 cells and may explain their enhanced metastatic potential.

Surprisingly, the most significant upregulations associated with metastasis were observed for axon guidance–associated genes. Interestingly, semaphorins and their receptors, the plexins, which play a major role in the directed migration of axons, have previously been implicated in the migration of tumor cells (24). Navigating axons are equipped with a terminal growth cone containing F-actin–rich filopodia and lamellipodia, which extends and retracts to explore the environment and allows directional migration (25). Adaptation of such mechanisms by tumor cells might be advantageous during invasion, intravasation, and extravasation.

The transitions between different states of the cellular model resulted in relatively stable cell derivatives, implying that genetic or epigenetic mechanisms are involved. Because the transitions occurred in a short period of time, genetic alterations are unlikely to cause them. Moreover, the cells partially reversed to their initial state during metastatic colonization, suggesting that epigenetic mechanisms are likely to be involved in the maintenance of the states as they would confer the necessary degree of plasticity. Recent studies showed that epigenetic modifications can be maintained by transcription factors or noncoding RNAs. For example, the HOTAIR lncRNA, which shows elevated expression in SA and SPH cells, recruits PcG (polycomb group) and LSD1 complexes, which facilitates the methylation of H3K27 and demethylation of H3K4, respectively (26). In breast cancer, HOTAIR is strongly overexpressed in metastatic tumors, which facilitates genome-wide retargeting of PcG factors and establishment of new H3K27me3 histone marks in hundreds of genes, including various metastasis-suppressors (19). Thus, the transcriptomic and epigenetic alterations identified here are likely to include factors, which may form epigenetic switches (27) that regulate the EMT/MET transitions during metastatic progression.

Global loss of DNA methylation is one of the characteristic features of cancer cells that is found early in carcinogenesis, but is also often associated with tumor progression (28–31). For example, global methylation of long interspersed nuclear element-1 regions is significantly lower in liver metastases than in primary colon cancers (32). Our results show that METS cells display the lowest global methylation as well as the lowest methylation in gene regulatory regions. Furthermore, the number of genes with a loss in methylation was higher during the transition from SA/SPH to METS than during the transition from DLD1 to SA/SPH cells. Altogether, these results suggest that most changes in DNA methylation during cancer progression occur during metastatic colonization.

Deregulated expression of H19 lncRNA has been implicated in various malignancies (33). However, H19 seems to exert both oncogenic and tumor-suppressive functions (34). In our cellular model, H19 expression decreased from the parental DLD1 cells to the mesenchymal SA and SPH derivatives. However, in METS cells, H19 expression was highly induced and accompanied by lower DNA methylation and increased H3K4me3 and H3K79me3 histone marks in H19 gene region. Because the conversion from DLD1 to SA/SPH represents early stages of metastasis and the conversion from SA/SPH to METS late stages, our findings indicate that H19 may function as a tumor suppressor in the early stages of metastasis, but has oncogenic properties in the late stages, such as promoting colonization.

In conclusion, our genome-wide transcriptomic and epigenetic characterization of the cellular model of colon cancer metastasis introduced here represents a comprehensive resource for enhancing our understanding of colon cancer metastasis. In the future, further functional characterization and validation of the leads identified here will allow to identify clinically relevant biomarkers and intervention points for the treatment of colon cancer.

No potential conflicts of interest were disclosed.

Conception and design: M. Rokavec, H. Hermeking

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M. Rokavec, D. Horst

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M. Rokavec, H. Hermeking

Writing, review, and/or revision of the manuscript: M. Rokavec, H. Hermeking

Study supervision: H. Hermeking

Other (performed experiments and analyses): M. Rokavec

We thank the microarray unit of the DKFZ Genomics and Proteomics Core Facility for providing the miRNA arrays and Illumina Human Methylation arrays and related services, Lin Zhang and Xiaohui Yan for TCGA COAD–derived lncRNA expression profiles. We are also grateful to Rene Jackstadt and Stefanie Jaitner for help with the analysis of NOD/SCID mice, and to Andreas Jung and Sabine Sagebiel-Kohler for assistance with the analysis of ChIP-seq library quality.

This work was supported by grants from the Else Kröner-Fresenius-Stiftung, the Rudolf-Bartling-Stiftung, and the DKTK.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.