Hypermethylated Colorectal Cancer Patients Present a Myc-Driven Hypermetabolism with a One-Carbon Signature Associated to Worsen Prognosis

1. Introduction

Colorectal cancer (CRC) is the second cause of cancer-related death worldwide, responsible for nearly a million death in 2020 [1]. CRC tumours can arise from two main mechanisms of carcinogenesis: the conventional pathway, associated with mutations in APC and p53, and the serrated pathway that is associated with mutations in KRAS or BRAF and a hypermethylated phenotype (CpG island methylator phenotype or CIMP). Serrated tumours, although less frequent, have been shown to be the most aggressive subtype of CRC [2]. Nevertheless, the strong heterogeneity in their pathological presentation and mutational pattern makes them very hard to diagnose. Moreover, multiple methylation panels are used to diagnose the CIMP phenotype, leading to heterogeneity in the diagnosis and new means of identifying these higher risk tumours are needed.

The CIMP subtype of CRC represent around 20% of all CRC patients and are the ones with the worst survival and resistance to classical treatment like 5-Fluorouracil [3]. Moreover, several meta-analysis have highlighted the disparities in the diagnosis of CIMP patients: there is no consensus on the gene panels to be used or in the experimental method of choice to test the methylation of these genes leading to 16 different definitions of CIMP in the literature [4]. There is hence a need to homogenize our definition of CIMPs and find reliable phenotypic markers for this group of patients.

Since KRAS and BRAF are known regulators of metabolism in multiple cancer, driving its rewiring in tumour cells, CIMP tumours could harbor a particular metabolic signature differentiating them from other types of CRCs. Indeed, it has previously been shown that CIMP have upregulation of enzymes in glycolysis and related processes [5]. Moreover, since CIMP tumours are characterized by a hypermethylation phenotype, it is very likely that the metabolism of these tumours also needs to adapt to accommodate for the higher demand for methyl groups.

Methylation is one of the main epigenetic ways to regulate gene expression: it consists on the addition of a methyl group to the cysteine in position 5 of histones and usually results in the downregulation of the methylated gene. The one-carbon (1-C) metabolism is the central provider of methyl groups needed for methylation processes and has been shown to be associated with poorer prognosis in colorectal and other cancers [6,7,8,9,10,11]. Indeed, 5-Fluorouracil, an anti-metabolite chemotherapy targeting 1-C metabolism, is the first line of treatment for CRC patients [12] and one of the most widely used chemotherapy [13], indicating the importance of this pathway for cancer cell growth in many different types of cancers. Finally, it has recently been shown that the mitochondrial part of the 1-C metabolism is not required for cancer growth but fosters migration of cancer cells [14]. Altogether these data suggest that 1-C metabolism is important for both cancer growth and dissemination. Nevertheless, its contribution to CIMP specifically and its value as a prognostic factor have never been investigated. Furthermore, targeting of the 1-C metabolism pathway is difficult because it is so essential to normal cells and hence generates a lot of resistance to treatment. To date, most of the targeting is done through targeting the Dihydrofolate reductase (DHFR) or the Thymidylate synthase (TYMS) but more targets for therapeutic intervention are required.

The aim of the present study is to evaluate if CIMP tumours present a distinct metabolic signature, if their levels of methylation could be predicted using transcriptomic data and if new therapeutic targets can be identified.

2. Materials and Methods

2.1. TCGA RNA-Sequencing Dataset

TCGA consortium RNA-sequencing data matrices of Z-scores from the CRC-2012 cohort [15] were downloaded with their corresponding clinical data on Cbioportal website [16]. Supplemental biological information was also provided for these patients such as: methylation group status, MLH1 status, hyper-mutation and MSI status. This cohort is composed of 223 patients and stratification on methylation status showed 80 patients with positive methylation status (CIMP low and CIMP high) versus 143 patients with methylation status negative (cluster3 and cluster4) (Table 1).

The main clinical parameters provided in this study such as primary site of tumour localization, oncotree code, cancer type were found significant on this methylation stratification (Table 1). Metabolism pathway enrichment was performed with Geneset enrichment Analysis software version 4.0.1 [17] with KEGG and REACTOME databases implemented in MSigDB version 7.0 [18]. Bioinformatics analyses were performed in software environment version 3.5.3. Unsupervised principal component analysis was performed with Factominer R-package [19]. Boxplot were drawn with ggplot2 graph definition [20] and expression heatmap were done with pheatmap R-package with the option of Euclidean distances. Machine learning random forest was performed on selected transcriptome features during unsupervised analysis. This learning was perfomed with Randomforest R-package [21] and the model was built with 150 trees and optimized mtry parameter after tuning with rftune function. Variables importance were kept in the model to interpret gene priority on their accuracy.

2.2. Overall Survival Analysis on TCGA RNA-Sequencing with Multi-Omics Integration in Colorectal Cancer

Using the Cbioportal web application [16], RNA-seq data from the TCGA 2018 colorectal cancer cohort (our validation cohort) was analyzed with Z-score diploid V2 matrix normalization. This dataset composed of 379 patients also had overall survival information [22]. This interface was used to investigate one carbon metabolism gene candidates, defined by AMIGO web database of gene ontology information [23] and available in Supplementary Table S3, on patient survival.

2.3. ChIP-Sequencing Analysis

Results from ChIP sequencing experiments performed in the colorectal cell line LoVo [24] were downloaded on the Cistrome project website. Gene promoter annotation of MYC epigenetics intervals were performed with BETA minus algorithm version 1.0.7 on human genome version HG38. Proximal peaks were filtrated around Transcription starting sites with 5kb upstream and 0.5kb downstream. Functional enrichment on MYC epigenetics genomic intervals (HG38) was predicted with GREAT website application with Gene Ontology Biological Process database [25]. Functional enrichment network of GREAT metabolism connections was drawn with Cytoscape software version 3.6.0 [26]. Motif prediction was performed with RSAT genomic application on human HG38 with Jaspar vertebrae 2018 motif collection [27]. Promoter heatmap were performed with deeptools version 3.3.1 in Linux MINT 19 operating system with custom wrapper in BASH script (https://github.com/cdesterke/chip2heat). Full results are available on Supplementary Table S2.

2.4. Integrative Analysis

In order to integrate Myc epigenetics analyses with RNA-sequencing transcriptome of CRC patients, genomic coordinates of CHIP-sequencing initially aligned on HG38 human transcriptome were liftover to HG19 genomic coordinates. This operation was done with liftover algorithm on UCSC website (https://genome.ucsc.edu/cgi-bin/hgLiftOver). Circosplot was drawn with OmicCircos R bioconductor package for multi-omics integration [28].

2.5. Deep Learning

Keras version 2.7.0, Scikit learn version 1.0.1, tensorflow 2.7.0, pandas 1.3.5 and matplotlib version 3.5.1 python libraries (python version 3.7.6) were used with Jupyter notebook version 6.4.6 to implement deep-learning neuronal network for validations of MYC targets predicting clinical data: methylation status, hyper-mutated phenotype, MLH1 silencing and MSI status. Neuronal network was built with three sequential neuron dense layers: first one with 12 neurons and “Rectified Linear Unit” (ReLU) activation function, second one with 8 neurons with ReLU activation function and third one with one neuron and sigmoid activation function. The model was compiled with “adam” optimizer and loss function based on binary cross entropy and fit with 150 epochs and 10 batch size. Corresponding python code was deployed at the following web address: https://condescending-saha-2dfa16.netlify.app/ (accessed on 2021, December 20th).

2.6. Multivariable Model Built on Methylation Status Outcome

In order to test dependency between qualitative variables of the study “chi2loop” R-package was developed to perform iteration chi.square test between character variables from the cohort dataset. This package is available at the following address: https://github.com/cdesterke/chi2loop (accessed on 2021, nov 30th). This package takes as input a dataset of qualitative variables. To run iterative chi.square tests with “cltest” function, the dataset need to be imported with parameter “stringsAsFactors = FALSE” because input qualitative variables need to be in character format and not as factors. On dataset “cltest” function could be applied to perform chi.square iteration between character columns and output results could be graphically represented with the “nlpplot” function (NLP: negative log10 of chi.squre test p-values). Subsequently “chinet” function could be applied also on “cltest” results to detect variable communities with a Louvain classification algorithm. After selection of associated clinical and expression markers, a multivariable model was built with generalized linear model (GLM) R base function for a logistic regression using binomial family as parameters (methylation status as binary outcome).

Graphical output of multivariable logistic regression model censored on binomial status of methylation was drawn with the R-packages combination: broom version 0.7.10, broom.helpers version 1.4.0 and GGally version 2.1.2. In the multivariable logistic model, a p-value p<0.0.5 for included variables in the model was taken as threshold of significance for their independency.

2.7. Statistical Analyses

Statistical analyses were performed in R software environment version 4.1.0. With Z-scores from RNA-seq data (Gaussian transformation), multi-group comparison was performed with one Way ANOVA Fisher statistical test and dual group comparison was performed with bilateral ttest with Welch correction. For both approach, statistical significance was assessed for a p-value inferior to 0.05.

3. Results

3.1. Hypermetabolism in CIMP CRC Transcriptome

To investigate the metabolic signature of CIMP tumours, we used the CRC-2012 TCGA cohort [15], comprising 222 patients, where methylation levels were measured and CIMP statuses determined. Patients were hence divided into unmethylated (clusters C3 and C4) versus CIMP patients (clusters CIMP_low and CIMP_high). REACTOME and KEGG databases were used to perform Geneset enrichment analysis (GSEA) to compare the transcriptome of methylated and unmethylated patients, with an emphasis on metabolic pathways. Surprisingly, we found that all metabolic pathway significantly modulated in CIMP patients were up-regulated. This hypermetabolism associated to the CIMP phenotype comprised 169 enzymes (Supplementary Table S1), involved in pathways of the biochemistry of carbohydrates, nucleotides and sphingolipids (Figure 1A). We then performed an unsupervised principal component analysis (PCA) using this metabolic signature that confirmed its ability to separate methylated and unmethylated patients (p-value = 3.35e-19, Figure 1B and Supplementary Table S2) but was also able to show a progressive stratification of methylation axis of the unsupervised analysis (p-value = 6.48e-25, Figure 1C). These results suggest a close association between the CIMP phenotype and this hypermetabolism in CRC patients. Indeed, feeding our 169 metabolic markers to a supervised machine learning algorithm was able to discriminate methylated from unmethylated tumors with a global efficiency of 77.58% after 150 trees of random forest learning (Figure 1D). This efficiency was confirmed by unsupervised classification, performed with Euclidean distances, showing on the expression heatmap a smaller cluster comprising a majority of CIMP tumor and a larger cluster comprising a majority of unmethylated patients (Figure 1E). A random Forest variable importance study further confirmed the importance of glycolysis, nucleotide- and one carbon- metabolism to predict CIMP status, with GAPDH being the best metabolic marker (Figure 1F, Supplementary Table S2).

3.2. Myc Regulates one Third of the CIMP-CRC Metabolic Program

We next wondered if there could be a master regulator orchestrating this hypermetabolism, characteristic of CIMP tumours. Myc is a strong promoter of metabolism to foster cell proliferation and cell fitness [29] but also to control stem cell fate decisions [30]. It was shown that its expression is changed in nearly all colorectal tumours, suggesting a strong role in CRC pathogenesis [15]. We asked if the metabolic rewiring observed in CIMP could be driven by Myc by using publicly available ChIP-seq data from the colorectal cell line LoVo [24]. Promoter heatmap of the 38663 MYC peaks confirmed that the signal is well centered on Transcription Starting sites (TSS) (Figure 2A). A phast-conservation study for Myc confirmed that its chromatin binding events were found well conserved on the core mammalian promoter database (Figure 2B). Myc chromatin bindings were then mapped on HG38 human genome promoters and proximal events were filtered around the TSS (5kb upstream and 0.5kb downstream). The resulting 5000 gene promoter prediction was crossed, in the LoVo cell line, with genes from the CIMP metabolic signature and we found 53 genes with Myc binding sites in their promoters (Figure 2C, Supplementary Table S2). These results indicate that at least a third of the metabolic signature observed in CIMP patients could be driven by Myc. A genomic Circosplot on the Myc signature revealed that sexual chromosomes, chromosomes 6, 13, 16 and 20 were not involved in the CIMP-specific metabolic program (Figure 2D). Chromosome 12 contained the most gene from the signature with eight Myc targets, especially GAPDH that we found as best predictive marker of the CIMP phenotype (Figures 1F and 2D).

3.3. Genes from the Myc Transcriptional Program Also Have Binding Sites for Other Transcription Factors

The observed mean size of Myc chromatin binding peaks in the CIMP-specific metabolic program was 651 pb (SD: 259 pb, n=54 peaks for 53 genes) and the identified peaks were narrow and centered on Transcription Starting Sites (TSS) (Figure 3A). This proximal transcriptional program was compatible to perform motif prediction with the RSAT web application, using the JASPAR core promoter motif database. After picking corresponding nucleotide sequences on HG38 human genome, we performed a transcription factor binding motif identification and confirmed MAX binding sites, amongst other Myc known co-factors (Figure 3B). Interestingly, we also found 182 compatible sites for homeobox transcription factors like POU6F1, DLX1 and EMX2 – 142, 156 binding sites for Kruppel like factors KLF16 and KLF5 respectively and finally 111, 73 and 78 binding sites for HOX family homeobox transcription factors HOXA13, HOXA1913, HOXB13 (Figure 3B). Myc metabolic chromatin binding marks and found activated in transcriptome of CIMP-CRC, were conjointly analyzed with LOVO CHIP sequencing of DLX1, KLF5 and HOXA13. The resulting proximal promoter heatmap analyses revealed that myc metabolic program found CIMP CRC shared chromatin proximal binding sites with these transcription factors previously predict by RSAT: KLF5, DLX1 and HOXA13 (Figure 3C). Also, at transcriptome level, MAX, DLX1 and SP8 transcription factors were confirmed as being progressively overexpressed with level of methylation. These results suggested that myc metabolic program found CIMP CRC could be active in a chromatin context implicating homebox and Kruppel like factors.

Functional enrichment, performed on genomic intervals from the Myc signature on Gene Ontology Biological Process database (GO-BP), showed that the principally altered pathways were the carboxylic acid, tetrahydrofolate interconversion, nucleobase containing small molecule and oxoacid pathways (Figure 4A). This metabolism enrichment highlighted enzymes at the interface between glycolysis (GAPDH, ENO1, LDHA) and one-carbon metabolism (TYMS, SHMT2, MTHFR, MTHFD1), further confirming the importance of these two pathways in CIMP tumours (Figure 4B).

3.4. Metabolism Targets in the Myc Signature Are Associated to Worst Clinical Group in CRC

Several clinical parameters are available in the CRC-2012 cohort: methylation status, expression subtype (clustering into 3 subtypes, according to their mRNA profile), hypermutation and MLH1 silencing. We next evaluated the correlation between these parameters and our CIMP-specific metabolic signature. Based on the classification from the Random Forest prediction factor analysis (Figure 1F), the Myc epigenetic score from the LoVo ChIP-seq (Supplementary Table S2) and also the methylation predictive score, the best candidates to discriminate the pathological parameters from our TCGA cohort were GAPDH, TYMS, DHFR and TK1. All these markers had a specific increased expression in CIMP (Figure 4C, first two left panels) and were correlated to hypermutation and MLH1 silencing (Figure 4C, last two right panels). To go further, we built a multivariate deep-learning model by neural network using keras tensorflow Python libraries and these four Myc-specific targets. We found that this model was very efficient in predicting the Hypermutation and MLH1 silencing status of tumours (95 and 96% respectively) but also the CIMP and MSI status (82 and 85% respectively, Table 2). Predicting the MSI and MLH1 status is of particular interest, since it is one of the only parameters used to stratify patient care and decide on the use of immunotherapies [31]. Very interestingly, three of our four markers (TYMS, DHFR and TK1) are also part of the one-carbon metabolism pathway.

3.5. Overexpression of One Carbon Metabolism Enzymes Are Independent Markers of Methylation Status, MLH1 Silencing, Hypermutation and MSI in Colorectal Cancer

To study in more detail the prognostic value of one carbon metabolism in CRC, we selected the enzymes significantly associated to the methylation status (TYMS, TK1, SHMT2, MTHFD1, MTHFD2, DHFR, Figure 4B) and their normalized Z-score were transformed in quantile five (Q5) classes. This allowed us to perform iterative Chi-square tests with the different clinical parameters (for details see Supplementary Table S1) inputted as character string for Natural Language Processing (NLP) (see Materials and methods for more details, Figure 5A). Amongst one carbon metabolism genes, the expression of TYMS had the highest association to the hypermutation status, iCluster classification, methylation subtypes, MLH1 silencing and MSI status (Figure 5A). Consistent with these proteins being part of the same pathway, the same associations, although with less significance, were found for TK1, SHMT2 and MTHFD2 (Figure 5A). MTHDF1 did not show significant association with MSI status and was not able to discriminate between CIN/CIMP (Expression_subtype variable on Figure 5A). DHFR showed the weakest association with clinical parameters and was not associated with methylation or iCluster variables (Figure 5A). Interestingly, none of the one-carbon metabolism genes were able to discriminate for tumour stage (“Tumor stage” variable). This was consistent with the fact that stratifying patients per methylation status gave similar distribution of tumour stages (Table 1). We next performed a variable community detection using a Louvain algorithm on the results from the Chi.square tests (Figure 5A) This analysis confirmed the association of TYMS, TK1, SHMT2, MTHFD2 and DHFR expressions with the main prognosis markers (hypermutation, MSI, iCluster, expression subtype and tumor stage) while MTHFD1, and not MTHFD2, was more strongly associated to the methylation status (Figure 5B). Based on these observed associations, we finally built a methylation status multivariable model with minimum confounding parameters based on the expression of these one-carbon metabolism genes and clinical parameters. This logistic multivariate model particularly confirmed that overexpression of TYMS and MTHFD1 are strong, independent, one-carbon metabolism markers able to discriminate for the CIMP phenotype (Figure 5C). It also confirmed that CIMP tumours are strongly associated to MSI, as previously published [15].

3.6. Activation of 1-C Metabolism Genes Predicts Worst Prognosis Colorectal Cancer Patients

In order to validate the association of one-carbon metabolism and the hypermethylation status as well as bad prognosis for colorectal cancer patients, a multi-omics investigation was performed in the 2018 augmented CRC dataset (TCGA-COAD, our validation cohort) from the TCGA consortium [22]. Using the Gene Ontology database, through the AMIGO website, genes from the one-carbon metabolism were selected (Supplementary Table S3). We looked at the overexpression of these one-carbon metabolism genes in the CRC tumour RNAseq data and selected the most changed to establish a gene signature (Figure 6A).

Overall survival analysis revealed that patients presenting alteration in this one-carbon signature (i.e., increased expression of the genes in the signature) had worsen survival (Figure 6B). As previously reported [2], tumours from the right colon (comprising Ascending, Transverse colon, hepatic flexure and cecum) showed a significant association with increased one-carbon metabolism (Figure 6C) and were associated to a higher age at diagnosis (data not shown). Looking at HM450 methylation profiles, we found that the presence of this altered one-carbon signature was associated to an increase in the global methylation of the genome (Figure 6D and Supplemental Table S5 for a detailed list of genes with differentially methylated promoters).

We then performed a functional enrichment analysis on the differentially methylated genes using the KEGG or Gene Ontology database (Figure 6E,F). This revealed that most changes were seen in primary bile acid biosynthesis, but also affected vitamin absorption, PPAR signaling as well as stem cell functions (WNT signaling, pathways of pluripotency regulation). We also used the Gene Ontology database to perform the functional enrichment and found the most changes in gene implicated in the regulation of SMAD protein phosphorylation.

Altogether these results confirmed that the activation of one-carbon metabolism in CRC is associated to a worst prognosis for patients and a hypermethylation of gene promoters that could impact important functionalities in intestinal cells and have causative effects in cancer progression and dissemination. This strongly suggests that targeting one-carbon metabolism is of high interest for CRC patients.

4. Discussion

Despite the role of KRAS and BRAF in driving the CIMP phenotype and rewiring cancer metabolism [2,32,33,34], few studies have aimed at establishing if CIMP tumours have a different metabolism [5]. Using transcriptomics data and focusing only on metabolic pathways, we have found that CIMP tumours indeed present a hypermetabolism, with surprisingly no down-regulated but several up-regulated metabolic pathways. We were also able to determine a metabolic signature specific to CIMP tumours, able to discriminate hypermethylated tumours and even separate CIMP-L from CIMP-H patients with a high accuracy. This suggests the possibility of identifying these higher risk tumours, with high accuracy, by using transcriptomics data. This provides an alternative to the current techniques using the methylation of a panel of genes, with several gene panels being used over the world and no consensus and which panel and which technique to test for methylation [4]. Furthermore, we were also able to show that one-carbon metabolism enzymes were amongst the most efficient predictors of the CIMP phenotype, confirming the importance of this pathway, most likely in providing methyl groups to support the high methylation of the genome. Indeed, a background of high genome methylation dramatically fosters the ability of mutant BRAF to drive carcinogenesis, since high expression of mutated BRAF on its own rather triggers senescence and death [35]. This is what is thought to foster the development of BRAF-driven CIMP tumours [36].

Myc’s expression is changed in almost all colorectal cancer, underlying its role in tumorigenesis [15]. Amongst the genes in our metabolic signature, one third could be under the control of Myc in colorectal patients. Indeed we confirmed using ChIP-seq data that Myc physically binds to the promoter of these genes in the LoVo colorectal cancer cell line. Pardoxically, Myc expression is lowest in CIMP tumours, compared to the other tumour subtypes, where the signature is observed. We also found that binding sites for Max, one of the best characterized co-factors of Myc, are present within the promoters of the metabolic signature and that Max physically binds to the promoter of these genes in the LoVo Chip-seq data. This suggests that the specificity of the signature could be coming from the expression of co-factors rather than of Myc itself. Such a phenomenon has already been observed in Small Cell Lung Cancer, where it has been shown that Max expression regulates the 1C-metabolism pathway in a context-dependent manner [37].

Thanks to several meta-analysis of CRC cohorts, it is now becoming clear that the CIMP phenotype is associated to a worsen prognosis and that the benefit of an adjuvant fluorouracil (5-FU, one of the classical chemotherapy protocol for CRC patients) treatment after surgery is limited for CIMP-High patients, especially those with advanced stage III-IV tumours [3,38,39]. Given the important toxicity of these treatments, identifying those patients is very important to improve their quality of life and prevent their exposure to a medication with limited interest and high toxicity for them. Counter-intuitively we found that one of the gene most associated and predictive of the CIMP phenotype is TYMS (thymidylate synthase), which has been shown to be up-regulated by Myc and to drive sensitivity to 5-FU treatments [40]. This means that other mechanisms are at play and suggests that other targets in the 1C-metabolism pathway should be explored. Indeed, our results suggest that MTHFD1, 2 could be very interesting targets for the treatment of colorectal cancer. MTHFD2 is of particular interest since it has been shown to be expressed only in embryonic tissues but becomes re-activated in cancer tissues [41,42]. It is hence a tumour-specific marker that has been shown to be overexpressed and associated with a worse prognosis in many cancer tissues [42,43,44]. It has also been recently proposed to play a role in the regulation of the immune system, fostering cancer immune evasion [45]. Interestingly, in this study the authors showed that MTHFD2 promotes O-glycosylation which increases Myc’s stability and PD-L1 transcription. A recent study has shown that inhibitors of MTHFD2 lead to the specific death of MTHFD2-expressing cells, suggesting that these therapeutic approaches could be very specific to cancer cells and have minimal impact on normal cells.

5. Conclusions

Our study shows that transcriptomics holds great power towards unbiasedly identifying hypermethylated tumours. This approach could provide an alternative to the actual disparities in the methods used to assess methylation and allow more consistency and robustness. Additionally, we show that transcriptomics identified tumours with high MLH1 silencing, which is an indication of good responders to immunotherapy and these patients should receive this treatment rather than 5-FU [31]. This is in accordance with recent studies showing that transcriptomics is becoming more and more relevant in the clinical practice and provides a great clinical benefit to orient treatment protocols, especially for patients without tractable DNA mutations in their tumours [46].