Transcriptome-wide Landscape of RNA Modication N6-Methyladenosine in Chromophobe Renal Cell Carcinoma

As the most abundant internal mRNA modication, N 6 -methyladenosine (m 6 A) is associated with various cancers. However, RNA modication m 6 A has not been studied in chromophobe renal cell carcinoma (chRCC). The present study aimed to comprehensively analyze the global m 6 A modication pattern in chRCC. Three subjects with chRCC were enrolled in our study. Transcriptome-wide m 6 A methylone and transcriptome analysis in chRCC and tumor-adjacent normal tissues were detected via m 6 A-SEAL-seq and RNA-sEq. m 6 A-modied mRNAs were further validated by m 6 A-immunoprecipition followed by quantitative real-time PCR (m 6 A-IP-qPCR). The least absolute shrinkage and selection operator (LASSO) Cox regression and multivariate Cox proportional hazards regression analysis were used to determine the candidate gene. two hypomethylated and FGFR1, which respectively a oncogene in chRCC. Three m 6 A-dependent signatures were identied using Cox regression screen and Based on the signicant prognostic signatures, we build a m 6 A-dependent prognostic model, with the Concordance index (C-index) = 0.96. The Kaplan-Meier survival curve and log-rank between the high-risk and low-risk group showed signicant difference.

end processing [15,16], alternative splicing [17,18], and translation e ciency at the post-transcriptional level [14,19,20]. m 6 A is a dynamically reversible RNA modi cation, which is regulated by "writers" (methyltransferases), "erasers" (demethylases), and "readers" (binding proteins). The majority of m 6 A modi cations are installed through the methyltransferase complex containing key catalytic subunits METTL3-METTL14 heterodimer [21,22] and other subunits like WTAP [23][24][25], and removed by demethylases like FTO and ALKBH5 [26,27]. m 6 A modi cation is recognized by m 6 A binding proteins, such as YTH domain family proteins, for regulation of RNA processing and metabolism [28][29][30]. The discoveries of these different m 6 A regulators contribute to a better understanding of the physiological functions of m 6 A.
Numerous studies have demonstrated a close relationship between m 6 A modi cation and tumor progression [31][32][33][34]. For instance, in bladder cancer, METTL3 installs m 6 A in pri-miR221/222 and accelerates the maturation of miRNAs, leading to the proliferation of cancer cell [35]. In pancreatic cancer, ALKBH5 inhibits cancer cell growth and progression through increasing PER1 mRNA levels by demethylating m 6 A modi cation in PER1 and subsequently escaping from YTHDF2-mediated mRNA decay [34]. These ndings suggest that m 6 A modi cation plays vital roles in carcinogenesis through the regulation of RNA processing and metabolism and provide new molecular mechanisms of cancer progression. Similarly, m 6 A regulators also have signi cant impacts on ccRCC [36]. METTL14 is downregulated in ccRCC tissue and patients with lower METTL14 expression tend to have worse prognoses [37], and the alteration of m 6 A regulators is associate with worse clinical characteristics [38].
In RPCC, a prognostic risk signature model with three m 6 A regulatory genes, IGF2BP3, KIAA1429 and HNRNPC, could predict survival outcomes accurately [39]. With the development of high-throughput sequencing, transcriptome-wide pro ling of m 6 A distribution in multiple different human carcinomas becomes available, which gives a way to interpret the molecular mechanisms between m 6 A modi cation and RCC. In 2020, transcriptome-wide m 6 A mapping in ccRCC were reported, and the unique m 6 A-related genes in ccRCC are associated with cancer-related pathways, providing a possible mechanism of m 6 Amediated gene regulation [40]. However, the transcriptome-wide distribution of m 6 A in chRCC has not been gured out yet. Here, we report the transcriptome-wide m 6 A pro ling in human chRCC by the use of antibody-free methods m 6 A-SEAL-seq [41], which is a FTO-assisted chemical labeling m 6 A sequencing method and has good sensitivity, speci city and reliability. This study will be helpful for providing a basis for more in-depth studies of the biological functions of m 6 A in pathogenesis of human chRCC.

Patients and specimens
A total of six patients with chRCC were involved in our study. chRCC tissues and corresponding tumoradjacent normal tissues were collected at the time of surgery from urology department, Peking University Third Hospital. All specimens were immediately separated into 1.5 ml RNase-free centrifuge tubes and stored at -80 • C before RNA isolation.

RNA Preparation
Total RNA was extracted from tissue specimens using TRIzol reagent (Magen) and poly(A) + RNA was isolated from total RNA using oligo(dT) 25 Dynabeads (Thermo Fisher Scienti c). RNA concentration was determined using a Nanodrop ultraviolet-visible light spectrophotometer (Thermo). m 6 A-SEAL-seq and library construction Poly(A) + RNA isolated from each sample was fragmented by a magnesium RNA fragmentation module (NEB) and subjected in FTO-assisted m 6 A oxidation step. In FTO-assisted m 6 A oxidation step, the reaction was performed in 300 μl aliquots of aqueous solution containing 300 μM of (NH 4  After ethanol precipitation, DTT-treated RNA was washed by 75% ethanol and dissolved in 200 μl biotinylation buffer that contained 100 μM of MTSEA-XX-biotin (Biotum), 100 mM HEPES (pH 7.0), 1 mM EDTA, and 20% DMF. The reaction was performed at 25 °C and 800 rpm in a ThermoMixer for 1 h. The biotinylated RNA was puri ed by phenol-chloroform extraction. 50 ng of the biotinylated RNA was saved as input, and the rest was proceeded to a nity enrichment. 20 μl Dynabeads MyOne Streptavidin C1 (Invitrogen) was washed twice by 200 μl 0.1 M NaOH to remove RNase contamination, and then washed with diethyl pyrocarbonate water to a neutral pH. The beads were resuspended in 100 μl binding solution containing 10 μl of high-salt wash buffer (100 mM Tris pH 7.5, 10 mM EDTA, 1 M NaCl, 0.05% Tween 20) and 90 μl diethyl pyrocarbonate water, and incubated with the biotinylated RNA for 1 h. The biotinylated RNA on beads was washed three times with 1 ml high salt wash buffer. 50 μl of 100 mM DTT was used to release the biotinylated RNA at 37 °C for 15 min on a ThermoMixer (800 rpm.). After collecting the supernatant, the second elution was performed with 50 μl of 100 mM DTT at 50 °C for 5 min to completely release the RNA. The twice-eluted RNA was combined and puri ed by ethanol precipitation. Library construction was performed using NEBNext Ultra II Directional RNA Library Prep Kit for Illumina according to the manufacturer's protocol. Libraries were sequenced on the Illumina HiSeq XTen platform with a paired-end model (PE150).

Analysis of m 6 A-seq data
Sequencing reads were trimmed and mapped to the reference genome (GRCh38) by using Cutadapt (v1.18) [42] and HISAT2 (v2.1.0) [43], respectively. The m 6 A-enriched regions in chRCC and normal tissues were identi ed using the MACS2 [44] peak-calling algorithm based on enrichment criteria (IP/Input) ≥ 2 and FDR < 0.05. Con dent m 6 A peaks were subjected to Hypergeometric Optimization of Motif EnRichment tools (HOMER) [45] for Motif Discovery. Genes with differentially methylated m 6 A sites were identi ed by MeTDiff [46] based on enrichment criteria fold change ≥ 2 and FDR < 0.05. Tissue analysis, Gene ontology (GO) and pathway enrichment analyses were performed by using DAVID.

Analysis of RNA-seq data
Adapter and low-quality reads were trimmed by using Cutadapt (v1.18) [42], and trimmed reads were aligned to the reference genome (GRCh38) using HISAT2 (v2.1.0) [43]. The differential expression genes between chRCC and adjacent normal tissues were screened by R package (DEseq2) [47] based on a cutoff criterion of fold change ≥ 2 and FDR < 0.05.

Risk strati cation and survival analysis
A cohort of 65 chRCC cases from The Cancer Genome Atlas (TCGA) database was used to illustrate the relationship between the differential expressed DMMGs and chRCC patients. We randomly chose 40 samples from 65 cases as a training set to predict signature model and the rest samples form a testing set to verify the model (make sure the training set and testing set both contain tumor samples and normal samples). Firstly, we used least absolute shrinkage and selection operator (LASSO) to select candidate genes (glmnet package) in training set. Secondly, we performed the multi-variates cox regression and removed genes not supported by PH hypothesis using the selected candidate genes in training set. Thirdly, we performed the second regression (survival package) to calculate the coe cients between candidate genes and 5-years survival using the remaining candidate genes in training set to build a Cox model. The concordance index (C-index) was calculated to evaluate the prognostic power.
Risk score of each sample was calculated through the sum of the product of each candidate gene fpkmuq and its coe cient in training set. The patients were then classi ed into high-risk or low-risk group using the risk score where the difference value of true positive and false positive reaches to the maximum as the cutoff value. The Kaplan-Meier survival curve (survminer package) was performed to evaluate the 5-years survival, and the sensitivity and accuracy of the cox model to predict clinical outcome were evaluated by the area under curve (AUC) of the receiver operating characteristic (ROC) curve (survival ROC package). At last we test the signature model in the testing set, ccRCC dataset (a cohort of 602 cases from TCGA database) and PRCC dataset (a cohort of 318 cases from TCGA database).

The aberrant expression of several m 6 A regulators in chRCC tissue
Hematoxylin-eosin (HE) staining indicated that chRCC tissues composed of large vegetable-like polygonal cells with eosinophilic cytoplasm, irregular nuclei, perinuclear clear halo, and prominent cell membrane (Fig. 1b). In order to determine whether m 6 A modi cation functions in chRCC, we rst analyzed the expression levels of 8 m 6 A regulators, including 3 key writer subunits, 3 readers, and 2 erasers (m 6 A writer subunits: METTLE14, METTL3, and WTAP; m 6 A readers: YTHDF1, YTHDF2, and YTHDF3; m 6 A erasers: ALKBH5 and FTO) in six patients. The qPCR results showed that the expression levels of WTAP, YTHDF2, FTO, and ALKBH5 were downregulated markedly in chRCC tissues compared with corresponding tumor-adjacent normal tissues (termed normal tissues) (Fig. 1c). The aberrant expression of these m 6 A regulators in chRCC indicated that m 6 A modi cation might play a crucial role in the progression of chRCC.

Overview of m 6 A methylation feature in normal and chRCC tissues
To investigate whether m 6 A methylation landscape changes between the normal and chRCC tissues, we performed m 6 A-SEAL-seq [41] using chRCC tissues and normal tissues from three subjects. Approximately 87.1-13.2 million reads were generated from each library and 82.8-12.8 million reads were mapped to GRCh38 genome (Additional le 1: Dataset S1). m 6 A peaks were called in each sample using the published m 6 A peak caller MACS2 algorithm [44] (fold enrichment (IP/input) ³ 2 and false discovery rate (FDR) £ 0.05). The m 6 A peaks identi ed in all three replicates were classi ed as "con dent m 6 A peaks". We identi ed 15,024 con dent m 6 A peaks corresponding to 11,396 transcripts/genes in normal tissues, and 12,841 con dent m 6 A peaks corresponding to 10,102 transcripts/genes in chRCC tissue (Fig.   2a). To evaluate the reliability and performance of m 6 A-SEAL-seq, we compared our con dent m 6 A peaks identi ed in normal tissues with the published m 6 A peaks identi ed in normal kidney tissues by MeRIP-seq (GSE122744) [48]. 5686 out of 7261 (78.3%) of the m 6 A peaks in MeRIP-seq were overlapped with our identi ed m 6 A peaks in m 6 A-SEAL-seq (Additional le 2: Fig. S1a). We then plotted the distribution of distance between m 6 A peaks from m 6 A-SEAL-seq and MeRIP-seq, and found that our con dent m 6 A peaks were highly enriched around MeRIP-seq peaks (Additional le 2: Fig. S1b). We further calculated normalized read coverages from m 6 A-SEAL-seq around MeRIP-seq peaks by deepTools [49]. Their co-enrichment was further shown in Supplementary Figure 1C and 1D. These results suggest that m 6 A-SEAL-seq is accurate and reliable.
We next investigated the m 6 A distribution across transcripts in normal and chRCC tissues. The metagene pro les was used to display the distribution of m 6 A peaks across transcripts. The results showed that con dent m 6 A peaks in normal and chRCC tissues were both highly located within coding sequences (CDS) and 3′ untranslated region (3′UTR) (Fig. 2b), which was consistent with the previous observation [40]. To further locate con dent m 6 A peaks, we divided the transcripts into ve non-overlapping regions and assigned the con dent m 6 A peaks into these regions. The fraction of con dent m 6 A peaks of normal and chRCC tissues in these ve regions showed that they were dominantly enriched in 3′UTR (40.46%, 40.55%), CDS (29.01%, 26.54%) and stop codon (17.32%, 18.59%) (Fig. 2c). We clustered the con dent Further we asked which RNA molecules prefer to contain m 6 A modi cation. We assigned con dent m 6 A peaks to GRCh38 genome and found that 76.24% and 77.78% were mRNA, 18.42% and 17.16% were long non-coding RNA (lncRNA) in normal and chRCC tissues, respectively (Fig. 2d). We noticed that the con dent m 6 A peaks number in chRCC tissues were less than that in normal tissues (Fig. 2a). We subsequently assigned m 6 A peaks to chromosome, genes and the ve non-overlapping regions. We found that the number of con dent m 6 A peaks in chRCC tissues were decreased globally among each chromosome except chromosome Y, which didn't have m 6 A peaks (Fig. 2e). By analyzing the distribution of m 6 A peaks per gene, we found that most of m 6 A-moti ed mRNAs contained one or two m 6 A peak, while a small number of them contained three or more (Fig. 2f), consistent with previous studies such as ccRCC [40]. In each group of m 6 A peaks per gene, chRCC tissues always contain less gene number than normal tissues (Fig. 2f). We also counted the number of m 6 A peaks among the ve non-overlapping regions in normal and chRCC tissues and found that chRCC tissues contain less m 6 A peak numbers in 3′UTR, 5′UTR, CDS, and stop codon compared with normal tissues (Fig. 2g). These results suggest that m 6 A modi cation level decreased in chRCC tissues. Differentially methylated m 6 A genes (DMMGs) participate in multi-cancer related pathways To dissect the role of m 6 A modi cation, we subsequently identi ed differentially methylated m 6 A peaks (DMMPs) between normal and chRCC tissues using the MeTDiff R package software [46] (p £ 0.05). We identi ed 644 hypermethylated m 6 A peaks representing 593 transcripts and 1,304 hypomethylated m 6 A peaks representing 1,137 transcripts in chRCC tissues of all three subjects compared with normal tissues (Fig. 3a). Recall that our qPCR results in six subjects of chRCC showing that m 6 A writer subunit WTAP and m 6 A erasers FTO and ALKBH5 were signi cantly downregulated in chRCC (Fig. 1c), the identi ed hypomethylated m 6 A and hypermethylated m 6 A sites in chRCC could be directly induced by the aberrant expression of WTAP and FTO/ALKBH5, respectively. The identi ed hyper-and hypomethylated m 6 A peaks in chRCC tissues were respectively regarded as hyper and hypo group. Motif search analysis using HOMER revealed one overrepensented motif GGACH (H=U>C/A) in hypermethylated m 6 A peaks and two highly enriched motifs GGAC and GAACU in hypomethylated m 6 A peaks; all these three identi ed motifs resemble the canonical m 6 A motif RRACH sequence (Additional le 2: Fig. S3).
We performed metagene pro ling to examine the distribution of DMMPs within transcriptomes and found that both hyper-and hypomethylated m 6 A peaks were highly enriched around the stop codon (Fig. 3b).
Further examination of m 6 A fraction in the ve non-overlapping segments of transcripts revealed that the m 6 A peaks in both hypo and hyper groups were dominantly enriched within 3′UTR (54.6% for hypo and 56.68% for hyper), CDS (25.92% and 23.14%) and around stop codon (13.19% and 11.02%) ( Fig. 3c and Additional le 2: Fig. S2a). The majority of hypo-and hypermethylated transcripts were mRNAs (80.44% for hypo and 70.61% for hyper), and ~14-18% were lncRNA and the rest were other types of RNAs ( Fig. 3d and Additional le 2: Fig. S2b).
The hypomethylated m 6 A peak number in chRCC is 2-fold more than the hypermethylated peaks (Fig. 3a), in line with the nding that total identi ed m 6 A peaks in chRCC are less than those in normal tissues (Fig.   2a). Therefore, we rstly focused on the hypomethylated m 6 A peaks in chRCC. To explore the potential role of hypo-methylated m 6 A peaks in chRCC, we took advantage of the algorithm DAVID to examine the most preferential expression tissues of m 6 A hypomethylated genes. The results showed that m 6 A hypomethylated genes were preferentially expressed in epithelium, followed by brain, placenta, and renal cell carcinoma (RCC) (Fig. 3e), indicating the correlation between these genes and RCC. We performed Gene Ontology (GO) enrichment analysis to uncover the functions of these genes. The results revealed that m 6 A hypomethylated genes were enriched in many biological processes involved in kidney development and cancer pathogenesis, including transcription, androgen receptor signaling pathway, GTPase activity, and cell-cell adhesion (Fig. 3f). Pathway analysis showed that m 6 A hypomethylated genes were mainly enriched in cancer-related pathways (Fig. 3g). These results suggested that the hypomethylated m 6 A genes may participate in various pathophysiologic aspects of chRCC through different pathways.
We also explored the potential function of hypermethylated m 6 A peaks in chRCC. The results showed that m 6 A hypermethylated genes were preferentially expressed in brain, followed by epithelium, duodenum, fetal kidney, and ovary (Additional le 2: Fig. S2c). GO biological process analysis revealed that m 6 A hypermethylated genes were signi cantly associated with protein phosphorylation, positive regulation of cholesterol e ux, ubiquinone biosynthetic process, regulation of mitophagy, and so on (Additional le 2: Fig. S2d). Pathway analysis showed that m 6 A hypermethylated genes were mainly enriched in ubiquitin mediated proteolysis, metabolic pathways, and adherent junction (Additional le 2: Fig. S2e). Collectively, the results reveal both m 6 A hypomethylated and hypermethylated genes are involved in many regulatory pathways, especially hypomethylated genes directly enriched in cancer pathways, indicating that dysregulation of m 6 A could be a regulatory factor in the pathogenesis of chRCC.
The expression dysregulated genes in chRCC impair the normal functions of kidney We next investigated the global mRNA expression patterns in normal and chRCC tissues by using the RNA-seq dataset (m 6 A-SEAL-seq input library). The results showed that a total of 3,911 mRNAs were signi cantly dysregulated in chRCC of three subjects compared with normal tissues, including 2,344 downregulated mRNAs and 1,567 upregulated mRNAs (fold change ≥ 2, p < 0.05) (Fig. 4a). Hierarchical clustering depicted differential expression pro les in all the samples. (Fig. 4b).
We examined the preferentially expressed tissues of the dysregulated genes (3,911) using the algorithm DAVID. The result showed that the dysregulated genes were preferentially expressed in kidney, followed by liver and plasma (Fig. 4c), suggesting our data was reliable and these genes may participate in kidney development. We further performed GO analysis and KEGG pathway analysis. GO analysis revealed that the dysregulated genes were signi cantly enriched in metabolic process, transmembrane transport including sodium ion transport, oxidation-reduction process, excretion, kidney development, angiogenesis, and P450 pathway (Fig. 4d). In line with the result of GO, the KEGG analysis result also revealed that the dysregulated genes were signi cantly associated with multiple metabolic pathways including many amino acid metabolism or degradation, fatty acid degradation, cytochrome P450-related drug metabolism and xenobiotics metabolism (Fig. 4e). These results showed that the expression levels of around four thousand genes were dysregulated in chRCC, preliminarily illustrating that the dysregulated genes in chRCC impair the normal function of kidney, especially metabolic function.

New m 6 A regulatory signature in the pathogenesis of chRCC
Considering that m 6 A can either destabilize m 6 A-modi ed transcripts through the recognition of YTHDF2 or stabilize m 6 A-modi ed transcripts through the recognition of IGF2BP [14], we investigated the correlation between m 6 A methylation levels and transcript levels. We overlapped the m 6 A hypermethylated and hypomethylated genes with the differentially expressed genes. In 593 hypermethylated genes, 44 genes were downregulated and 68 genes were upregulated in chRCC tissues (Fig. 5a left). In 1,137 hypomethylated genes, 123 genes were downregulated and 51 genes were upregulated (Fig. 5a right). We next took the m 6 A hypermethylated and hypomethylated genes as two groups to analyze their transcript accumulation in chRCC and normal tissues. The result showed that the m 6 A hypermethylated genes (ie. transcripts with higher m 6 A levels) tended to preferentially exhibit upregulated transcription levels in chRCC (Fig. 5b), revealing the positive correlation between m 6 A methylation levels and transcript levels in chRCC. Note we found that the expression of YTHDF2 transcript was downregulated in chRCC (Fig. 1c). The positive correlation between m 6 A levels and transcript levels in chRCC suggests m 6 A in chRCC tends to affect gene regulation positively, potentially stabilizing m 6 A-modi ed transcripts or escaping from YTHDF2-mediated mRNA decay pathway due to the lower expression YTHDF2 in chRCC.
To further analyze the role of DMMGs in cancers, we intersected DMMGs with Cancer Gene Census (CGC) database [50], a database consists of genes with strong indications of a role in cancer, and found that 73 and 23 genes were annotated in the hypomethylated and hypermethylated genes separately (Additional le 1: Dataset S2 and Dataset S3). Among these cancer-related genes, 10 genes were differentially expressed genes, including NOCH1 and FGFR1. NOTCH1 plays distinct roles in different cancers: NOTCH1 functions as a tumor suppressor gene in mouse skin and oral squamous cell carcinoma (OSCC) [51], and the expression of NOTCH1 is decreased signi cantly in these cancer tissues [52]; whereas in Glioma and colon cancer, NOTCH1 functions as an oncogene and its expression level is increased in these two cancer tissues [53][54][55]. According to Integrative Genomics Viewer (IGV) software, the m 6 A modi cation level in NOTCH1 transcript was decreased signi cantly in chRCC compared to normal tissues (Fig. 5c left) and qPCR results in six patients showed the transcript expression level of NOTCH1 was reduced signi cantly in chRCC tissue (Fig. 5c right). The downregulated expression of NOTCH1 in chRCC is consistent with its expression pattern in mouse skin and OSCC cancers, suggesting NOTCH1 might act as a tumor suppressor in chRCC. Fibroblast growth factor receptor 1 (FGFR1) is a known oncogene. In breast cancer, the expression level of FGFR1 shows positive correlation with the ampli cation of cancer cell [56][57][58]. In lung cancer models, activation of FGFR1 promotes proliferation and migration of tumor cell while inhibition of FGFR1 suppresses tumor growth [59]. We found the m 6 A methylation level of FGFR1 was signi cantly decreased and the transcript expression level of FGFR1 was increased in chRCC tissue (Fig. 5d right), which strongly suggests FGFR1 is an oncogene in chRCC similar to that in breast cancer and lung cancer.
We conducted principal component analysis of the differential expressed DMMGs in 65 chRCC cases from The Cancer Genome Atlas (TCGA) database. Based on the expression of these genes, we could completely distinguish chRCC samples from normal samples (Fig. 5e). Cox regression screen and least absolute shrinkage and selection operator (LASSO) identi ed three m 6 A-dependent signatures (Additional le 1: Dataset S4) and de ned a m 6 A-dependent cox model in the training set (Additional le 2: Fig.   S4a). Concordance index (C-index = 0.96) showed that the proposed model has a high prognostic power. In this model, we separated patients into high-risk or low-risk group according to their risk score, and patients with different 5-years survival could be distinguished completely between the two groups (Logrank p < 0.0001) (Fig. 5f). AUC of ROC curve also con rmed the prognostic power of the m 6 A-dependent model (Additional le 2: Fig. S4b). Then the proposed model was applied to the testing set for prediction.
We calculated risk score of each patient in the testing set and assigned to high-risk or low-risk group according to the cut off value in the training set. The Kaplan-Meier survival curve and log-rank between the two groups showed signi cant difference (Fig. 5g), which demonstrated the high predictive ability of the m 6 A-dependent model. Furthermore, we also test the m 6 A-dependent model in ccRCC samples (a cohort of 602 cases from TCGA database) and PRCC dataset (a cohort of 318 cases from TCGA database). The results indicated that the m 6 A-dependent model is also suitable for ccRCC, but not PRCC (Additional le 2: Fig. S4c, d).

Discussion
Although the incidence and mortality of chRCC is the least among the three types of kidney cancer, there are still nearly thousands of people died owing to chRCC all of the world according to GCO statistic in 2020. Thereby, the underlying mechanism for adjuvant therapy is needed. Over the past few years, an increasing amount of efforts have been tried to illustrate the mechanism of m 6 A modi cation in RCC, leading to extensive accumulation about the correlation between m 6 A modi cation and RCC. However, these researches mainly focused on ccRCC and RPCC, while chRCC were rarely studied. Here, we demonstrated that m 6 A writer WTAP and m 6 A erasers FTO and ALKBH5 were downregulated in chRCC tissues. Consistently, our m 6 A-SEAL-seq in three subjects of chRCC and corresponding tumor-adjacent normal tissues identi ed 644 hypermethylated m 6 A peaks representing 593 transcripts and 1,304 hypomethylated m 6 A peaks representing 1,137 transcripts in chRCC tissues, which could be directly induced by the aberrant expression of WTAP and FTO/ALKBH5. Further functional studies showed that genes with hyper-or hypomethylated peaks were mainly enriched in kidney development and cancer pathogenesis related pathway, which is a further proof of the fundamental role of m 6 A modi cation in chRCC. We also found m 6 A reader YTHDF2 were downregulated in chRCC tissues. YTHDF2 is a wellstudied m 6 A reader and functions in m 6 A-dependent gene regulation by affecting RNA stability [12,[60][61][62][63], suggesting there is an m 6 A-dependent RNA degradation and gene dysregulation in chRCC.
Cumulative fraction of RNA transcript accumulation between the m 6 A hypermethylated genes and hypomethylated genes revealed the positive correlation between m 6 A methylation levels and transcript levels in chRCC, which may be caused by m 6 A-mediated mRNA stabilization or escaping from YTHDF2mediated mRNA decay pathway due to the lower expression YTHDF2 in chRCC.
Our RNA-seq results identi ed 2,344 downregulated mRNAs and 1,567 upregulated mRNAs. The GO and KEGG analysis of these dysregulated genes revealed that chRCC impairs the normal function of kidney, especially in metabolism. By CGC database analysis [50], we found 96 genes with differential m 6 A methylation levels (DMMGs) were causally implicated in cancers, including 73 hypomethylated genes and 23 hypermethylated genes in chRCC. Among them, ZEB1 as an oncogene and PBRM1, ASXL2 and SETD2 as tumor suppressor genes. Combined with our RNA-seq data, we found 10 cancer-related DMMGs were differential expressed in chRCC, including NOTCH1 and FGFR1. The m 6 A methylation level at 3′UTR of NOTCH1 and the transcript expression level of NOTCH1 were signi cantly reduced in chRCC. The m 6 A methylation level at 3′UTR of FGFR1 was also decreased, but the expression level of FGFR1 was increased signi cantly. Considering previous functions of NOTCH1 and FGFR1 in other cancers, our results suggest that NOTCH1 and FGFR1 might respectively act as a tumor suppressor and an oncogene in chRCC. Base our preliminary results, we proposed a potential m 6 A regulatory role in pathogenesis of chRCC, where the downregulated m 6 A writer subunit WTAP in chRCC reduces m 6 A levels in NOTCH1 and FGFR1, leading to the decreased NOTCH1 expression and the upregulated FGFR1 expression through m 6 A-mediated post-transcriptional gene regulation. Further functional studies are needed to clarify the molecular mechanisms of above-mentioned genes in the development of chRCC.
Based on the signi cant prognostic signatures identi ed from the cox regression and LASSO analysis in training set, we build a m 6 A-dependent prognostic model. With benign performance in testing data, this prognostic model shows excellent predictive ability of the survival outcome, even in ccRCC. These results indicated the m 6 A-dependent prognostic model could be used to comprehensively evaluate the survival outcome of chRCC patients and ccRCC patients in clinical practice.

Conclusion
This study presented the rst m 6 A transcriptome-wide pro le of human chRCC, which provide a potential link between abnormal m 6 A RNA modi cations and cancer-related gene expressions. We hope it can be helpful for the m 6 A methylation-based research on chRCC epitranscriptomic etiology and pathogenesis.

Declarations
Ethics approval and consent to participate The study has been approved by the ethics committee of the Peking University Third Hospital, and each participant has signed written informed consent. Hematoxylin and eosin (HE) staining and relative expression level of known m6A-related genes in normal and chRCC tissues by qPCR. a HE staining in normal tissues. b HE staining in chRCC tissues. c Relative expression level of METLL3, METTL14, WTAP, YTHDF1, YTHDF2, YTHDF3, FTO and ALKBH5. Data are presented as means ± SE, n = 6 biological replicates. *p < 0.05 (two-sided), * * * p < 0.005 (two-sided).  Pie chart presenting the fraction of the con dent hypomethylated m6A peaks across ve non-overlapping transcript segments. d Pie chart presenting RNA types (that is, transcript species) of the con dent hypomethylated m6A peaks identi ed in chRCC. e-g UP_TISSUE (e), Gene ontology (GO) (f) and KEGG (g) analysis of the hypomethylated m6A genes.