Expression of m6A RNA Methylation Regulators and Their Clinical Predictive Value in Intrahepatic Cholangiocarcinoma

Background : N6-methyladenosine (m 6 A) RNA methylation regulators have been implicated in the carcinogenesis and progress of a variety of cancers. Until now, the effects of them on intrahepatic cholangiocarcinoma (ICC) have been poorly understood. Methods : We used the GEO databases to systematically evaluate the expression profiles of 36 m 6 A RNA methylation regulators in ICC patients and produced a signature to assess its prognostic values. In vitro experiments were implemented to confirm the expression level. Results : Compared to normal intrahepatic bile duct tissues, more than half of these 36 genes showed different levels of expression in ICC tissues. Two groups emerged from the consensus cluster analysis of these 36 genes. The two cluster of patients had significantly different clinical outcomes. In addition, we created a m 6 A-related prognostic signature that performed exceptionally well in the prognostic categorization of ICC patients, based on the ROC curves, Kaplan-Meier curves, and univariate and multivariate Cox regression analyses. Further research showed that there was a significant association between the m 6 A-related signature and the manifestations of tumor immune microenvironment in ICC. The expression level and biological effect of METTL16, one of the two m 6 A RNA methylation regulators incorporated in the signature, were confirmed and explored by using in vitro experiments. Conclusions : This analysis revealed the predictive roles of m 6 A RNA methylation regulators in ICC.


Introduction
Intrahepatic cholangiocarcinoma (ICC) is one of the two most common subtypes of primary liver cancer, second only to hepatocellular carcinoma (HCC) [1,2].Worldwide, the incidence of ICC is gradually increasing [3].The different subtypes of liver cancer show obvious heterogeneity in tumor biological behavior.Studies have shown that ICC tumors are significantly more malignant than other types of HCC [4].Tumor recurrence and malignant progression are higher in ICC patients, even with conventional treatments such as surgical resection, local interventional therapy, liver transplantation and chemo-radiotherapy [5].Substantial progress in exploring the genetic pattern of ICC have been made [6], but beneficial treatment options are still insufficient.Therefore, predictive biomarkers or therapeutic targets are of vital importance for ICC patients.
Epigenetics, known as RNA modification, has recently been introduced in clinical practice.It has been widely reported that N6-methyladenosine (m 6 A) plays an crucial biological function as one of the most universal RNA modification methods in eukaryocytes [7,8], especially when it comes to the regulation of gene expression [9,10].m 6 A RNA methylation is a invertible and dynamic process conducted by m 6 A RNA methylation regulators, which include "writers", "erasers" and "readers", that is, methyltransferases, demethylases and binding proteins.A growing body of data suggests that m 6 A mutations have a role in the biological activity of a variety of malignancies [11,12].The clinical role of m 6 A RNA methylation regulator in many tumors has also been established, including gastric cancer [13], head and thyroid carcinoma [14], head and neck squamous cell carcinoma [15], breast cancer [16], and hepatocellular carcinoma [17].
However, few studies have looked at the prognostic significance of m 6 A RNA methylation regulators in the patients of ICC.The Gene Expression Omnibus (GEO) is a public repository of genomic data uploaded by researchers worldwide, and is available for free download.We used GEO data to comprehensively analyzed the relevance between the expression levels and clinical outcomes of 36 m 6 A RNA methylation regulators in ICC.Furthermore, we constructed a m 6 A-related prognostic signature that has the potential to predict survival in ICC patients.

Data Source
Data was draw from 30 ICC patients and 27 normal intrahepatic bile duct tissues with data of RNA-seq transcriptome and matching clinical characteristics in the GSE107943 from GEO dataset.The general clinical characteristics of the ICC patients who were included in the study are summarized in Table 1.In the GEO datasets, 36 m 6 A RNA methylation regulators with attainable expression data were discovered (Supplementary Material 1) [18,19].The expression profiles of these 36 genes were compared in ICC and normal intrahepatic bile duct tissues.

Consensus Clustering
Consensus Clustering is widely used in cancer classification.The overall survival rates (OS) related m 6 A RNA methylation regulators were screened using univariate Cox regression analysis.The STRING database (https: //cn.string-db.org/) and the R package "corrplot" were used to investigate interactions.Based on the expression levels of these 36 genes, the R package "ConsensusClusterPlus" was used to explore clear and appropriate clusters of ICC patients [20].The package "survminer" in R software (version 4.1.2,Lucent, New Jersey, USA) was utilized to plot the Kaplan-Meier curve to compare the OS and progression free survival (PFS) of the different clusters.

Development of the m 6 A-Related Prognostic Signature
The least absolute contraction and selection operator (LASSO) Cox regression was used to study all 36 m 6 A RNA methylation regulators [21].Two m 6 A RNA methylation regulators determined by LASSO coefficients to effectively calculate risk scores.The calculation formula is: risk score = S (βi × Expi) (i = the number of m 6 A RNA regulators).The median of risk score was employed as a cut-off to divide into the high-/low-risk groups.And the R software package "timeROC" was used to draw the ROC curves, which assess the accuracy of predicting survival of the m 6 A-based signature.The uni-and multi-variate Cox regression analyses were performed to study the effects of the signature in ICC patients.

Nomogram Development
The nomogram is a common tool for tumor prognosis assessment.The typical practice is to first screen the biological characteristics and clinical indicators of patients to build a prognosis model, and then to visualize the model.It can facilitate clinical decision making.A nomogram was created with the R packages "regplot", "rms" and "Hmisc".The signature and clinical factors (age, sex, AJCC stage, grade, levels of CEA and CA19-9, tumor size, vascular invasion) were covered to construct the nomogram to evaluate the survival of 1-, 3-, and 5-year OS for ICC patients.

Evaluation of the Immune Landscape
Gene Ontology (GO) and Kyoto Encyclopedia of genes and Genomes (KEGG) are two major databases for gene function and structure enrichment.To functionally annotate the differentially expressed genes, the GO analysis and the KEGG analysis were performed in different groups of signature risk scores.Then we used Immune Cell Abundance Identifier (ImmuCellAI) to assess the infiltrates of 24 immune cells in different groups [22].

Cell Culture
The two ICC cell lines, HCCC-9810 and RBE cells was obtained from Institute of Biochemistry and Cell Biology, Shanghai, China.All newly purchased cell lines have been undergone mycoplasma testing and cell identification of Short Tandem Repeat (STR) at the first time.They were cultured for subsequent in vitro validation experiments.The culture conditions were set to 37 °C and 5% CO 2 concentration.The cell lines were transfected with plasmids using duo transfection reagent.Shanghai GeneChem (Shanghai, China) provided vectors expressing targeted short hairpin RNA (shRNA) sequences.The targeting sequences for METTL16 control and the non-targeting disruptive RNA sequence are 5′-AGGGAGTAAACTCACGAAATCCT-3′ (shMETTL16) and 5′-TTCTCCGAACGTGTCACGT-3′ (shNC), respectively.

CCK8 and EdU Cell Proliferation Assay
In the CCK8 assay, cells were seeded into 96-well plates at 2 × 10 3 cells/well.On days 1, 2, 3, 4, and 5, 10 µL of CCK8 reagent (Dojindo, Kumamoto Ken, Japan) were applied to each well after the cells had adhered, and the specific absorbance was determined using spectrophotometry at 450 nm.In the EdU assay, ICC cells knock-down METTL16 and the corresponding control cells were seeded at a rate of 1.5 × 10 5 cells/well on a 24-well culture plate and incubated for 24 hours.To analyze and assess cell proliferation, the EdU Staining Proliferation Kit (abcam, iFluor 647) was utilized, and the particular step was carried out in accordance with the manufacturer's instructions.A Leica fluorescence microscope was used to examine the samples.

Western Blotting Assay
The protein lysis buffer containing the protease inhibitor was added to the cell suspension to extract the total protein.After protein concentration quantification, gel electrophoresis was uesd to separate the proteins that were rapidly transferred to a PVDF membrane.After blocking, the membranes were incubated with the primary antibody anti-METTL16 (Sigma-Aldrich, St. Louis, USA) overnight at 4 °C.The next day, we used a secondary antibody to incubate with the membranes for 1.5 hour at room temperature.Finally, Enhanced Chemiluminescence (ECL) chemiluminescence was used to demonstrate the immune response.

Statistical Analysis
The R software (version 4.1.2) was utilized in this study.When comparing continuous variables with normal distribution, the student t-test or one-way ANOVA test was used in their corresponding pairs.The Wilcoxon test was used to deteremine the gene expression levels in different subgroups.The associations between the m 6 A RNA methylation regulators and immunological checkpoints were calculated using the "Spearman" method.We defined statistically significant as two-tailed p < 0.05.

Expression Profile of m 6 A RNA Methylation Regulators in ICC
Firstly, a heatmap and violin plot were created to determine the expression levels of 36 m 6 A RNA reg- ulators between ICC and normal intrahepatic bile duct tissues.Compared with normal tissues, 10 upregulated genes and 13 downregulated genes were confirmed in ICC tissues.Among these, 36 m 6 A RNA methylation regulators, METTL14, ZC3H13, ALKBH5, TRMT112, SETD2, SRSF10, XRN1, FMR1, YTHDC1, YTHDC2, YTHDF2 and YTHDF3 were down-regulated in ICC, and RBM15B, METTL16, CPSF6, RBMX, IGF2BP1, IGF2BP2, IGF2BP3, SRSF3, EIF3B and EIF3H were upregulated in ICC (Fig. 1).The protein-protein interaction (PPI) network and co-expression analysis were implemented to demonstrate the relationships between 36 m 6 A RNA methylation regulators."Writers" interacted with a wide spectrum of additional m 6 A RNA methylation regulators, but "erasers" interacted with fewer regulators (Fig. 2A).The majority of m 6 A RNA methylation regulators were favorably connected to others, and negative coexpression interactions were underrepresented (Fig. 2B).

Two Clusters of ICC Patients
ICC patients was divided into different subgroups using Consensus Clustering based on K-means via the Con-sensusClusterPlus package.Consistency clustering verifies the rationality of clustering (i.e., finding a suitable K value) through resampling based method.A common criterion is to select K value with small CDF descent slope.But the specific determination of K value can also be judged according to cumulative distribution function (CDF) and consensus matrix.Here, we selected k = 2 as the most appropriate value to categorize ICC patients into two clusters.Consensus matrix was established with the value of [0,1] (Fig. 3A-C).The majority of m 6 A RNA methylation regu-lators had higher expression levels in cluster 2 than in cluster 1 (Supplementary Material 2).The survival analysis showed that compared with cluster 1, ICC patients in cluster 2 had shorter OS (p = 0.017) and PFS (p = 0.0074) which were statistically significant (Fig. 3D,E).

Construction of a m 6 A-Related Prognostic Signature
To initially uncover the key m 6 A RNA regulators related to OS in ICC patients, we first performed univariate Cox regression analysis.WTAP (p = 0.035), CBLL1 (p = 0.037), METTL16 (p = 0.004), FTO (p = 0.011), IGF2BP2 (p = 0.003), SRSF10 (p = 0.021), NXF1 (p = 0.032) were shown to be strongly associated with OS.WTAP, IGF2BP2 and SRSF10 were identified as risk genes with HR greater than 1.CBLL1, METTL16, FTO, and NXF1 were protective genes with HR less than 1 (Fig. 4A and Supplementary Material 3).A m 6 A-based signature was created using the LASSO Cox regression analysis that ultimately included two m 6 A RNA regulators, METTL16 and FTO.According to the above calculation formula, that is, risk score = S (βi × Expi) (i = the number of m 6 A RNA regulators).The risk score for ICC patients was obtained by calculation: risk score = (-0.3417× METTL16) + (-0.0158 × FTO).We selected the median as the cut-off point and classified ICC patients into two groups, low-/high-risk (Fig. 4B,C).As was depicted by the curves of Kaplan-Meier analysis, those in the high-risk group had decreased OS and PFS when compared with patients in the low-risk group (Fig. 4D,E).The area under the ROC curve validated the precision of this m 6 A-related prognostic signal in predicting OS.In the ICC patients, the AUCs of OS were 0.950, far better than other clinical factors (Fig. 4F).

The m 6 A Based Signature in the Prognosis of ICC Patients
As is shown in the univariate Cox analysis, CEA (p = 0.01), AJCC Stage (p = 0.006), vascular invasion (p = 0.018), and risk score (p = 0.004) were significantly linked with the OS of ICC patients (Fig. 5A,B).Cox regression analysis indicated that CEA was also significantly corelated with OS (p = 0.010 and 0.042, respectively).Therefore, to develop a practical clinical tool for predicting the prognosis of ICC patients, we established a nomogram based on CEA with the risk score based on the m 6 A-related signature.The prognostic nomogram was proved to accurately predict the 1-, 3-, and 5-year OS of ICC patients (Fig. 5C).

GO and KEGG Signaling Pathway and Immune Landscape of the m 6 A-based Signature
To explore the mechanism for these pathways, 110 differentially expressed genes (DEGs) were screened out in two different risk groups and the functional annotation of these DEGs was performed by the KEGG and GO pathway analyses.As is shown in Fig. 6A, DEGs were involved in the regulation of the inflammatory, humoral immune, and acute inflammatory responses, and the lipid transport and protein activation cascade, which are immunity-related biological processes (Fig. 6A).The consistent results showed that inflammation-related pathways were significantly enriched in the KEGG pathway analysis, including the complement and coagulation cascade, the PPAR signaling pathway, retinol metabolism, and mineral absorption (Fig. 6B).
In view of these results, we continued to explore the association between the m 6 A-based signature and the immunophenotype of the ICC.As is shown in Fig. 6C, the risk score was significantly related to the expression of immune checkpoints, covering PD-1, PD-L1, B7H3, LAG-3, CTLA-4, IDO1 and TIM-3.The following trends appeared in the samples included in this study.The highrisk group tended to have varying degrees of immune cell infiltration when compared with the low-risk group, with lower abundance of CD4+ T cells, CD8+ T cells, Tfh cells, CD8 naive cells, and MALT cells but a higher abundance of macrophages and neutrophils (Fig. 6D,E).Immune landscape of the m 6 A-based signature need more specific and detailed experiments to describe.

Validation of the Expression of METTL16 by in Vitro Experiments
Based on the above bioinformatics analysis, we obtained a prognostic signature based on FTO and METTL16.We found that the expression of FTO and METTL16 was inversely proportional to the prognosis of ICC patients.It has been reported that down-regulation of FTO might build a gene network that facilitated ICC progression [23].Therefore, METTL16 is regarded as an important candidate gene for further research.After transfecting HCCC-9810 and RBE cells with METL16-specific shRNA, we first verified the transfection efficiency (Fig. 7A), and then explored the effect of knocking out METTL16 on ICC proliferation.As is shown in the CCK8 assay, compared with control cells, down-regulation of METTL16 can considerably promote cell proliferation (Fig. 7B,C).The BrdU assay also confirmed the inhibiting role of METTL16 in the proliferation of HCC cells (Fig. 7D,E).Therefore, METTL16 may emerge as a treatment for ICC in the future.We continue to conduct research on the specific and in-depth mechanism of METTL16 in ICC.

Discussion
In this research, we used GEO database to construct a m 6 A-based prognostic signature in ICC after rigorously evaluating the expressions of 36 m 6 A RNA methylation regulators.In 23 of them, the expression was altered in ICC patients, and WTAP, CBLL1, METTL16, FTO, IGF2BP2, SRSF10, and NXF1 were found to be significantly related with the prognosis.Using consensus cluster analysis, ICC patients were allocated into two cohorts based on 36 m6A RNA regulators.Clusters 1/2 were shown to strongly link with the OS of ICC patients in the Kaplan-Meier analysis.We eventually screened two m 6 A RNA regulators (METTL16 and FTO) through Lasso regression and incorporated them into a signature.High risk scores were observed to have a positive association with poor OS in ICC.Accordingly, the signature of m 6 A RNA regulators could be regarded as a new biomarker to evaluate the outcome of ICC patients.It was demonstrated that this signature can be an independent prognostic marker for ICC in both univariate and multivariate analyses.Furthermore, the outcome of ICC patients can be monitored more accurately by using the nomogram which was based on the m 6 A based signature and clinical factors.Through GO and KEGG pathway ICC is the second most common intrahepatic malignancy.ICC is characterized by a high degree of malignancy, unknown cause and nonspecific clinical presentation.There are no obvious clinical symptoms in the early stage of ICC, and most patients have lost the opportunity to operate when they are found.Therefore, higher requirements are put forward for the early diagnosis and timely treatment of ICC.At present, there are still insufficient studies on the clinical prognosis of ICC patients, and the exploration of specific mechanisms is extremely desired.The pathogenesis of ICC is a complicated process, which involves numerous abnormal gene expression profiles in diverse signaling pathways [4].There are currently no effective therapies for ICC patients, especially for those with advance stage or metastatic disease.Furthermore, biomarkers that can robustly predict the outcomes of ICC patients are scarce.Given the significance of m 6 A RNA methylation regulators in the development of cancers, it is necessary to explore m 6 A RNA regulator-based biomarkers for ICC prevention and therapy.In the m 6 A modification process, m 6 A methyltransferase, as a "writer", is an important family of molecules [10].So far, two types of RNA methyltransferases have been found, including METTL16 which is capable of modifying RNA with m 6 A [24].The impact of METTL16 in malignant tumors is now being studied in hepatocellular carcinoma and gastric cancer.Xiao-Kun Wang et al. [25] found that METTL16-mediated m 6 A methyla- tion increased the proliferation of gastric cancer via increasing cyclin D1 expression.Pei Wang et al. [26] revealed that down-regulation of METTL16 was an independent risk factor for disease free survival and suggests poor OS and disease free survival in patients with hepatocellular carcinoma.However, the role of METTL16 in ICC remains unknown.
Currently, METTL16 is regarded as a type of "writer" that does not depend on the METTL3/METTL14 complex to regulate pre-mRNA alternative splicing RNA stability and protein translation efficiency [27].METTL3 prefers to methylate RNAs with the RRACH motif (H = A, C or U; R = A or G), while substrates catalyzed by METTL16 usually has the unique motif UACAGAGAA [28].METTL16 differs from METTL3 in that it may methylate RNA substrate without interacting with other proteins or components.The function of METTL3 is dependent on forming a methyltransferase complex with ZC3H13, METTL14, VIRMA, RBM15 and WTAP.U6 snRNA, MAT2A mRNA and lncRNA MALAT1 and XIST are currently recognized METTL16 substrates [29].It has been shown that METTL16 recognizes its RNA substrate via a mixture of sequence and structure, and METTL16 is also involved in the pre-mRNA splicing process.Therefore, METTL16 is both the "writer" and the "reader" of m 6 A [30].Our study is a preliminary exploration of the role of METTL16 in the progress of ICC and its possible biological mechanism, but specific verification is still needed, which is the direction of our future research.
Our study still has some limitations.First, the specific small number of patients included may vibrate the reliability of the conclusions, although we have included ICC data with the largest amount of patients that have detailed genetic sequencing and follow-up data.Second, the GEO database lacks crucial information regarding patient prognosis, such as centralized pathologic review, quality of the surgery, details of chemoradiotherapy, and further treatment of complications.Third, our verification results are only in vitro, and the specific mechanism is still not very clear.Therefore, multi-center RCT studies with long-term follow-up are needed to check on the conclusions in our study, and we continue to collect ICC specimens for further research.

Conclusions
We constructed and verified a signature based on m 6 A RNA methylation regulators, which can be applied as an independent biomarker for evaluating the prognosis of ICC patients.Additional biochemical experiments and multicenter clinical trials are needed to verify our conclusions.

Fig. 1 .
Fig. 1.The expression profiles of 36 m 6 A RNA methylation regulators between ICC and normal intrahepatic bile duct tissues.(A) The heatmap revealed the expression levels of these 36 correlated genes.(B) The violin plot revealed the expression comparison of these 36 correlated genes between ICC and the control samples.*p < 0.05, **p < 0.01, ***p < 0.005.

Fig. 3 .
Fig. 3. Consensus clustering of ICC patients.(A) When k = 2, we divided ICC patients into two distinct clusters.(B) Cumulative distribution function of consensus clustering for k = 2 to 9. (C) Relative change in area under CDF curve for k = 2 to 9. (D,E) Kaplan-Meier curve for OS and PFS of ICC patients in cluster 1/2.

Fig. 5 .
Fig. 5.The signature based on m 6 A in ICC.(A,B) Univariate and multivariate Cox regression analysis on the OS of ICC patients.(C) The nomogram was made based on different characteristics.

Fig. 6 .
Fig. 6.GO and KEGG Signaling Pathway and Immune Landscape of the m 6 A-based Signature.(A) GO biological processes of DEGs in two different risk cohorts.(B) KEGG pathway analysis in two different risk cohorts.(C) Correlation between risk score and immune checkpoints.(D) The heatmap showed the abundance of immune cells.(E) The difference of abundance of immune cells in different groups.

Fig. 7 .
Fig. 7.The expression levels and the role in proliferation of METTL16 by in Vitro Experiments.(A) The expression of METTL16 in hcc-9810 and RBE cells was significantly reduced after transfection with METTL16 -specific shRNA specific shRNA.(B,C) In the CCK8 assay, knocking out METTL16 promoted cell proliferation.(D,E) In the BrdU assay, knocking out METL16 promoted cell proliferation both in the hcc-9810 and RBE cell lines.*p < 0.05, **p < 0.01, ***p < 0.005.