Profile and diagnostic value of N6-methyladenosine-related mRNA regulators in colorectal adenoma and colorectal cancer.

Introduction
Colorectal cancer (CRC) is one of the most common malignancies of the digestive tract and is the second leading cause of cancer death worldwide [1]. Approximately 85% of CRC cases originate from colorectal adenomas (CRAs), 80% of which have adenomatous polyposis coli (APC) mutations and multiple gene mutations, such as those in KRAS, p53, and SMAD4; on the basis of these findings, these patients could further develop cancer [2]. It usually takes 5–15 years for CRAs to develop into carcinomas [3]. These 5–15 years were also the best window periods for the treatment of CRA and clinical prevention of CRC. Therefore, preventing CRA and blocking its progression to adenoma in a timely and effective manner are essential for reducing the incidence of CRC.
CRA can be classified as low-risk and high-risk adenomas based on their malignant potential, with high-risk adenomas including more than two villous or tubulovillous adenomas (with or without dysplasia), no less than 10 mm in diameter and considered precancerous lesions of CRC, while low-risk adenomas including 1 to 2 tubular adenomas less than 10 mm in diameter [4]. Epidemiological studies have shown that approximately 10% of adenomatous polyps develop into CRC, and approximately 25% of high-grade adenomas develop into colorectal cancer [5]. Colonoscopy screening and endoscopic adenoma eradication are currently the most effective methods for detecting and treating adenomas, but adenomas tend to recur after resection. Lengthy bowel preparation, unbearable colonoscopy pain, intraoperative perforation, and the surgical risk of bleeding and infection of postoperative wounds greatly reduce the motivation of patients to be examined and treated. Therefore, identification of key genes, exploration of the potential diagnostic value and prognostic effect of these genes in CRA, and prevention of progression from CRA to CRC are urgently needed.
Epigenetic modifications, such as posttranscriptional RNA and DNA methylation, are thought to be key regulatory mechanisms in many biological processes [6]. N6-methyladenosine (m6A) RNA methylation is an important epigenetic modification and is ubiquitous in eukaryotic cells. Additionally, m6A RNA modifications regulate RNA biogenesis, degradation, trafficking and cellular localization [7]. In recent years, m6A and its related regulators have been shown to be involved in the occurrence and development of tumors. For example, in bladder cancer, the m6A regulator METTL3 promoted bladder cancer progression through the AFF4/NF-κB/MYC signaling network [8]. In terms of chemocarcinogenicity, the METTL3-m6A-CDCP1 axis has synergistic effects on promoting the malignant transformation of urothelial cells and the occurrence of bladder cancer both in vitro and in vivo [9]. These studies fully support the present hypothesis that m6A regulators might serve as potential biomarkers for the early diagnosis of CRA or CRC. However, comprehensive analysis of the expression of m6A regulators in the mucosa-adenoma-adenocarcinoma process is lacking.
In this study, we aimed to explore the role of the m6A regulatory mechanism and disease diagnostic value in CRA and CRC patients based on mRNA expression profiles from the GSE41657, GSE100179, and GSE117606 datasets of the Gene Expression Omnibus (GEO). First, we analyzed the differentially expressed m6A (DE-m6A) regulators in CRA and CRC patients and divided CRA and CRC patients into two subgroups by consensus clustering and analyzed the differentially expressed genes (DEGs) between the two subgroups. Enrichment analysis was subsequently performed on the DEGs to explore m6A-related functions and pathways, and protein‒protein interaction (PPI) network analysis was subsequently performed on the DEGs to screen for hub genes related to m6A. Finally, important DE-m6A regulators were used to identify and verify CRA and CRC diagnostic markers.

Materials and methods

Data collection and processing
We downloaded the microarray data plus clinical information from the GSE117606, GSE41657, and GSE100179 datasets of the GEO database. In total, twenty m6A regulators were extracted from the dataset, including two erasers (ALKBH5 and FTO), five writers (METTL3, METTL14, WTAP, ZC3H13, and KIAA1429) and thirteen readers (HNRNPA2B1, FMR1, IGF2BP3, IGF2BP1, IGF2BP2, LRPPRC, RBMX, YTHDC1, YTHDC2, YTHDF1, YTHDF3, YTHDF2, and HNRNPC).

Consensus clustering analysis
To identify different m6A patterns, we performed a consensus clustering method using the “ConsensusClusterPlus” package in R software based on DE-m6A regulators to provide quantitative evidence for determining the number and members of possible clusters in a dataset [10]. Principal component analysis (PCA) was used to distinguish the distances between the m6A patterns [11].

Screening DEGs between m6A patterns
The “limma” package of R software was used to analyze and screen DEGs between two m6A patterns. P < 0.05 and |log FC| > 1 were used as screening criteria, as described in the previous report [12].

Gene set enrichment analysis
Gene set enrichment analysis (GSEA) was performed using the “clusterprofiler” R package, which was used for Gene Ontology (GO) functional enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses to explore enrichment pathways and functional annotations, genes were screened by P < 0.05 and normalized enrichment scores |NES| > 1 signaling pathway. Bioinformatics (http://www.bioinformatics.com.cn/) was used to integrate and visualize the above results.

Construction of protein‒protein interaction networks
We used the String database to construct a protein‒protein interaction (PPI) network [13]. We selected genes with a combined score ≥ 0.4 and removed the noninteracting proteins. To assess potential PPI relationships, we sent the results to Cytoscape 3.10 software (https://cytoscape.org/) and constructed the PPI network. The degree of expression of each protein was subsequently calculated, and the genes in the top 5 genes according to the degree score were screened as hub genes.

LASSO regression and SVM-RFE analysis
The least absolute shrinkage and selection operation (LASSO) and support vector machine recursive feature elimination (SVM-RFE) algorithms were used to analyze important biomarkers in CRA and CRC patients. LASSO is a regression analysis algorithm that filters variables to prevent overfitting. Based on the expression of DE-m6A regulators, LASSO regression analysis was performed using the “glmnet” package, and the calculation was performed by 10-fold cross-validation. The SVM-RFE algorithm was executed by the “e1071” software package and was used to screen the best variables.

RNA extraction and quantitative real-time polymerase chain reaction (qRT-PCR)
Sixty tissues were collected from the groups (20 in each group) and independently diagnosed by two pathologists from July 2020 to January 2022 in the Department of Gastroenterology of the Affiliated Hospitals of Anhui Medical University. This study was approved by the Ethics Committee of Anhui Medical University (Approval Number: 20190292).
Total RNA was extracted from the biological samples with TRIzol reagent (TransGen Biotech, Beijing, China). The cDNA was obtained by reverse transcription according to the kit instructions (Sparkjade, Shandong, China). The amplification reaction mixture was prepared using a real-time PCR kit (Sparkjade, Shandong, China), and the CT value of the target gene was set to be detected using a LightCycler 96 (Roche Diagnostics, USA). The sequences of primers used are shown in Supplementary Tables 1, and the primers used were synthesized by Shanghai Sangon Bioengineering Co., Ltd.

Dot blotting
The 50 ng/µL RNA was denatured at 95 °C for 3 min and then cooled on ice. The RNA was spotted sequentially on labelled nylon membranes and irradiated with a UV lamp at 254 nm for 30 min. The membranes were subsequently washed with TBST, blocked with 5% nonfat milk, incubated overnight with anti-m6A (1:1000, Abclonal, Wuhan, China), incubated with a secondary antibody (1:10000, Abclonal, Wuhan, China) and subjected to chemiluminescence (Tanon, Shanghai, China). Afterwards, the membranes were stained with methylene blue and imaged under white light.

Statistical analyses
R 4.1.3 software was used for statistical analysis and graphing. Normally distributed continuous data are expressed as the mean ± standard deviation, and comparisons between groups were analyzed using the t test. Count data are expressed as rates, and the chi-square test was used for comparisons between groups. Skewed distribution data were compared with Wilcoxon signed-rank tests. P < 0.05 was considered to indicate statistical significance. The data were plotted using GraphPad Prism 8.0.1 software.

Results

Identification and correlation between DE-m6A regulators
Based on the GSE41657 dataset, GEO2R, a built-in online tool from the GEO (https://www.ncbi.nlm.nih.gov/), was used to identify DE-m6A regulators between CRA patients and controls and between CRC patients and controls (P < 0.05, |log FC| > 0.35). Fifteen DE-m6A regulators, namely, IGF2BP1, IGF2BP3, YTHDC2, IGF2BP2, YTHDF2, YTHDF1, ZC3H13, HNRNPA2B1, WTAP, ALKBH5, METTL3, YTHDF3, RBMX, HNRNPC and LRPPRC, were found to be differentially expressed between CRA patients and controls (Supplementary Table 2). Thirteen DE-m6A regulators, including IGF2BP3, METTL14, IGF2BP2, YTHDF2, LRPPRC, HNRNPC, RBMX, YTHDF1, ZC3H13, IGF2BP1, ALKBH5, KIAA1429 and YTHDF3, were found to be differentially expressed between CRC patients and controls (Supplementary Table 3). To verify the stability of DE-m6A regulators in CRA and CRC, we analyzed the DE-m6A regulators in CRA and CRC in the GSE117606 and GSE100179 datasets. GSE41657 included most of the DE-m6A regulators (Supplementary Table 4), and we used the GSE41657 dataset for subsequent analysis. We used the GSE41657 dataset to evaluate the expression of twenty m6A regulators in CRA patients, CRC patients, and controls (Fig. 1a). Compared with controls, HNRNPA2B1, METTL3, and WTAP were all low-expressed in CRA, YTHDC2 was only overexpressed in CRA, KIAA1429 was high-expressed in CRC, METTL14 was low-expressed in CRC. In addition, two regulators were downregulated in both CRA and CRC patients, and nine regulators were upregulated in both CRA and CRC patients.
To further understand the interaction between m6A methylation regulators and possible intrinsic characteristics, two correlation heatmaps with signs were performed on the DE-m6A regulators of CRA and CRC (Fig. 1b-c) using the OmicStudio tools at https://www.omicstudio.cn. The results showed that, in CRA patients, IGF2BP3 was negatively correlated with WTAP and ALKBH5. IGF2BP2 was negatively correlated with WTAP and positively correlated with YTHDF2. YTHDF1 was significantly positively correlated with HNRNPC, ZC3H13, etc. HNRNPA2B1 and WTAP were positively correlated with ALKBH5, METTL3, RBMX, LRPPRC, etc. The interaction of m6A regulators in CRC was not obvious (Fig. 1c). HNRNPC was positively correlated with RBMX, ZC3H13 was positively correlated with YTHDF1, and KIAA1429 was positively correlated with YTHDF3. To explore the relationship between m6A regulators and the prognosis of CRC, survival analysis revealed that YTHDF3 and LRPPRC were significantly correlated with overall survival in CRC patients (Fig. 1d, and Supplementary Fig. 1a-m) based on the survival data of CRC patients in The Cancer Genome Atlas (TCGA).

Consensus clustering analysis of two m6A patterns
Based on the identification of fifteen DE-m6A regulators in CRA and thirteen DE-m6A regulators in CRC, we used the “ConsensusClusterPlus” package to cluster CRA and CRC, and when k = 2, CRA and CRC were clustered into two subgroups (Fig. 2a and c) with a more stable grouping effect. Twenty-five CRA patients were clustered into Cluster 1, and twenty-six CRA patients were clustered into Cluster 2. For CRC, Cluster 1 had eleven samples, and Cluster 2 had fourteen samples. PCA clearly revealed that fifteen DE-m6A regulators could distinguish CRA patients into two m6A groups (Fig. 2b). The effect of grouping CRC patients based on DE-m6A regulators was not good (Fig. 2d). A heatmap was constructed to show the expression of the fifteen m6A regulators and the clinical information of the two subgroups of CRA patients (Fig. 2e). Pathological features and other clinical information were not significantly different between cluster 1 and cluster 2 (P > 0.05) in patients with CRA. The expression and clinical information of the thirteen DE-m6A regulators in CRC are shown in Fig. 2f. There was a significant difference in sex between the two CRC subgroups (P < 0.05). N stage refers to the stage with or without lymph node metastasis. There were significant differences in the degree and number of lymph node metastases between the two groups (P < 0.05). These results suggested that the DE-m6A regulators of CRC were associated with the prognosis of CRC patients.

GSEA of DEGs between the two m6A patterns
To explore the potential function of these m6A-related DEGs in CRA and CRC, the “limma” package was used to analyze DEGs in the two subgroups, and GSEA was used to enrich DEGs. A total of 6,917 GO entries and 184 KEGG pathways were annotated in the CRA cohort. GO functional enrichment analysis revealed that the m6A-related DEGs were involved mainly in ribosomal structural components, signal receptor regulation, and mitochondrial structure and function (Fig. 3a-c). The first ten terms are significantly enriched in the Supplementary Table 5. KEGG pathway analysis revealed that these genes were significantly associated with the TCA cycle, pyrimidine metabolism, ubiquitin-mediated proteolysis, etc. (Fig. 3d), and the first ten terms are significantly enriched in the Supplementary Table 6. In CRC, the six genes associated with m6A include ZNF682 and ZNF439 from the ZNF family; KLK10, KLK8 and KLK7 from the KLK family; and HS6ST1 (Supplementary Table 7).

Construction of the PPI network and selection of hub genes
A PPI network of the m6A-associated DEGs in CRA patients was constructed using the String tool (https://cn.string-db.org/), and five hub genes were selected from the PPI network using the “CytoHubba” plugin of Cytoscape software (Fig. 4a): TP53, score = 4; EEF1A1, score = 20; TBP, score = 15; MARS1, score = 15; and CCT7, score = 15. The expression of the hub genes in tumor tissues was analyzed using the GSE41657 mRNA expression profile, and TBP, MARS1, and CCT7 were found to be highly expressed in CRC (P < 0.05) (Fig. 4b-f). Survival analysis of CRC patients in the TCGA cohort revealed that TBP and TP53 were associated with overall survival (P < 0.05) (Supplementary Fig. 2a-e). We also found that the hub genes had high value as diagnostic markers for CRA and CRC (Supplementary Fig. 2f-g).

Construction of SVM-RF and LASSO regression models for screening marker genes
We constructed LASSO regression and SVM-RFE models to screen important m6A regulators among the 15 DE-m6A regulators in CRA (Fig. 5a-b). Based on the minimum mean square error, 10-fold cross-validation was used to select nine m6A regulators (YTHDF1, YTHDF2, YTHDF3, YTHDC2, IGF2BP1, LRPPRC, METTL3, HNRNPC, and HNRNPA2B1) as candidate marker genes by LASSO. In addition, eight key DE-m6A regulators (METTL3, HNRNPA2B1, WTAP, IGF2BP1, YTHDF3, YTHDF2, RBMX, and YTHDF1) were identified as candidate marker genes via SVM-RFE (Fig. 5c). The candidate marker genes common to both models were selected using a Venn diagram (METTL3, HNRNPA2B1, IGF2BP1, YTHDF3, YTHDF2, and YTHDF1) (Fig. 5d). Then, we constructed a nomogram model based on six candidate DE-m6A regulators to predict the prevalence of CRA patients using the “rms” package (Fig. 5e-f). Moreover, we screened a key DE-m6A regulator, METTL14 (Supplementary Fig. 3a-d). Finally, we detected the predictive power of seven candidate genes in GSE117606 and GSE100179 (Fig. 5g-h; Supplementary Fig. 3e-f), in addition to YTHDF1 in CRA (AUC = 0.723 and AUC = 0.805) and CRC (AUC = 0.849 and AUC = 0.873); moreover, HNRNPA2B1 had high diagnostic efficacy only for CRC (AUC = 0.769 and AUC = 0.737).

qRT‒PCR verification of DE-m6A regulators
The expression levels of significant m6A regulators were verified. Samples from twenty CRA patients, twenty CRC patients, and twenty normal controls were randomly selected from our perspective observation cohort and subjected to qRT‒PCR, and there was no significant difference in age or sex among the three groups. Compared with those in the controls, the expression of YTHDF1, YTHDF2, and YTHDF3 was significantly greater in CRC patients; that of IGF2BP1 and HNRNPA2B1 was greater in CRA patients; and that of METTL3 was elevated but not significantly elevated in adenoma patients, which was consistent with the above results (Fig. 6a-f). In addition, dot blot analysis was used to measure the m6A modification level in the samples, and the results (Fig. 6g) showed that the m6A modification level was high from CRA to CRC and was significantly greater than that in the control group. Furthermore, m6A modification dysregulation may be involved in the occurrence of CRC, and disordered expression of m6A regulatory factors may be involved in the occurrence of CRC by affecting the m6A modification level. The ROC curve of CRA showed that the AUC of IGF2BP1, HNRNPA2B1, YTHDF2 and METTL3 were 0.890, 0.874, 0.871 and 0.792, respectively. The ROC curve of CRC showed that the AUC of YTHDF1, YTHDF2, HNRNPA2B1 and YTHDF3 were 0.958, 0.912, 0.883 and 0.825, respectively. ROC results indicated that the above gene expression may be valuable biomarkers in CRA and CRC (Supplementary Fig. 3a-b).

Discussion
CRC is a global public health problem with increasing morbidity and mortality worldwide, with 1.4 million newly diagnosed cases and 600,000 deaths each year [1]. Adenomatous polyps, which are intestinal polyps, easily develop into adenocarcinomas within an average period of 8–10 years. However, rapid neoplastic transformation of adenomatous polyps to adenocarcinoma was also observed within one year [14]. Due to the insidious symptoms of intestinal polyps and the low diagnostic rate, in addition to the high recurrence rate, early screening and early treatment of gastrointestinal polyps are particularly important for the prevention of CRC.
Compared with those in normal controls, we found fifteen DE-m6A regulators in CRA patients and thirteen DE-m6A regulators in CRC patients. On the basis of the findings of a large number of previous studies, we found that these DE-m6A regulators are clearly related to the cell cycle and tumor cell migration and progression. However, few reports have evaluated the relationship of these genes with CRA and CRC [15, 16]. According to these DE-m6A regulators, both CRA and CRC patients were divided into two subgroups according to the two m6A modes. Next, we analyzed DEGs between the two m6A molecular patterns. We then enriched the functional pathways in which these DEGs participated. Furthermore, we found that the PPI network revealed five hub genes in CRA, and the expression of EEF1A1, TBP, MARS1 and CCT7 was greater in CRC than in CRA. TBP proteins in mammals are divided into three families, TRF1, TBPL1, and TBPL2, according to the conserved nature of their C-terminal core domain, and the complex role of TBP-like factor complexity in transcriptional regulation is worth exploring. Ultimately, based on the DE-m6A regulators, we constructed a regression model and an algorithm to identify and verify the diagnostic markers of CRA and CRC. Using qRT‒PCR, we also confirmed the expression of the m6A genes METTL3, HNRNPA2B1, IGF2BP1, YTHDF3, YTHDF2, and YTHDF1 in CRA and CRC cases, which was consistent with our bioinformatics results and previous published results of those in CRC cases [17, 18].
Epigenetic modifications such as DNA methylation, histone acetylation, and RNA modification are involved in the occurrence and development of tumors and have been recognized as new therapeutic and prognostic targets. To date, a growing number of posttranscriptional modifications of RNA have been recognized [19]. m6A modifications are common in mRNAs, lncRNAs, and miRNAs [20] and are the most common reversible mRNA modifications in eukaryotes. These modifications are involved in many biological processes. Its modification process is mediated by the methylase “Writer”, which can be reversed by the demethylase “Eraser”, and the above process involves a series of “readers” that recognize and bind the target RNA. METTL3 and METTL14 act as “writers” of m6A, forming heterodimeric complexes through the support of WTAP and fulfilling the role of methylases [21]. FTO and ALKBH5 act as m6A “erasers” that not only reverse methylations but also modulate the stemness of cancer stem cells [22, 23]. m6A regulators play indispensable roles in numerous biological processes. However, the role of m6A regulators in CRA and CRC remains unclear.
These m6A-related DEGs were involved in ribosomal structural components, signal receptor regulation, and mitochondrial structure and function in CRA, and the pathways were related to the TCA cycle, pyrimidine metabolism, and ubiquitin-mediated proteolysis. Others included chromosomal organization, cell cycle regulation, intracellular protein transport, cellular macromolecule catabolic processes, organelle envelopes, envelopes, hydrolase activity, action on anhydrides, triphosphatase activity, and ATPase activity. Previous studies have shown that organelle ribosomes are involved in transcription and translation [24]. During this process, changes in the quantity and quality of ribosomes could promote tumorigenesis [25]. Proliferating cancer cells are known to exhibit abnormal glucose metabolism, which is characterized by high glucose absorption and high glucose solubility [26]. The metabolic flux increases through glycolysis and metabolic inhibition of the TCA cycle, which provides a biosynthetic precursor for the rapid synthesis of polymers and maintains the cellular redox balance to increase the survival rate [27]. O-GlcNAcylation may coordinate important signals of glycolysis and the TCA cycle to promote tumorigenesis [28]. The progression of CRC is thought to be dependent on glycolysis and dysfunctional mitochondria [20]. Recent studies have shown that most people with CRC have mitochondrial DNA mutations [2]. However, the effect of these mutations on the respiration capacity of patients with mitochondrial DNA mutations is uncertain. These findings suggested that mitochondria may be involved in the progression of CRC.
The polyglutamine region of the TBP protein contains structural polymorphisms that could be associated with a variety of hereditary disorders when folded and dilated, and its amplification coexists with spinocerebellar ataxia 17 (SCA17), a Huntington’s disease-like phenotype, and Parkinson’s disease [29–31]. TBP-related factor 1 (TRF1) is expressed in Drosophila and Anopheles mosquitoes but not in humans [32]. TBPL1 may be involved in the transcriptional program responsible for the construction of the pericentromeric region and may contribute to the formation of new karyotypes, particularly by facilitating centromeric fusion [33]. TBPL2 plays a very important role in gene expression and regulation in growing oocytes [34]. MARS1 is a tRNA synthase gene, and mutations in MARS1 may lead to ISR activation, which drives neurodegenerative diseases [35]. The chaperone protein (CCT) of TCP-1 promotes the folding of intracellular proteins (mainly cytoskeletal proteins, such as tubulin and actin) in the cytoplasm and participates in the development of cancer, and the η subunit of CCT is encoded by CCT7 [36]. CCT7 overexpression is associated with the prognosis of HCC patients [37]. Although few studies have explored the role of these hub genes in CRA and CRC, we found that these hub genes have high value as markers for the diagnosis of CRA and CRC, providing clues for the study of diagnostic markers of CRA and CRC.
Still, there were some limitations to the present study. The data set lacked complete clinical characteristics of patients with CRA and CRC, as well as adenoma location, number of adenomatous polyps, and hemafecia conditions, which were closely related to CRA and CRC. Although our study effectively divided CRA and CRC into two subgroups, namely, the expression patterns of two m6A regulatory factors, However, the association between the clinical features of the two subgroups and the expression pattern of m6A regulatory factors could not be provided. In addition, there is a continuous process from colorectal adenoma to carcinogenesis, and the gene sequencing data of CRA and CRC come from two different populations respectively, which cannot better observe the continuous changes of markers in CRA and CRC, and the small sample size may lead to unstable differences between groups, so prospective investigation and large sample detection need to be solved in future studies.

Conclusion
Taken together, our findings suggest that m6A modulators are expressed specifically in CRA and CRC, suggesting a potential mechanistic role for m6A modulators in these molecular pathways, which provides new insights into pathological analysis and diagnostic biomarker exploration in the future.

Supplementary Information