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 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.

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 ($\beta$i $\times$ 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.

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.

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].

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.

In the CCK8 assay, cells were seeded into 96-well plates at 2 $\times$ 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 $\times$ 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.

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.

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.

Firstly, a heatmap and violin plot were created to determine the expression levels of 36 m${}^6$A RNA regulators 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 up-regulated 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 co-expression interactions were underrepresented (Fig. 2B).

Fig. 1.

The expression profiles of 36 m${}^{\mathrm{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. 2.

Interrelationships among 36 m${}^{\mathrm{6}}$A RNA methylation regulators. (A) Protein-protein interaction network. (B) Correlation of 36 m${}^6$A RNA methylation regulators.

ICC patients was divided into different subgroups using Consensus Clustering based on K-means via the ConsensusClusterPlus 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 regulators 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).

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.

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 ($\beta$i $\times$ Expi) (i = the number of m${}^6$A RNA regulators). The risk score for ICC patients was obtained by calculation: risk score = (–0.3417 $\times$ METTL16) + (–0.0158 $\times$ 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).

Fig. 4.

The m${}^{\mathrm{6}}$A-related prognostic signature. (A) Univariate Cox regression analysis of OS related m${}^6$A RNA methylation regulators with statistical significance in ICC patients. (B,C) LASSO analysis with minimal lambda value. (D,E) The Kaplan-Meier curve for OS and PFS of ICC patients in high/low risk. (F) ROC curves revealed that the AUCs of risk score, age, sex, AJCC stage, CEA, CA19-9, grade, size, and vascular invasions were 0.950, 0.511, 0.405, 0.629, 0.658, 0.726, 0.532, 0.493, 0.649.

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).

Fig. 5.

The signature based on m${}^{\mathrm{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.

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).

Fig. 6.

GO and KEGG Signaling Pathway and Immune Landscape of the m${}^{\mathrm{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.

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 high-risk 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.

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.

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.