The YTH Domain Family of N6-Methyladenosine “Readers” in the Diagnosis and Prognosis of Colonic Adenocarcinoma

To profile the landscape of methylation N6 adenosine (m6A) RNA regulators in colonic adenocarcinoma (COAD) and to explore potential diagnostic and prognostic biomarkers, we assessed the differential expression patterns of m6A RNA methylation regulators between 418 COAD patients and 41 controls based on profiling from The Cancer Genome Atlas (TCGA) database. We plotted the receiver operating characteristic (ROC) curves and calculated the area under the curve (AUC) values to estimate the discrimination ability. The relationship between the expression of m6A RNA methylation regulators and clinicopathological characteristics was explored. Kaplan-Meier plotter, log-rank test, and Cox regression were used and a nomogram was created to explore the prognostic significance of m6A-related genes in overall survival at the mRNA level. Pathway analysis was performed by gene set enrichment analysis (GSEA) using TCGA dataset, and a coexpression network was built based on the STRING database. We observed that YTHDF1, METTL3, and KIAA1429 were significantly upregulated, while YTHDF3, YTHDC2, METTL14, and ALKBH5 were significantly downregulated in COAD samples compared to normal samples. YTHDF1 had the highest diagnostic value. Low expression of YTHDF3 predicted a poor survival rate in COAD patients. YTHDC2 was related to sex and showed a downward trend as clinical stage increased. Our results indicate that the YT521-B homology (YTH) domain family (“readers”), especially YTHDF1, YTHDF3, and YTHDC2, might play a significant role in the detection, progression, and prognosis of COAD, indicating that they are promising cancer biomarkers.


Introduction
Colon cancer is the third most common cancer and the second leading cause of cancer-related death worldwide [1]. Colonic adenocarcinoma (COAD) is the typical type that accounts for 98% of new cases [2]. Despite advances in diagnosis and treatment [3], the prognosis of patients remains poor. Chemotherapy has shown significant value, but surgery is the only method of curative treatment. Thus, it is crucial to discover novel diagnostic and prognostic biomarkers for COAD.
Epigenetic modification has gained increasing attention in research on carcinogenesis and progression. Methylation is a common epigenetic trait and a simple biochemical process-it is the transfer of four atoms, one carbon atom and three hydrogen atoms (CH3), from one substance to another [4]. The N6-methyladenosine (m 6 A) modification of mRNA/lncRNA is the most common RNA methylation, and it can regulate cell fate determination, self-renewal, and cancer development, indicating it as a promising therapeutic target for cancer [5]. Three major types of enzymes participate in m 6 A methylation: "writers," "erasers," and "readers." "Writers" mainly include METTL3, METTL14, KIAA1429, and WTAP. METTL3 and METTL14 form a structure with other enzymes (e.g., WTAP) to methylate the sixth N of adenosine during transcription from DNA to RNA [6]. The m 6 A-modified bases can be demethylated under the action of "erasers" enzymes, such as FTO and ALKBH5, making RNA methylation a reversible reaction [7]. These methylated   2 BioMed Research International RNA base sites finally require specific enzymes ("readers") to recognize. YTHDF1, YTHDF2, and YTHDF3 are primary members of the "readers" family, devoted to recognizing bases that undergo m 6 A methylation, participating in downstream translation, mRNA degradation, and accelerating the rate at which mRNA leaves the nucleus [8]. Aberrant expression of m 6 A writers, erasers, or readers can lead to the deregulation of m 6 A modification, thus triggering the abnormal translation of specific mRNAs and promoting or inhibiting carcinogenesis and cancer progression [9]. Accumulative evidence has supported the correlation between aberrant cellular m 6 A and human cancers [10]. For example, the overexpression of YTHDF1 is related to poor prognosis in patients with liver cancer [11]; downregulation of ALKBH5 in MDA-MB-231 human breast cancer cells can lead to a decrease in the number of breast cancer stem cells, resulting in a significant reduction in tumorigenic capacity [12].
Although traditional prognostic factors for COAD, such as age, tumor stage, surgical margins, and tumor grade, have produced significant improvements in predicting patient clinical outcomes [13], their limitations in distinguishing the prognostic value of molecular heterogeneity should not be ignored. With deepening research on RNA methylation regulators, the roles of m 6 A in diagnosis and treatment have gradually become valued. This study is aimed at profiling the landscape of m 6 A RNA methylation regulators in COAD and at exploring their potential value as diagnostic and prognostic biomarkers.

Materials and Methods
2.1. Data Processing. We downloaded the expression data and clinicopathological data of 418 COAD patients and 41 con-trols after the deletion of missing records from TCGA database (https://tcga-data.nci.nih.gov/tcga/) [14]. Age at diagnosis, sex, follow-up days, and clinical data were retrospectively extracted. All patients were staged using the 2009 TNM classification system recommended by the American Joint Committee on Cancer. We used the edgeR package to normalize the expression data.

2.2.
Portraying the Landscape of m 6 A RNA Methylation Regulator Expression. The expression levels of six m 6 A RNA methylation regulators between tumor and normal tissues were compared using the t-test. The expression pattern of m 6 A RNA methylation regulators in COAD samples was denoted by the Pearson correlation matrix, and a correlation coefficient (r) was calculated. We plotted the ROC curve and calculated the area under the curve (AUC) to estimate the discrimination ability. The m 6 A RNA methylation regulators and clinicopathological features were analyzed using the t-test or one-way ANOVA.

Survival Analysis.
A total of 418 patients were followed, with a maximum follow-up period of 4502 days. The median survival time was 670.5 days. All patients had no chemotherapy history and underwent the same type of radical surgery and postoperative adjuvant chemotherapy. The overall survival (OS) was calculated from the date of diagnosis to the date of death. We categorized the expression of m 6 Arelated mRNA into two groups using the lower quartile as the cutoff point. Kaplan-Meier analysis was used to portray the survival curves of m 6 A RNA methylation regulators. The log-rank test was used to compare the survival distributions between groups. We performed univariate and multivariate analyses to determine the independent prognostic factors using the Cox proportional hazard model. A nomogram was generated based on the significant prognostic factors in the Cox regression model to predict the 1-year, 3-year, and 5-year survival of COAD patients.

Biological Function Analysis.
We further used gene data in TCGA database to search for pathways related to YTHDF1, YTHDF3, and YTHDC2 by gene set enrichment analysis (GSEA) using GSEA v2.2.2 software (http://www .broadinstitute.org/gsea). GSEA revealed significant differ-ences in the enrichment of genes in the MSigDB Collection (c2.cp.kegg. v7.2. symbols) [15]. The high-and lowexpression phenotype groups were divided according to the lower quartile expression level of candidate genes. Gene set permutations were performed 1000 times for each analysis. The expression of selected genes was used as a phenotype label. We selected the most significantly enriched signal pathway using the criterion of NOM P value < 0.05, FDR q value < 0.25, and high normalized enrichment score (NES). The Retrieval of Interacting Genes (STRING v10) (http:// string-db.org/) tool was used to analyze the interactive  5 BioMed Research International relationships to construct a protein-protein interaction (PPI) network. Experimentally validated interactions with a combined score > 0:4 were regarded as significant. We also selected the YTHDF1, YTHDF3, and YTHDC2 modules and constructed a PPI network of these selected m 6 A RNA regulators.

Statistical Analysis.
A comparison of normalized data between two groups was conducted using t-test or Mann-Whitney U test by GraphPad Prism 7 (GraphPad Software, La Jolla, CA). Comparisons among multiple groups were performed using one-way ANOVA. Other analyses were visualized using R software (v.3.4.3). The significant level was set at 0.05.

m 6 A RNA Methylation Regulators and Clinicopathological
Features. Compared to younger patients, patients older than 65 years had a lower expression level of YTHDC1 (P < 0:0001) and a higher expression level of ALKBH5 (P = 0:011) (Figure 2(a)). YTHDC2 (P = 0:028) was highly expressed in female patients (Figure 2(b)). Clinical stage was associated with YTHDF1 ( Figure 3 Figure 5) using the stepwise forward method. YTHDF3, together with age, N stage, and M stage, was an independent prognostic factor for COAD (Table 1).

Nomogram to
Predict the Survival of COAD. Data on YTHDF3 expression, age, and clinical stages from TCGA dataset were used to establish a prognostic nomogram predicting overall survival based on the stepwise Cox regres-sion model. Total score was obtained by adding up the individual contributions of the corresponding covariates on the points scale; the total score was used to predict 1-year, 3-year, and 5-year related survival probabilities ( Figure 6).

GSEA Identified Related Signaling Pathways.
Considering the critical associations of YTHDF1 with diagnosis, YTHDF3 with prognosis, and YTHDC2 with clinicopathological features observed by the above analysis, we further used the GSEA approach to explore the potential biological processes related to these genes. YTHDF1 mainly participates in basal transcription factors and spliceosomes (Supplementary Figure S1A Figure S1E, F, and G).
3.6. PPI Network Construction. Based on the STRING database, a coexpression network was constructed to explore interactions between m 6 A RNA methylation regulators (Supplementary Figure S2). We found that the writers (WTAP, KIAA1429, METTL3, and METTL14) closely interacted with each other (Supplementary Figure S2A). Further, we explored three module networks of selected m 6 A regulators (YTHDF1, YTHDF3, and YTHDC2) and broader genes. YTHDF1 was firmly related to YTHDF3 and shared many interacting proteins. Moreover, the writers (WTAP, KIAA1429, METTL3, and METTL14) also played essential roles in this network (Supplementary Figure S2B). In addition, many other genes were involved in forming the network, showing ample room for exploration.

Discussion
COAD is one of the most popular malignant gastrointestinal tract tumors and the second leading cause of cancer-related death in adults in Western countries [16]. During the past decades, many efforts have been made to search for biomarkers of COAD. Studies have shown the diagnostic value of DCTN1, DCTN2, and DCTN4 [17] and the prognostic role of LAYN [18], KCNQ1OT1 [19], and PYK2 [20] in COAD, but with limited sample size or inconsistent results. In this study, we analyzed the expression patterns of m 6 A RNA methylation regulators in COAD and their relationship with clinical features by using cohort data from TCGA. Most m 6 A-related proteins were dysregulated in COAD samples, and some were associated with clinical features. Among them, YTHDF1 had the highest diagnostic value for COAD. Survival analysis confirmed that high   10 BioMed Research International expression levels of YTHDF2 and YTHDF3 predicted a favorable prognosis. Moreover, YTHDF3 expression levels were independent prognostic factors of 5-year overall survival in COAD patients. YTHDC2 was highly expressed in female patients and showed a significant relationship with tumor progression. YTHDC1 was typically upregulated in younger patients. These results underscore the significant roles of the YTH domain family of readers, especially YTHDF1, YTHDF3, and YTHDC2, in the detection, progression, and prognosis of COAD.
In this study, we explored the roles of the YTH domain family of readers in COAD. m 6 A exerts many of its functions through reader proteins in the cytoplasm and nucleus that selectively bind to m 6 A-containing transcripts [21]. The "readers" are the bridge between "writers" and "erasers", and the ability of "readers" to block the accessibility of RNA demethylase may be critical to determining the m 6 A status of target transcripts [10], which may explain the critical roles of the YTH family in COAD. YTH domain family proteins include YTHDF1, YTHDF2, and YTHDF3 in the cytoplasm and YTHDC1 and YTHDC2 in the nucleus [22,23]. YTHDF1 and YTHDF3 work in concert to affect the translation of m 6 Acontaining mRNAs, YTHDF2 expedites mRNA decay, YTHDC1 affects the nuclear processing of its targets, and YTHDC2 selectively binds m 6 A at its consensus motif and enhances the translation efficiency of its targets and decreases their mRNA abundance [24][25][26]. Consistent with our results, Bai et al. reported that YTHDF1 could regulate tumorigenicity and cancer stem cell-like activity in colorectal carcinoma [27], and Tanabe et al. reported that YTHDC2 was upregulated in colon cancer and had a positive correlation with tumor stages [28]. Although YTHDF2, YTHDF3, and YTHDC1 showed essential roles in multiple cancer types [29][30][31][32], to date, no research on these regulators in COAD has been conducted. Therefore, the role of YTHDF2, YTHDF3, and YTHDC1 in COAD needs to be further studied. Our study showed a negative relationship between METTL14 and clinical stage, which implied its critical role in COAD progression. Ma et al. reported that METTL14 inhibits the metastatic potential of hepatocellular carcinoma by modulating N6methyladenosine-dependent primary microRNA processing [33]. Further research is needed to explore the function of METTL14.
Moreover, our study showed a high correlation between the expression of YTHDF3 and KIAA1429 in COAD. We also observed a moderate correlation between METTL14 and YTHDF2, YTHDC1, and YTHDC2. The PPI network showed the interactions among "writers", "readers", and "erasers". The cross-talk among different categories of m 6 A regulators may influence cancer growth and progression. Chen et al. found that METTL3 promoted liver cancer progression through YTHDF2-dependent posttranscriptional silencing of SOCS2 [34]. Panneerdoss et al. reported that METTL14 and ALKBH5 inhibited the expression of YTHDF3, reversely block RNA demethylase activity, and altered the m 6 A status of target transcripts in cancer cells [10]. These studies have emphasized interactions between writers, readers, and erasers. Moreover, the interactions of members in the same category are also frequent. For example, YTHDF3 initiates mRNA translation and methylated mRNA decay through cooperation with YTHDF1 and YTHDF2, affecting the functions of the other two YTHDFs [25]. Therefore, the interactions of these regulators may play essential roles in COAD and should be further explored.
GSEA revealed that YTHDF1, YTHDF3, and YTHDC2 were enriched in several meaningful pathways. YTHDF3 and YTHDC2 share the same pathways of ubiquitinmediated proteolysis, which further illustrates their functional biological connection. Taniue et al. found that lncRNA UPAT and UHRF1 may be promising molecular 11 BioMed Research International targets for the therapy of colon cancer by regulating protein ubiquitination and degradation and thereby play a critical role in the survival and tumorigenicity of cancer cells [35]. We suggest that the ubiquitin-mediated proteolysis affected by YTHDF3 and YTHDC2 may play a crucial role in tumor progression. The spliceosome can be induced for the activation of splicing and mRNA production in the carcinogenic process. Takayama et al. suggested that targeting spliceosome proteins by RNA-binding protein to promote AR splicing and expression could be a therapeutic possibility for hormone-refractory prostate cancer [36]. We hypothesize that YTHDF1 could also be a potential therapeutic target by modulating spliceosome and other RNA activities. TGFβ promotes tumor growth and metastasis by inducing angiogenic factors and facilitating EMT [37]. TGF-β signalingassociated genes are particularly sensitive to changes in m 6 A levels [10]; thus, the malfunction of the TGF-β signaling pathway in response to YTHDF3 dysregulation may be the main factor affecting the prognosis of COAD.

Conclusions
Most m 6 A-related proteins are dysregulated in COAD samples compared to normal samples, and some are associated with clinical features. Our results indicated a significant role of the YTH domain family ("readers") in the diagnosis, progression, and prognosis of COAD. The exact mechanism of the YTH domain family of N6-methyladenosine readers is worth further study.

Data Availability
All data generated or analyzed during this study are included in this published article.

Ethical Approval
This study was approved by the ethics committee of Nanjing Medical University.

Consent
All data were downloaded from the public database and followed the data access policy. All subjects were anonymous. So, there was no requirement for informed consent.

Disclosure
The funding agencies had no role in the study design, data collection, analysis, decision to publish, or preparation of the manuscript.