Here, we developed and comprehensively characterized a cellular model of colon cancer progression consisting of four defined derivatives of a colon cancer cell line that resulted from consecutive epithelial–mesenchymal and mesenchymal–epithelial transitions (EMT/MET) and phenotypically recapitulate the metastatic cascade. Initial EMT was induced by prolonged exposure to IL6, a cytokine also generated by the tumor-stroma. Genome-wide characterization of transcriptional (mRNA, miRNA, and lncRNA) and epigenetic (DNA methylation, H3K4me3, H3K79me3, and H3K27me3 histone modifications) profiles of the cell derivatives, combined with correlative analyses of expression, methylation, and clinical data from the TCGA-COAD database gave insights into the molecular basis of their phenotypic changes. The signatures characterizing invasive, mesenchymal-like cell states as well as the metastases-derived epithelial-like state showed significant association with metastasis, positive nodal status, and poor survival of colon cancer patients. Global hypomethylation of gene-regulatory regions was observed during tumor progression, with the lowest degree of methylation present in cells isolated from metastases. Upregulation of an axon-guidance–related gene signature was the most significant feature of metastatic tumor cells and was also found in primary tumors from colon cancer patients with distant metastases. Furthermore, the microRNAs miR-99a, miR-100, and miR-125b showed elevated expression in mesenchymal-like cells, associated with poor survival, and promoted migration and invasion. Finally, elevated expression of H19 lncRNA due to promoter demethylation was observed in cells isolated from metastases and was associated with poor survival of colon cancer patients. In the future, our results may be further exploited for the discovery and evaluation of novel metastasis-associated mechanisms and biomarkers. Cancer Res; 77(8); 1854–67. ©2017 AACR.

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.

Cell lines, cell culture, and reagents

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.

Soft agar colony assay

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.

Colonosphere assay

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).

Statistical analysis

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.

Data and materials availability

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.

Establishment of a cellular model of colon cancer metastatic progression

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).

Figure 1.

Establishment of a cellular model that recapitulates the EMT/MET switching during colon cancer metastasis. A, qPCR analysis of the indicated mRNAs in DLD1 cells treated with IL6 for 3 weeks, followed by 3 weeks of IL6 withdrawal. B, Soft agar colony formation assay. DLD1 cells were cultured in soft agar and treated with vehicle or 20 ng/mL IL6 for 3 weeks; scale bar, 2 mm; left, representative picture showing soft agar colonies; right, quantification of colonies larger than 50 μm. C, Colonosphere formation assay. DLD1 cells were cultured in colonosphere medium containing vehicle or 20 ng/mL IL6. Each generation represents 7 days of culture; left, representative picture showing spheres; right, quantification of colonospheres; scale bar, 50 μm. D, Phase contrast microscopy of the indicated derivates of DLD1 cells (top part); scale bar, 100 μm. Summary of the cellular model of colon cancer progression consisting of DLD1 cell derivatives that recapitulate EMT/MET transitions during metastasis (bottom part). E, Lung metastasis formation after tail vein injection of indicated DLD1 derivatives in NOD/SCID mice (top part). Percentage and number of mice with metastases/total number of analyzed mice is indicated. Hematoxylin and eosin staining of lungs from mice nine weeks after injection (pictures); scale bars, 200 μm. F, qPCR analyses of the indicated mRNAs in DLD1 cell derivatives. G, Immunofluorescence detections of E-cadherin in DLD1 cell derivatives. Nuclear DNA was detected by staining with DAPI; scale bar, 50 μm. H, Invasion assay in a modified Boyden-chamber. Cells were allowed to migrate through Matrigel-coated filter for 48 hours. Mean values ± SD (n = 3) are provided with *, P < 0.05; **, P < 0.01; ***, P < 0.001.

Figure 1.

Establishment of a cellular model that recapitulates the EMT/MET switching during colon cancer metastasis. A, qPCR analysis of the indicated mRNAs in DLD1 cells treated with IL6 for 3 weeks, followed by 3 weeks of IL6 withdrawal. B, Soft agar colony formation assay. DLD1 cells were cultured in soft agar and treated with vehicle or 20 ng/mL IL6 for 3 weeks; scale bar, 2 mm; left, representative picture showing soft agar colonies; right, quantification of colonies larger than 50 μm. C, Colonosphere formation assay. DLD1 cells were cultured in colonosphere medium containing vehicle or 20 ng/mL IL6. Each generation represents 7 days of culture; left, representative picture showing spheres; right, quantification of colonospheres; scale bar, 50 μm. D, Phase contrast microscopy of the indicated derivates of DLD1 cells (top part); scale bar, 100 μm. Summary of the cellular model of colon cancer progression consisting of DLD1 cell derivatives that recapitulate EMT/MET transitions during metastasis (bottom part). E, Lung metastasis formation after tail vein injection of indicated DLD1 derivatives in NOD/SCID mice (top part). Percentage and number of mice with metastases/total number of analyzed mice is indicated. Hematoxylin and eosin staining of lungs from mice nine weeks after injection (pictures); scale bars, 200 μm. F, qPCR analyses of the indicated mRNAs in DLD1 cell derivatives. G, Immunofluorescence detections of E-cadherin in DLD1 cell derivatives. Nuclear DNA was detected by staining with DAPI; scale bar, 50 μm. H, Invasion assay in a modified Boyden-chamber. Cells were allowed to migrate through Matrigel-coated filter for 48 hours. Mean values ± SD (n = 3) are provided with *, P < 0.05; **, P < 0.01; ***, P < 0.001.

Close modal

Characterization of DLD1 derivatives by expression profiling

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.

Figure 2.

High expression of SA, SPH, and METS cell-associated mRNAs correlates with metastasis and poor survival in patients with colon adenocarcinoma. A, Heatmap and unsupervised hierarchical clustering analysis of mRNA expression in DLD1 cell derivatives. The top 10% of mRNAs with most variable expression values are shown. Six main gene signatures are indicated. B, Principal component analysis of mRNA expression in DLD1 cell derivatives. C, GSEA of the MES mRNA signature in the epithelial- versus mesenchymal-like colon cancer cell lines from the CCLE database. D and E, Association of the MES mRNA signature with overall (D) and relapse-free (E) survival in the TCGA COAD dataset. F and G, Association of the MES mRNA signature with metastasis (F) and nodal status (G) in the TCGA COAD dataset. Expression in primary tumors of patients with (M1) or without (M0) distant metastasis, or different levels of lymph node infiltration (N0–N1) is shown in H–L same as C–G, but for the SSM mRNA signature. M–R, Same as C–G, but for the METS mRNA signature. *, P < 0.05; ***, P < 0.001.

Figure 2.

High expression of SA, SPH, and METS cell-associated mRNAs correlates with metastasis and poor survival in patients with colon adenocarcinoma. A, Heatmap and unsupervised hierarchical clustering analysis of mRNA expression in DLD1 cell derivatives. The top 10% of mRNAs with most variable expression values are shown. Six main gene signatures are indicated. B, Principal component analysis of mRNA expression in DLD1 cell derivatives. C, GSEA of the MES mRNA signature in the epithelial- versus mesenchymal-like colon cancer cell lines from the CCLE database. D and E, Association of the MES mRNA signature with overall (D) and relapse-free (E) survival in the TCGA COAD dataset. F and G, Association of the MES mRNA signature with metastasis (F) and nodal status (G) in the TCGA COAD dataset. Expression in primary tumors of patients with (M1) or without (M0) distant metastasis, or different levels of lymph node infiltration (N0–N1) is shown in H–L same as C–G, but for the SSM mRNA signature. M–R, Same as C–G, but for the METS mRNA signature. *, P < 0.05; ***, P < 0.001.

Close modal

miRNA profiling of the colon cancer progression model

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).

Figure 3.

miRNAs of miR100HG and miR99aHG clusters are associated with mesenchymal-like tumor cell state and poor survival in patients with colon adenocarcinoma. A, Heatmap and unsupervised hierarchical clustering analysis of miRNA expression in DLD1 cell derivatives. miRNAs with more than 1.5-fold difference in expression are shown. B, Principal component analysis of differentially expressed miRNAs in DLD1 cell derivatives. C, qPCR analyses of indicated miRNAs in DLD1 cell derivatives. D, Correlations between expression of indicated miRNAs in the TCGA COAD dataset. E and F, Association of the MES miRNA signature with overall (E) and relapse-free (F) survival in the TCGA COAD dataset. G and H same as E and F, but for the DM miRNA signature. I and J same as E and F, but for the METS miRNA signature. K, Invasion assay in a modified Boyden-chamber. SA and SPH cells were transfected with indicated miRNA inhibitors for 48 hours and allowed to migrate through Matrigel-coated filter for another 48 hours. L, Wound-healing assay of SA and SPH cells transfected with indicated miRNA inhibitors for 72 hours. Quantification of the percentage of wound closure 24 hours after the scratch is shown. M, qPCR analyses of ZEB1 expression in SA and SPH cells transfected with indicated miRNA inhibitors for 72 hours. N, Correlation of expression of MES and DM mRNAs with expression of MES and DM miRNAs in the TCGA COAD dataset. In K-M, mean values ± SD (n = 3) are provided with *, P < 0.05; **, P < 0.01; ***, P < 0.001.

Figure 3.

miRNAs of miR100HG and miR99aHG clusters are associated with mesenchymal-like tumor cell state and poor survival in patients with colon adenocarcinoma. A, Heatmap and unsupervised hierarchical clustering analysis of miRNA expression in DLD1 cell derivatives. miRNAs with more than 1.5-fold difference in expression are shown. B, Principal component analysis of differentially expressed miRNAs in DLD1 cell derivatives. C, qPCR analyses of indicated miRNAs in DLD1 cell derivatives. D, Correlations between expression of indicated miRNAs in the TCGA COAD dataset. E and F, Association of the MES miRNA signature with overall (E) and relapse-free (F) survival in the TCGA COAD dataset. G and H same as E and F, but for the DM miRNA signature. I and J same as E and F, but for the METS miRNA signature. K, Invasion assay in a modified Boyden-chamber. SA and SPH cells were transfected with indicated miRNA inhibitors for 48 hours and allowed to migrate through Matrigel-coated filter for another 48 hours. L, Wound-healing assay of SA and SPH cells transfected with indicated miRNA inhibitors for 72 hours. Quantification of the percentage of wound closure 24 hours after the scratch is shown. M, qPCR analyses of ZEB1 expression in SA and SPH cells transfected with indicated miRNA inhibitors for 72 hours. N, Correlation of expression of MES and DM mRNAs with expression of MES and DM miRNAs in the TCGA COAD dataset. In K-M, mean values ± SD (n = 3) are provided with *, P < 0.05; **, P < 0.01; ***, P < 0.001.

Close modal

Genome-wide epigenetic changes distinguish DLD1 derivatives

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.

Figure 4.

Metastatic colonization is associated with strong DNA hypomethylation. A, Box and whisker plot of global DNA methylation in DLD1 cell derivatives. From bottom to top, minimum, first quartile, median, third quartile, and maximum is shown. B, Principal component analysis of DNA methylation in DLD1 cell derivatives. C, Heatmap and unsupervised hierarchical clustering analysis of DNA methylation in DLD1 cell derivatives. The top 10% of methylation sites with most variable methylation are shown. D, Correlation plot of expression versus average DNA methylation in gene regulatory region of the corresponding gene in DLD1 cells. Expression/methylation of 14,881 genes was analyzed. E, Box and whisker plot of DNA methylation in gene regulatory regions in DLD1 cell derivatives. Indications are same as in A. F, box and whisker plot of DNA methylation in gene-regulatory regions in primary tumor and the corresponding liver metastasis of the NH-A8F7 patient from the TCGA COAD cohort. Indications are same as in A. G, Box and whisker plot of DNA methylation in gene-regulatory regions in primary tumors of TCGA COAD patients with (M1) or without (M0) distant metastasis. Indications are same as in A. H, Correlation plots of changes in expression (Δexpression) versus changes in DNA methylation in the regulatory regions (Δmethylation) of the corresponding gene between indicated DLD1 cell derivatives. Genes with Δexpression (log2 RPKM) < −1 or >1 and Δmethylation (β) < −0.1 or >0.1 are shown in red, genes with Δexpression (log2 RPKM) < −0.6 or >0.6 and Δmethylation (β) < −0.05 or >0.05 are shown in blue, and other genes are displayed in gray. ***, P < 0.001.

Figure 4.

Metastatic colonization is associated with strong DNA hypomethylation. A, Box and whisker plot of global DNA methylation in DLD1 cell derivatives. From bottom to top, minimum, first quartile, median, third quartile, and maximum is shown. B, Principal component analysis of DNA methylation in DLD1 cell derivatives. C, Heatmap and unsupervised hierarchical clustering analysis of DNA methylation in DLD1 cell derivatives. The top 10% of methylation sites with most variable methylation are shown. D, Correlation plot of expression versus average DNA methylation in gene regulatory region of the corresponding gene in DLD1 cells. Expression/methylation of 14,881 genes was analyzed. E, Box and whisker plot of DNA methylation in gene regulatory regions in DLD1 cell derivatives. Indications are same as in A. F, box and whisker plot of DNA methylation in gene-regulatory regions in primary tumor and the corresponding liver metastasis of the NH-A8F7 patient from the TCGA COAD cohort. Indications are same as in A. G, Box and whisker plot of DNA methylation in gene-regulatory regions in primary tumors of TCGA COAD patients with (M1) or without (M0) distant metastasis. Indications are same as in A. H, Correlation plots of changes in expression (Δexpression) versus changes in DNA methylation in the regulatory regions (Δmethylation) of the corresponding gene between indicated DLD1 cell derivatives. Genes with Δexpression (log2 RPKM) < −1 or >1 and Δmethylation (β) < −0.1 or >0.1 are shown in red, genes with Δexpression (log2 RPKM) < −0.6 or >0.6 and Δmethylation (β) < −0.05 or >0.05 are shown in blue, and other genes are displayed in gray. ***, P < 0.001.

Close modal

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.

Figure 5.

Alterations in H3K4me3 and H3K79me3 histone modifications are associated with changes in expression in DLD1 derivatives. A, Correlation plots of changes in expression (Δexpression) versus changes in H3K4me3 mark (ΔH3K4me3) in the corresponding gene between indicated DLD1 cell derivatives. B, Correlation plots of changes in expression (Δexpression) versus changes in H3K79me3 mark (ΔH3K79me3) in the corresponding gene between indicated DLD1 cell derivatives. Genes with Δexpression (log2 RPKM) < −1 or >1 and Δhistone modification (log2 reads) < −1 or >1 are shown in red, genes with Δexpression (log2 RPKM) < −0.6 or >0.6 and Δhistone modification (log2 reads) < −0.6 or >0.6 are shown in blue, and other genes are displayed in gray.

Figure 5.

Alterations in H3K4me3 and H3K79me3 histone modifications are associated with changes in expression in DLD1 derivatives. A, Correlation plots of changes in expression (Δexpression) versus changes in H3K4me3 mark (ΔH3K4me3) in the corresponding gene between indicated DLD1 cell derivatives. B, Correlation plots of changes in expression (Δexpression) versus changes in H3K79me3 mark (ΔH3K79me3) in the corresponding gene between indicated DLD1 cell derivatives. Genes with Δexpression (log2 RPKM) < −1 or >1 and Δhistone modification (log2 reads) < −1 or >1 are shown in red, genes with Δexpression (log2 RPKM) < −0.6 or >0.6 and Δhistone modification (log2 reads) < −0.6 or >0.6 are shown in blue, and other genes are displayed in gray.

Close modal

Association of H19 lncRNA expression with metastasis

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).

Figure 6.

The expression of H19 lncRNA is associated with colon cancer metastasis. A, Heatmap and unsupervised hierarchical clustering analysis of lncRNA expression in DLD1 cell derivatives. The top 10% of lncRNAs with most variable expression values is shown. Three main gene signatures are indicated. B, Principal component analysis of lncRNA expression in DLD1 cell derivatives. C and D, Association of the MES lncRNA signature with overall (C) and relapse-free (D) survival in the TCGA COAD dataset. E and F, Association of the MES lncRNA signature with metastasis (E) and nodal status (F) in the TCGA COAD dataset. G, qPCR analyses of H19 expression (left) and DNA methylation of H19 genomic region (right) in DLD1 cell derivatives. H,H19 expression (left) and DNA methylation of H19 genomic region (right) in primary tumor and the corresponding liver metastasis of the NH-A8F7 patient from the TCGA COAD cohort. I, qPCR analyses of H19 expression in indicated DLD1 cell derivatives after treatment with the DNA methylation inhibitor 5-Aza-2′-deoxycytidine for 1 week. J and K, Enrichment of H3K4me3 (J) and H3K79me3 (K) marks in the H19 genomic region of DLD1 cell derivatives. L and M, qPCR analyses of H19, CDH1, SLUG, and ZEB1 expression in DLD1 (L) and METS (M) cellstransfected with control of H19-specific siRNA for 72 hours. *, P < 0.05; ***, P < 0.01. In G, H, J, and K, hg19 genomic coordinates and the H19 gene structure with indicated exons (black boxes) and transcription start site (arrow) is shown at the bottom.

Figure 6.

The expression of H19 lncRNA is associated with colon cancer metastasis. A, Heatmap and unsupervised hierarchical clustering analysis of lncRNA expression in DLD1 cell derivatives. The top 10% of lncRNAs with most variable expression values is shown. Three main gene signatures are indicated. B, Principal component analysis of lncRNA expression in DLD1 cell derivatives. C and D, Association of the MES lncRNA signature with overall (C) and relapse-free (D) survival in the TCGA COAD dataset. E and F, Association of the MES lncRNA signature with metastasis (E) and nodal status (F) in the TCGA COAD dataset. G, qPCR analyses of H19 expression (left) and DNA methylation of H19 genomic region (right) in DLD1 cell derivatives. H,H19 expression (left) and DNA methylation of H19 genomic region (right) in primary tumor and the corresponding liver metastasis of the NH-A8F7 patient from the TCGA COAD cohort. I, qPCR analyses of H19 expression in indicated DLD1 cell derivatives after treatment with the DNA methylation inhibitor 5-Aza-2′-deoxycytidine for 1 week. J and K, Enrichment of H3K4me3 (J) and H3K79me3 (K) marks in the H19 genomic region of DLD1 cell derivatives. L and M, qPCR analyses of H19, CDH1, SLUG, and ZEB1 expression in DLD1 (L) and METS (M) cellstransfected with control of H19-specific siRNA for 72 hours. *, P < 0.05; ***, P < 0.01. In G, H, J, and K, hg19 genomic coordinates and the H19 gene structure with indicated exons (black boxes) and transcription start site (arrow) is shown at the bottom.

Close modal

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.

Axon guidance–associated genes as potential mediators of colon cancer metastasis

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.

Figure 7.

Axon guidance-associated genes are associated with colon cancer metastasis. A, GSEA analyses of mRNAs, which show more than 1.5-fold elevated expression in highly metastatic SA, SPH, and METS cells compared with low-metastatic DLD1 cells in KEGG pathway gene sets. The three most significant gene sets are shown. B, GSEA analyses of mRNAs, which display significantly elevated expression in TCGA COAD patients with metastasis compared with patients without metastasis in KEGG pathway gene sets. The three most significant gene sets are shown. C and D, Association of the KEGG axon guidance signature with overall (C) and relapse-free (D) survival in the TCGA COAD dataset. E and F, Association of the KEGG axon guidance signature with metastasis (E) and nodal status (F) in the TCGA COAD dataset. G, qPCR analyses of PLXNB3 expression in DLD1 cell derivatives. H and I, Association of PLXNB3 expression with overall (H) and relapse-free (I) survival in the TCGA COAD dataset. J and K, Association of PLXNB3 expression with metastasis (J) and nodal status (K) in the TCGA COAD dataset. L and M, Enrichment of H3K4me3 (L) and H3K79me3 (M) marks in the PLXNB3 genomic region in DLD1 cell derivatives. N, DNA methylation of the PLXNB3 genomic region in DLD1 cell derivatives. O, DNA methylation of the PLXNB3 genomic region in TCGA COAD patients with metastasis (M1) compared with patients without metastasis (M0). Mean values ± SD are provided with *, P < 0.05; **, P < 0.01; ***, P < 0.001. In L–O, hg19 genomic coordinates and the PLXNB3 gene structure with indicated exons (white boxes, untranslated region; black boxes, protein encoding region) and transcription start site (arrow) is shown at the bottom.

Figure 7.

Axon guidance-associated genes are associated with colon cancer metastasis. A, GSEA analyses of mRNAs, which show more than 1.5-fold elevated expression in highly metastatic SA, SPH, and METS cells compared with low-metastatic DLD1 cells in KEGG pathway gene sets. The three most significant gene sets are shown. B, GSEA analyses of mRNAs, which display significantly elevated expression in TCGA COAD patients with metastasis compared with patients without metastasis in KEGG pathway gene sets. The three most significant gene sets are shown. C and D, Association of the KEGG axon guidance signature with overall (C) and relapse-free (D) survival in the TCGA COAD dataset. E and F, Association of the KEGG axon guidance signature with metastasis (E) and nodal status (F) in the TCGA COAD dataset. G, qPCR analyses of PLXNB3 expression in DLD1 cell derivatives. H and I, Association of PLXNB3 expression with overall (H) and relapse-free (I) survival in the TCGA COAD dataset. J and K, Association of PLXNB3 expression with metastasis (J) and nodal status (K) in the TCGA COAD dataset. L and M, Enrichment of H3K4me3 (L) and H3K79me3 (M) marks in the PLXNB3 genomic region in DLD1 cell derivatives. N, DNA methylation of the PLXNB3 genomic region in DLD1 cell derivatives. O, DNA methylation of the PLXNB3 genomic region in TCGA COAD patients with metastasis (M1) compared with patients without metastasis (M0). Mean values ± SD are provided with *, P < 0.05; **, P < 0.01; ***, P < 0.001. In L–O, hg19 genomic coordinates and the PLXNB3 gene structure with indicated exons (white boxes, untranslated region; black boxes, protein encoding region) and transcription start site (arrow) is shown at the bottom.

Close modal

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.

1.
Tsai
JH
,
Donaher
JL
,
Murphy
DA
,
Chau
S
,
Yang
J
. 
Spatiotemporal regulation of epithelial-mesenchymal transition is essential for squamous cell carcinoma metastasis
.
Cancer Cell
2012
;
22
:
725
36
.
2.
Brabletz
T
,
Hlubek
F
,
Spaderna
S
,
Schmalhofer
O
,
Hiendlmeyer
E
,
Jung
A
, et al
Invasion and metastasis in colorectal cancer: epithelial-mesenchymal transition, mesenchymal–epithelial transition, stem cells and beta-catenin
.
Cells Tissues Organs
2005
;
179
:
56
65
.
3.
Ocana
OH
,
Corcoles
R
,
Fabra
A
,
Moreno-Bueno
G
,
Acloque
H
,
Vega
S
, et al
Metastatic colonization requires the repression of the epithelial-mesenchymal transition inducer Prrx1
.
Cancer Cell
2012
;
22
:
709
24
.
4.
Tsai
JH
,
Yang
J
. 
Epithelial-mesenchymal plasticity in carcinoma metastasis
.
Genes Dev
2013
;
27
:
2192
206
.
5.
Hay
ED
,
Zuk
A
. 
Transformations between epithelium and mesenchyme: normal, pathological, and experimentally induced
.
Am J Kidney Dis
1995
;
26
:
678
90
.
6.
Tam
WL
,
Weinberg
RA
. 
The epigenetics of epithelial-mesenchymal plasticity in cancer
.
Nat Med
2013
;
19
:
1438
49
.
7.
Vizoso
M
,
Esteller
M
. 
DNA methylation plasticity contributes to the natural history of metastasis
.
Cell Cycle
2015
;
14
:
2863
4
.
8.
Iliopoulos
D
,
Hirsch
HA
,
Struhl
K
. 
An epigenetic switch involving NF-kappaB, Lin28, Let-7 MicroRNA, and IL6 links inflammation to cell transformation
.
Cell
2009
;
139
:
693
706
.
9.
Rokavec
M
,
Oner
MG
,
Li
H
,
Jackstadt
R
,
Jiang
L
,
Lodygin
D
, et al
IL-6R/STAT3/miR-34a feedback loop promotes EMT-mediated colorectal cancer invasion and metastasis
.
J Clin Invest
2014
;
124
:
1853
67
.
10.
Thiery
JP
,
Acloque
H
,
Huang
RY
,
Nieto
MA
. 
Epithelial-mesenchymal transitions in development and disease
.
Cell
2009
;
139
:
871
90
.
11.
Mani
SA
,
Guo
W
,
Liao
MJ
,
Eaton
EN
,
Ayyanan
A
,
Zhou
AY
, et al
The epithelial-mesenchymal transition generates cells with properties of stem cells
.
Cell
2008
;
133
:
704
15
.
12.
Takeyama
Y
,
Sato
M
,
Horio
M
,
Hase
T
,
Yoshida
K
,
Yokoyama
T
, et al
Knockdown of ZEB1, a master epithelial-to-mesenchymal transition (EMT) gene, suppresses anchorage-independent cell growth of lung cancer cells
.
Cancer Lett
2010
;
296
:
216
24
.
13.
Harma
V
,
Virtanen
J
,
Makela
R
,
Happonen
A
,
Mpindi
JP
,
Knuuttila
M
, et al
A comprehensive panel of three-dimensional models for studies of prostate cancer growth, invasion and drug responses
.
PLoS One
2010
;
5
:
e10431
.
14.
Yu
M
,
Bardia
A
,
Wittner
BS
,
Stott
SL
,
Smas
ME
,
Ting
DT
, et al
Circulating breast tumor cells exhibit dynamic changes in epithelial and mesenchymal composition
.
Science
2013
;
339
:
580
4
.
15.
Kristensen
H
,
Haldrup
C
,
Strand
S
,
Mundbjerg
K
,
Mortensen
MM
,
Thorsen
K
, et al
Hypermethylation of the GABRE∼miR-452∼miR-224 promoter in prostate cancer predicts biochemical recurrence after radical prostatectomy
.
Clin Cancer Res
2014
;
20
:
2169
81
.
16.
Emmrich
S
,
Rasche
M
,
Schoning
J
,
Reimer
C
,
Keihani
S
,
Maroz
A
, et al
miR-99a/100∼125b tricistrons regulate hematopoietic stem and progenitor cell homeostasis by shifting the balance between TGFbeta and Wnt signaling
.
Genes Dev
2014
;
28
:
858
74
.
17.
Mazzoccoli
G
,
Pazienza
V
,
Panza
A
,
Valvano
MR
,
Benegiamo
G
,
Vinciguerra
M
, et al
ARNTL2 and SERPINE1: potential biomarkers for tumor aggressiveness in colorectal cancer
.
J Cancer Res Clin Oncol
2012
;
138
:
501
11
.
18.
Pavon
MA
,
Arroyo-Solera
I
,
Cespedes
MV
,
Casanova
I
,
Leon
X
,
Mangues
R
. 
uPA/uPAR and SERPINE1 in head and neck cancer: role in tumor resistance, metastasis, prognosis and therapy
.
Oncotarget
2016
;
7
:
57351
66
.
19.
Gupta
RA
,
Shah
N
,
Wang
KC
,
Kim
J
,
Horlings
HM
,
Wong
DJ
, et al
Long non-coding RNA HOTAIR reprograms chromatin state to promote cancer metastasis
.
Nature
2010
;
464
:
1071
6
.
20.
Dhamija
S
,
Diederichs
S
. 
From junk to master regulators of invasion: lncRNA functions in migration, EMT and metastasis
.
Int J Cancer
2016
;
139
:
269
80
.
21.
Kogo
R
,
Shimamura
T
,
Mimori
K
,
Kawahara
K
,
Imoto
S
,
Sudo
T
, et al
Long noncoding RNA HOTAIR regulates polycomb-dependent chromatin modification and is associated with poor prognosis in colorectal cancers
.
Cancer Res
2011
;
71
:
6320
6
.
22.
Takahashi
Y
,
Sawada
G
,
Kurashige
J
,
Uchi
R
,
Matsumura
T
,
Ueo
H
, et al
Amplification of PVT-1 is involved in poor prognosis via apoptosis inhibition in colorectal cancers
.
Br J Cancer
2014
;
110
:
164
71
.
23.
Stange
DE
,
Engel
F
,
Longerich
T
,
Koo
BK
,
Koch
M
,
Delhomme
N
, et al
Expression of an ASCL2 related stem cell signature and IGF2 in colorectal cancer liver metastases with 11p15.5 gain
.
Gut
2010
;
59
:
1236
44
.
24.
Rizzolio
S
,
Tamagnone
L
. 
Semaphorin signals on the road to cancer invasion and metastasis
.
Cell Adh Migr
2007
;
1
:
62
8
.
25.
Vitriol
EA
,
Zheng
JQ
. 
Growth cone travel in space and time: the cellular ensemble of cytoskeleton, adhesion, and membrane
.
Neuron
2012
;
73
:
1068
81
.
26.
Tsai
MC
,
Manor
O
,
Wan
Y
,
Mosammaparast
N
,
Wang
JK
,
Lan
F
, et al
Long noncoding RNA as modular scaffold of histone modification complexes
.
Science
2010
;
329
:
689
93
.
27.
Rokavec
M
,
Oner
MG
, 
Hermeking H. lnflammation-induced epigenetic switches in cancer
.
Cell Mol Life Sci
2016
;
73
:
23
39
.
28.
Ehrlich
M
. 
DNA hypomethylation in cancer cells
.
Epigenomics
2009
;
1
:
239
59
.
29.
Goelz
SE
,
Vogelstein
B
,
Hamilton
SR
,
Feinberg
AP
. 
Hypomethylation of DNA from benign and malignant human colon neoplasms
.
Science
1985
;
228
:
187
90
.
30.
Feinberg
AP
,
Vogelstein
B
. 
Hypomethylation distinguishes genes of some human cancers from their normal counterparts
.
Nature
1983
;
301
:
89
92
.
31.
Gama-Sosa
MA
,
Slagel
VA
,
Trewyn
RW
,
Oxenhandler
R
,
Kuo
KC
,
Gehrke
CW
, et al
The 5-methylcytosine content of DNA from human tumors
.
Nucleic Acids Res
1983
;
11
:
6883
94
.
32.
Hur
K
,
Cejas
P
,
Feliu
J
,
Moreno-Rubio
J
,
Burgos
E
,
Boland
CR
, et al
Hypomethylation of long interspersed nuclear element-1 (LINE-1) leads to activation of proto-oncogenes in human colorectal cancer metastasis
.
Gut
2014
;
63
:
635
46
.
33.
Matouk
I
,
Raveh
E
,
Ohana
P
,
Lail
RA
,
Gershtain
E
,
Gilon
M
, et al
The increasing complexity of the oncofetal h19 gene locus: functional dissection and therapeutic intervention
.
Int J Mol Sci
2013
;
14
:
4298
316
.
34.
Raveh
E
,
Matouk
IJ
,
Gilon
M
,
Hochberg
A
. 
The H19 Long non-coding RNA in cancer initiation, progression and metastasis—a proposed unifying theory
.
Mol Cancer
2015
;
14
:
184
.

Supplementary data