DNA methylation profile is remodeled in primary-resistant mCRCs.
Primary-resistant mCRCs were selected for this study as representative colorectal malignancies with poor prognosis and poor response to anti-cancer agents (15). Thus, in order to identify epigenetic alterations with prognostic relevance, global DNA methylation was assessed on 24 mCRCs primary-resistant to 1st-line FOLFOX (8 patients) or FOLFIRI (16 patients) chemotherapy combined or not with molecular targeted agents. Twelve drug-sensitive mCRCs (4 treated with FOLFOX and 8 with FOLFIRI combined with molecular targeted agents) were used as controls to obtain the differential methylation profile between primary-resistant and drug-sensitive tumors (in-house cohort; Additional Table 1). Differential methylation profiles were analyzed in a multistep process, as described in Additional Fig. 1. Indeed, 74843 and 36876 probes were significantly differentially methylated in FOLFOX and FOLFIRI datasets (p-value < 0.05), respectively (Fig. 1A-B) and these were widely distributed between different genomic regions (Fig. 1C-D).
Epigenetic alterations predict prognosis in human mCRCs.
Since it is well established that promoter hypo/hypermethylation is the main mark resulting in gene expression modifications (16), only genes with methylation modifications in promoter regions (with a p-value < 0.05 and an absolute difference of Beta value > 0.1) were used in subsequent analyses. Using this more restrictive cut-off, 7432 probes, corresponding to 4533 genes, for patients treated with 1st-line FOLFOX and 5005 probes, corresponding to 3803 genes, for patients treated with 1st-line FOLFIRI resulted differentially methylated between drug-resistant and drug-sensitive tumors. These two lists were named, respectively, FOLFOX_DMGs and FOLFIRI_DMGs (Additional Dataset). We next quested whether these DMGs were also fMET, with a methylation profile consistent with the gene expression profile. Since we could not obtain gene expression data from in-house colorectal tumor samples due to the poor amount and quality of RNA purified from paraffin-embedded specimens, this issue was addressed using a cohort of 33 mCRCs obtained from the TCGA COAD database, which provides gene expression, DNA methylation, DNA sequencing and clinical data for each patient (Additional Table 2). From the analysis of the TCGA COAD database, we obtained 741 fMET genes defined as COAD fMET genes. Among these 741 TCGA COAD fMET genes, 542 resulted DMGs in FOLFOX dataset and 248 in FOLFIRI dataset. Applying more stringent filters (p-value < 0.01 and absolute difference of Beta value > 0.2) on the FOLFOX_DMGs and the FOLFIRI_DMGs datasets and the COAD fMET genes (R2 > 0.25), this list was restricted to 8 fMET genes for FOLFOX and 7 fMET genes for FOLFIRI datasets, with 3 of them common to both datasets (Additional Table 4). Interestingly, hierarchical clustering on these sets of genes (using both expression and methylation data) allowed to separate the TCGA COAD samples into two quite homogeneous clusters characterized by over or under expression (Fig. 2A and D) and hypo or hyper methylation (Fig. 2B and E) of, respectively, the above 8 and 7 fMET genes. A similar separation was obtained in our in-house FOLFOX and FOLFIRI cohorts upon hierarchical clustering of methylation data using the same gene sets (Fig. 2C and F). In order to evaluate the prognostic relevance of these 8 and 7 fMET gene signatures, log-rank test was performed on both the TCGA COAD and the in-house cohorts using the two previously obtained clusters. Noteworthy, a significant (p-value < 0.05) separation for both Relapse Free Survival (RFS) and Overall Survival (OS) was observed between hypermethylated/underexpressed tumors, characterized by worst prognosis, and hypomethylated/overexpressed tumors using both FOLFOX (Additional Fig. 2A-D) and FOLFIRI (Additional Fig. 3A-D) 8 and 7 genes signatures. Log-rank test performed on in-house FOLFOX (Additional Fig. 4A-B) and FOLFIRI (Additional Fig. 4C-D) cohorts resulted in a significant (p-value < 0.05) separation of the two clusters considering RFS and a non-significant separation considering OS. Based on this evidence, the two cohorts were labeled as “good” and “poor” prognosis clusters.
As a next step, we combined the 8 and 7 fMET gene signatures from FOLFOX and FOLFIRI datasets obtaining a new signature of 12 fMET genes, being 3 of them common to both datasets (Fig. 3). Hierarchical clustering on such signature separated TCGA COAD patients into two well-defined cohorts (22 hypo and 11 hypermethylated tumors, Fig. 3A-B and Additional Table 2). A similar clustering was obtained in our in-house FOLFOX dataset (Fig. 3C) and partially in the FOLFIRI dataset (Fig. 3D). Furthermore, upon combination of in-house FOLFOX and FOLFIRI datasets, the 12-genes signature obtained a more defined separation between 6 hypo and 30 hypermethylated tumors (Fig. 3E). Log-rank test confirmed that the cohort with hypermethylation of the 12-genes signature is characterized by significantly shorter survival (both OS and RFS; p-value < 0.05) compared to the cohort with hypomethylation of these genes in both TCGA COAD and in-house cohorts (Fig. 3F-G). Altogether, these data suggest that this 12-fMET genes signature discriminates between mCRC patients with good (hypomethylated tumors) and poor (hypermethylated tumors) prognosis.
The poor prognosis hypermethylated cluster is characterized by a MSI-like phenotype and is enriched of CIMP-high tumors.
The poor prognosis hypermethylated and the good prognosis hypomethylated clusters were further characterized respect their clinical and biological profiles using gene expression and DNA sequencing and gene copy number data from the TCGA COAD database. No major differences were observed between the poor and good prognosis clusters respect to T and N categories and sites of primary tumor (right versus left colon) (Additional Table 2). Similarly, no major differences were observed respect to the tumor mutational load, with the exception of two hypermutated cases in the poor prognosis cluster (Additional Fig. 5, insert). Interestingly, specific gene mutations were differently distributed between the two subgroups, being mutations in SRGAP2B, AC007682.1, AC104820.2, AF121898.3 and MROH5 genes enriched in the good prognosis cluster and mutations in GRP98, NRXN2, HDN1 and TTC40 genes in the poor prognosis cluster (Additional Fig. 5). Consistently, several gene aberrations were significantly more abundant in the good prognosis cluster (Additional Fig. 6).
A differential gene expression comparison between the two clusters yielded 116 Differentially Expressed Genes (DEGs) (false discovery rate, FDR, adjusted p-value < 0.05 and abs(logFC) > 0.58). Among these, 80 genes were downregulated and 36 of them hypermethylated in the poor prognosis cluster and conversely upregulated/hypomethylated in the good prognosis subgroup (Additional Fig. 7). Gene Set Analysis (GSA) was performed on the gene set collection of the mSigDB repository, obtaining significant enrichments for signaling pathways and positional collections. Among different signaling pathways (Fig. 4A), GSA identified the Watanabe gene set, which includes genes discriminating between MSI and MSS (microsatellite instability/stability) colorectal cancers (17). The statistical analysis identified 7 genes in our list of DEGs, which enrich Watanabe gene dataset and whose expression profile is consistent with a separation of the TCGA cohort in good and poor prognosis clusters (Fig. 4B). Noteworthy, 5/7 genes enriching Watanabe gene set were significantly downregulated due to hypermethylation. Since this observation suggests that the 12-genes hypermethylated signature identifies a subgroup of mCRCs with a MSI-like phenotype, an independent MSI-like gene expression signature was evaluated for the capacity to reproduce the separation of the TCGA COAD cohort in the same good and poor prognosis clusters, according to the 12-genes signature. Noteworthy, the MSI-like gene expression signature of Pačínková et al (18) mirrored the separation of the 33 mCRCs TCGA cohort in the same clusters as obtained by our 12-genes signature (Fig. 4C). Consistently, five genes from the Pačínková signature (i.e., PLAGL2, ACSL6, ARID3A, NKD1, TNNT1) were characterized by an expression profile consistent with the expression profile of the MSI-like poor prognosis TCGA cluster.
Interestingly, 15–20% of human CRCs are characterized by the CpG-island methylator phenotype (CIMP), subdivided in CIMP-high (CIMP-H) and CIMP-low (CIMP-L) and this correlates with the MSI phenotype (19). Thus, the relationship between our 12-genes hypermethylated signature and CIMP status was evaluated in TCGA COAD dataset and in our in-house cohort according to Hinoue et al (20). Interestingly, the poor prognosis TCGA hypermethylated cohort was enriched of CIMP-H cases, being 5 out of 6 CIMP-H samples classified as poor prognosis patients, whereas the hypomethylated good prognosis cohort was enriched on no-CIMP tumors. CIMP-L cases were distributed between the two subgroups, being classified in 7 out of 11 cases as good prognosis patients (Fig. 4D). This difference between the groups was statistically significant by two-sided Fisher exact test (p-value < 1e-02). Consistently, the vast majority of CIMP-H tumors were characterized by the hypermethylation of the 12-genes signature in our in-house cohort (Additional Fig. 8A). Thus, the prognostic relevance of our 12-genes signature was compared to CIMP status and, noteworthy, the RFS and the OS analyses showed a better capacity of our 12-genes hypermethylated signature to predict poor prognosis in both the TCGA COAD (Fig. 4E-F) and the in-house cohorts (Additional Fig. 8B-C). Altogether, these observations strongly support the conclusion that the 12-genes hypermethylated signature correlates with a MSI-like phenotype and is characterized by a better capacity to predict prognosis compared to CIMP status.
Hypermethylated genes are enriched on arms q of chromosomes 13 and 20.
Among positional collections, GSA identified enrichments of chromosome 13 arm q and chromosome 20 arm q gene sets (Fig. 5A). Interestingly, 23/80 downregulated genes in our list of DEGs are located on chromosome 13 arm q and 8 of them are hypermethylated. Consistently, 12/80 downregulated genes are located on chromosome 20 arm q. Noteworthy, the expression profile of each of these gene sets (Additional Fig. 9A-B) and of their combination (Fig. 5B) mirrored the separation of the TCGA COAD cohort in the good and poor prognosis clusters, obtained according to the 12-genes signature. These data suggest an enrichment of methylation events in genes located in specific chromosomal regions in mCRCs with poor prognosis. In such a context, the enrichment analysis for signaling pathways identified the Ding lung cancer expression by copy number and the Nikolsky breast cancer 20q12-13 amplicon gene sets (Fig. 4A). These authors reported respectively a correlation between the copy number variation and the expression of 26 genes in lung cancers (21) and the identification of 149 genes in amplicon 20q12-13 in breast tumors (22). As expected, all the genes enriching the Nikolsky gene set overlap with our chromosome 20 arm q genes and 5 out of 6 Ding genes with our chromosome 13 arm q genes. In both cases, the expression profile of these genes reproduced the separation of the TCGA dataset in good and poor prognosis clusters (Additional Fig. 9C-D). Altogether, these data highlight the relevance of expression/methylation modifications of genes located on chromosomes 13 and 20 in human colorectal cancer.
Epigenetic modifications are reproduced in drug-resistant cell models.
To validate epigenetic data obtained from FOLFOX and FOLFIRI primary-resistant mCRCs, we generated in vitro drug-resistant cellular models chronically adapted to oxaliplatin (Oxa; HCT116-OxaR and HT29-OxaR cells) or irinotecan (Iri; HCT116-IriR and HT29-IriR cells). In preliminary experiments, apoptotic cell death was evaluated in drug-sensitive and drug-resistant cell lines in response to Oxa or Iri in combination or not with the demethylating agent 5-Aza-dC. These experiments confirmed that cell lines chronically exposed to chemotherapeutics are indeed poorly sensitive to Oxa and Iri and that drug resistance is reverted upon treatment with 5-Aza-dC (Fig. 6A-B). Since these data support the hypothesis that methylation modifications are responsible for resistance to Oxa and Iri in these CRC cell lines, in subsequent experiments, drug-resistant cell lines were used to validate the expression profiles of the 12-fMET genes signature. Real Time RT-PCR analysis of the 12-genes confirmed a significant downregulation of 10 genes in HCT116-OxaR and HT29-OxaR cell lines (Fig. 6C-D) and 9 and 7 genes in, respectively, HCT116-IriR and HT29-IriR cell lines (Fig. 6E-F) compared to the respective drug-sensitive cell lines. NROB2 was undetectable in both cell lines (data not shown). Noteworthy, the pretreatment of drug-resistant cell lines by 5-Aza-dC resulted in a significant upregulation of the majority of the downregulated genes: 7 genes in, respectively, HCT116-OxaR and HT29-OxaR cells (Fig. 6C-D) and 6 genes in, respectively, HCT116-IriR and HT29-IriR cell lines (Fig. 6E-F). These data suggest that the downregulation of genes belonging to the 12-genes signature correlates with onset of drug resistance in CRC cell lines and that this downregulation is likely mediated by methylation events.
In parallel analyses, starting from GSA results suggesting an enrichment of a MSI-like phenotype in hypermethylated poor-prognosis tumors, 7 genes enriching Watanabe pathway and 5 representative genes belonging to the mismatch repair system (18) were evaluated in drug-resistant cell lines (Additional Fig. 10). Indeed, PCR data confirmed the downregulation of the majority of these genes in drug-resistant cell lines with HT29-OxaR and HT29-IriR characterized by a more evident MSI-like phenotype (Additional Fig. 10).