Analysis of Multi-Layer RNA Modification Patterns for the Characterization of Tumor Immune Microenvironment in Hepatocellular Carcinoma
Jiyuan Xing1,2 
Shen Shen1,2 
Zihui Dong1 Xiaobo Hu1,2 Lixia Xu1,2 Xiaorui Liu1,2 Qinggang Li1,2 Yize Zhang1 
Gangying Cui1,2 
Zujiang Yu1,2*- 1Gene Hospital of Henan Province, Precision Medicine Center, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China
- 2Department of Infectious Diseases, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China
Background: RNA modifications have emerged as important posttranscriptional changes in multiple tumor cellular processes and tumorigenesis, including hepatocellular carcinoma (HCC). However, the potential roles and the interaction between regulators of RNA modifications and the tumor microenvironment (TME) are unclear in HCC.
Methods: The gene expression profiles of 26 RNA modification “writers” were investigated in the TCGA cohort. The unsupervised clustering approach was used to class these RNA modification regulators. The characteristics of immune cell infiltration from TME for each cluster was tested by the CIBERSORT method. Additionally, we established a scoring model to evaluate the RNA modification characteristics of individual tumors. The associations between the scoring model and genetic as well as clinical characteristics, drug sensitivity, and response to immunotherapy were also analyzed.
Results: We mapped the somatic mutations and somatic copy number variation of the RNA modification regulators. The expression of all selected regulators was detected, and two modification patterns were identified that featured distinct immune cell infiltration characteristics. Subsequently, we developed a score model (termed as WM-Score model). Furthermore, the survival analysis showed that the WM-Score value was associated with HCC patient prognosis. The results of the ROC curves analysis and multivariate analysis all confirmed that the WM-Score value was strongly associated with anti-cancer drug resistance and therapeutic efficacy of immunotherapy, thus could be used as an independent risk factor in HCC.
Conclusion: Our research identified two RNA modification patterns characterized by distinct TME, and the WM-Score model was developed that might serve as reliable prognostic and immunotherapeutic effect predictor of HCC.
Introduction
In 2018, hepatocellular carcinoma (HCC) was predicted to be the sixth most prevalent cancer worldwide (Bray et al., 2018), with a 5 years survival rate as low as 9.1%, and an overall median survival of 9 months (Giannini et al., 2015). Infection by HBV or HCV, chronic alcohol consumption, and obesity-related NASH are the principal causes of HCC (Llovet et al., 2021). The condition is usually diagnosed at an advanced stage, therefore, effective treatments for advanced metastatic HCC are limited. Although there are surgical and chemotherapy options, the mortality rate of HCC remains high. Forms of immunotherapy, such as immune checkpoint inhibitors (ICIs) have been used to capture the disease progression and to enhance adaptive immunity in advanced HCC (Ou et al., 2020). Meanwhile, only a subset of patients show therapeutic response to ICIs, and this response it is difficult to predict. Therefore, a deeper understanding of the molecular mechanism of HCC is necessary to improve patient survival.
Recently, RNA modifications, coined the “epitranscriptome”, have emerged as crucial posttranscriptional regulators of the gene expression process (Barbieri and Kouzarides, 2020). Increasing evidence has revealed that these modifications have huge implications for human pathophysiology, including cancer (Frye et al., 2016; Jonkhout et al., 2017; Nachtergaele and He, 2017; Ontiveros et al., 2019). Accordingly, over 170 different types of chemical modifications of cellular RNAs have been described, among which methylation modifications account for two-thirds and are widely present in various RNA types (Barbieri and Kouzarides, 2020). The most abundant and better characterized internal RNA modification is N6-methyladenosine (m6A) that regulates multiple aspects of RNA metabolism, such as RNA processing, RNA translation, and nuclear export (Roundtree et al., 2017; Sun et al., 2019). N1-methyladenosine (m1A) is an important post-transcriptional RNA modification that has been found in tRNA, rRNA, mitochondrial RNA and mRNA (RajBhandary et al., 1966; Peifer et al., 2013; Li et al., 2017; Safra et al., 2017). APA is an RNA-processing mechanism that generates distinct 3′ termini on mRNAs and other RNA polymerase II transcripts (Tian and Manley, 2017). RNA editing mediated by adenosine deaminase acting on RNA enzymes a well-documented post-transcriptional mechanism altering nucleotide in selected transcripts (Nishikura, 2010). RNA modification is catalyzed by RNA methyltransferases called “writers” (they add a specific modification), demethylases or “erasers” (they remove a specific modification), and m6A-binding proteins or “readers” (they recognize and bind modified nucleotides). The RNA modification is a dynamic process, and the interaction between each type of methylation modification has not yet been fully elucidated (Davalos et al., 2018; Xue et al., 2020; Nombela et al., 2021).
Accumulating evidence supports the prominent role of the complex and diverse tumor immune microenvironment (TIME), including cancer cells, locally infiltrating immune cells, stromal cells, and active medium, in tumor cell proliferation, invasion, and metastasis (Azambuja et al., 2019; Fu et al., 2019). Non-malignant cells are not only one of the major players of cancer progression, but also determine the immunotherapeutic response (Lu et al., 2019). Therefore, a comprehensive analysis of the diversity of TME and different immune phenotypes can guide and improve immunotherapeutic responsiveness (Binnewies et al., 2018).
In this study, we focused on the most heavily modified RNA types, including m6A, alternative polyadenylation (APA), m1A, and A-to-I RNA editing. Furthermore, we comprehensively analyzed the correlation between various types of RNA modification regulators and cell-infiltrating characteristics of TIME by integrating the genomic and transcriptomic alterations of samples from The Cancer Genome Atlas - Liver Hepatocellular Carcinoma (TCGA-LIHC) databases. Two distinct modification patterns with different immune cell characteristics were identified. In addition, we developed the WM-Score model to quantify the efficacy of “writers” in modifying individual tumors and to predict the prognosis and immunotherapeutic response of HCC patients.
Methods
Data Acquisition and Processing
The gene expression profiles and clinical annotations were downloaded from the Cancer Genome Atlas (TCGA) portal (http://cancergenome.nih.gov/). Data cohorts with missing information were removed. A total of 356 cases of TCGA-LIHC were used for further analysis. The R Bioconductor package and R (version 3.6.2) were employed for data analysis.
Drug sensitivity data were collected from The Genomics of Drug Sensitivity in Cancer (GDSC) database (www.cancerRxgene.org) (Yang et al., 2013). Spearman’s correlation analysis was utilized to evaluate the association between the scoring model and drug reaction, where |Rs| > 0.2, and FDR <0.05 was considered significant correlation.
The immunotherapy dataset IMvigor210 cohort was used to explore the immunotherapy response and prognosis of HCC patients with different WM-Score values. The standardized RNA-sequencing data of 1111 HCC patients with detailed clinicopathological data were downloaded from http://research-pub.gene.com. The data were analyzed using the IMvigor210CoreBiologies R package.
Unsupervised Clustering Analysis
In order to explore the robust clustering of HCC cases, we employed the unsupervised clustering approach to analyze the gene profiles of RNA modification writers. A total of 26 RNA modification regulators, including seven m6A modification enzymes (KIAA1429, METTL14, ZC3H13, METTL3, WTAP, RBM15B, and RBM15), 12 APA modification enzymes (CPSF1, CPSF2, CPSF3, CPSF4, CSTF1, CSTF2, CSTF3, CF1, PCF11, CLP1, NUDT21, and PABPN1), four m1A modification enzymes (TRMT10C, TRMT6, TRMT61A, and TRMT61B), and three A-I modification enzymes (ADARB1, ADARB2 and ADAR) were analyzed. An NMF-based consistent clustering algorithm was used to determine RNA modification patterns based on the mRNA expression of analyzed regulators. Unsupervised cluster analysis was performed by The Consensus Cluster Plus package as previously described (Wilkerson and Hayes, 2010).
Gene Set Variation Analysis (GSVA)
GSVA is a gene set enrichment method that provides increased power to estimate changes of subtle pathway activity over a sample population in an unsupervised manner (Hänzelmann et al., 2013). We conducted GSVA analysis to explore the association between RNA modifications and biological processes. The gene set “h.all.v7.2” and “c2.cp.kegg.v7.1”were derived from the MSigDB database (Zhu et al., 2020). The functional annotation of 26 “writer” genes was conducted by the clusterProfiler R package, with a cutoff value of FDR <0.05. An adjusted P with value < 0.05 was considered as indicative of statistical significance.
Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT)
CIBERSORT is a method that can accurately estimate the fraction of diverse cell subsets in gene expression profiles from complex tissues (http://cibersort.stanford.edu) (Newman et al., 2015). To predict the immune subset composition of HCC samples from gene expression profiles, CIBERSORT was used to estimate the relative abundance of 22 types of immune cells (model = absolute, permutation = 1,000, disable quantile normalization for RNA-Seq data as recommended).
Construction of the WM-Score Scoring System
Firstly, the RNA modification-related differentially expressed genes (DEGs) among distinct RNA modification clusters were collected using “limma” package of R software. Next, we performed univariate cox regression model to analyze the correlation of each gene with overall survival, and the significant prognosis DEGs were used for further analysis. Subsequently, distinct genomic subtypes were determined by unsupervised clustering analyses. In addition, the prognostic analysis was performed for each genomic subtype and extract principal component 1 and 2 as the signature scores. Finally, the RNA modification score was defined using a method similar to that used in analyzing gene-gene interactions (GGIs): WM-Score = Σ (PC1i + PC2i), which is defined as the expression of final RNA modification phenotype-related genes (Sotiriou et al., 2006; Zeng et al., 2019).
Statistical Analysis
A Wilcoxon rank-sum test was utilized to compare differences between two groups, and Kruskal-Wallis test was used for comparisons of multiple groups. The discrimination accuracy of the WM-Score model was described by receiver operating characteristic (ROC) analysis. Kaplan-Meier method estimate curves were generated for prognostic analysis, and the differences between groups were evaluated by a log-rank test. Univariate and multivariate analyses were further carried out to assess independent risk factors. All data were analyzed by the R 4.0.1 software. A two-tailed p < 0.05 was considered as statistically significant.
Results
Landscape of Genetic Alterations of 26 RNA Modification “Writers”
A total of 26 RNA modification “writers” were selected in this study, which included seven m6A modification “writers”, three A-I modification “writers”, 12 APA modification “writers”, and four m1A modification “writers” (Supplementary Table S1) (Li et al., 2016; Tang S. J. et al., 2020; Shen et al., 2021). To explore the genetic alterations in RNA modification writers, we examined the incidence of somatic mutations and somatic copy number variation (CNV) for all “writers” based on the TCGA database. Among 356 samples from TCGA-LIHC, 42 (11.8%) exhibited genetic changes of these writers, and the details was shown in the Figure 1A. The highest mutation frequency was presented in CPSF1, followed by ADARB2 and KIAA1429 (Figure 1A), while METTL3, METTL14, TRMT61A, TRMT61B, CSTF3, and NUDT2 did not show any mutations in tumor samples. Next, we used the hallmark gene set to perform gene set variation analysis (GAVA) to compare the mutation groups and those without mutation in “writers”. The GSVA indicated significantly enriched carcinogenic activation pathways in the mutation group, such as those of E2F targets, G2M checkpoint, MYC, and MTORC1 signaling pathway (Figure 1B). Furthermore, the investigation of CNV alteration in 26 regulators showed that ADAR, CPSF1, CPSF4, TRMT10C and KIAA1429 had a widespread frequency of CNV gain, while ZC3H13, CF1, METTL14, NUDT21, and WTAP had a significant CNV loss (Figure 1C). To explore whether the above CNV alterations affected the expression of the 26 RNA modification regulators, we compared the expression level of these regulators between tumor samples and paired normal samples. The results showed increased mRNA levels of most “writers” in tumor samples in comparison to normal samples (Figure 1D), suggesting that CNV might be the major factor leading to the aberrant expression of medication regulators. Notably, the mRNA levels of some “writers” were increased, while the frequencies of CNV loss for those were high. Therefore, further investigations were performed. According to the CNV value, patients were divided into 3 groups, including CNV amplification group, CNV deletion group, and normal group, and the mRNA expression of “writers” were compared between these groups (Figure 1E). The results showed mostly elevated expression for the group of patients with CNV amplification compared with the other groups with CNV deletion or normal CNV in these “writers”. Taken together, we mapped the genetic alterations of the 26 RNA modification “writers” between control tissues and tumor tissues, suggesting that these changes might play vital functions in HCC tumorigenesis and progression.
FIGURE 1. Expression pattern of 26 RNA modification “writer” genes in TCGA-LIHC (A) Frequency of mutations of the 26 identified regulator genes (B) Gene set variation analysis (GAVA) was used to compare the regulator mutation group and the non-mutation group (C) The CNV mutation frequency of 26 regulator genes in TCGA-LIHC (D) The expression of 26 RNA modification regulator genes between tumor tissues and control tissues (F) The mRNA expression of “writer” among three groups, including amplification group, CNV deletion group, and normal group.
The RNA Modification Patterns Are Characterized by Distinct TIME Cell Infiltration Characteristics
In order to further understand the role of RNA modification “writers” in HCC, we performed univariate analysis of the 26 regulators based on the TCGA-LIHC cohort. We found that 16 of 26 “writers” were markedly correlated with the OS of HCC patients (Figure 2A). Next, we explored the relationship among “writers” and found that most were positively or negatively correlated with each other (Figure 2B). Thus, it is suspected that the crosstalk between different “writers” may have a vital function in the different modification patterns of HCC.
FIGURE 2. RNA methylation modification pattern and related biological pathways (A) The univariate cox regression analysis shows that 16 of 26 regulators are associated with of patient prognosis in the TCGA-LIHC cohort (B) Heatmap of the Spearman’s correlation analysis presenting negative (blue) and a positive (red) correlation among the “writers” in HCC (C) Unsupervised cluster analysis of 26 “writers” in HCC. Blue represents low expression of “writers” genes; red represents high expression of these genes (D) GSVA enrichment analysis of KEGG pathway chances between cluster 1 and 2. Blue indicate not activated pathways, and red indicates activated pathways.
We performed consensus clustering to classify patients into distinct RNA modification patterns based on the mRNA expression of “writers” (Supplementary Table S2). Eventually, two RNA modification patterns with 204 cases were determined in pattern 1 (cluster 1), and 113 cases in pattern 2 (cluster 2) (Figure 2C). Subsequently, “GSVA” enrichment analysis was employed to further understand the biological behaviors between the distinct two clusters. Our results indicated that cluster 1 was significantly enriched in metabolism and drug metabolism pathways, such as sulfur metabolism, primary bile acid biosynthesis, tyrosine metabolism, tryptophan metabolism, drug metabolism cytochrome P450, drug metabolism other enzymes, renin angiotensin system, while cluster 2 enrichment pathways were mainly linked to proliferation and signal transduction, including cell cycle, calcium conduction, etc. (Figure 2D).
Emerging evidence suggests that RNA modifications interact with the tumorigenic environment, thus affecting tumor occurrence, development, and prognosis (Jiang et al., 2020; Chen H. et al., 2021; Chong et al., 2021). Therefore, the function of the RNA methylations in the TME were further explored. The association analysis using the CIBERSORT method revealed that the identified RNA modification regulators might have close links with immune cell infiltration from the TME (Figure 3A). For instance, METTL14, ZC3H13, CSTF3, and ADAR were markedly negatively associated with Mo macrophage differentiation, while their positive association was observed with METTL3, RBM15B, KIAA1429, TRMT61A, TRMT6, CPSF1, and NUDT21. Moreover, we analyzed the difference in immune cell infiltration from TME between cluster 1 and cluster 2. The results revealed that the infiltration of M2, T cells, mast cells, and monocytes was higher in cluster 1. Notably, though, the infiltration of M1, regulatory T cells and follicular helper T cells was higher in cluster 2 (Figures 3B,C). Overall, cluster 2 was usually enriched in immunosuppressive cells, indicating a poor prognosis, whereas cluster 1 was characterized by immune cell activity, indicating a beneficial prognosis. These findings suggested that RNA modification “writers” play crucial roles in immune cell infiltration and TME formation.
FIGURE 3. Tumor immune microenvironment characterization of the RNA modification patterns (A) The correlation between the 26 “writers” and TME in HCC was analyzed by CIBERSORT (B) The different characters of immune cell infiltration between clusters. Log(FC) > 0 represents that the immune cells were enriched in cluster 1 (C) The types of immune cells between distinct RNA modification patterns.
Generation of RNA Modification Signature Model
Our results above demonstrated the important role of RNA modification in TME formation and patient prognosis, while these findings were based on RNA modification patterns and could not accurately evaluate the capacity of the RNA modification as a prognostic predictor in individual HCC patients. The underlying genetic alterations in these two RNA modification patterns were still unclear. Based on these queries, we examined the transcriptional expression change between the two patterns. A total of 273 DEGs related to RNA modification patterns were identified, and the further enrichment analysis showed that these DEGs were enriched in many essential biological processes, including DNA-binding transcription activator activity, signaling receptor activator activity, and multicellular organismal response to stress (Figures 4A–C). Subsequently, according to unsupervised clustering analysis based on the 273 DEGs, patients were classified into two stable transcriptomic subtypes: cluster A and cluster B (Figure 4D), with 242 and 75 of the 317 HCC patients, respectively. The prognosis of patients in cluster B was poorer than those in gene cluster A (Figure 4E; p < 0.0001, log-rank test).
FIGURE 4. Construction of RNA modification model (A–C). GSVA enrichment analysis revealed the DNA transcription signaling pathways. Stress reception signaling pathways (A) and signal transduction activation signaling pathways (B) were correlated with (C) 273 DEGs between cluster 1 and cluster 2 (D). Unsupervised clustering of the 273 DEGs to identify two genomic subtypes (E) Survival analysis showing the poor prognosis of patients in cluster A group compared with those in cluster B group (p < 0.0001, Log-rank test) (F). The score of cluster 1 was higher than that of cluster 2. G. The score of cluster A was significantly higher than that of cluster B (H). The difference of immune cell infiltration abundance between WM-Score groups calculated by the CIBERSORT algorithm.
Furthermore, we developed a score model based on the DEGs between gene clusters. As described in the Methods section, a scoring model named writers of RNA modification-score (WM-Score) was constructed. We discovered that cluster 1 had a higher WM-Score value than cluster 2 (Figure 4F). Consistently with this, cluster A also showed a higher score value than cluster B (Figure 4G). To evaluate the association of WM-Score value with TME, we further calculated the abundance of immune cell infiltration for the low and the high WM-Score value groups. We found that the infiltration rate of M0 macrophages, monocytes, and TfCD8 was higher in the high WM-Score value group, and that of activated NK-activated cells and M1 macrophages was higher in the low WM-Score value group (Figure 4H).
Association Between WM-Score and Clinical Characteristics
After confirming the efficacy of the WM-Score model in predicting patient prognosis, we investigated whether this model could be applied to determine the tumorigenesis, progression, invasion and metastasis of HCC. The prognostic efficiency of the scoring model was explored through classifying patients into low and high score groups the using “survminer” package. As expected, patients with high score demonstrated a poorer prognosis than those with low score in the TCGA-LIHC cohort (Figure 5A). We used ROC curve analysis to determine the discrimination accuracy of the scoring model in predicting patient prognosis. The area under the ROC curves (AUCs) of WM-Score values were 0.84, 0.76 and 0.79 at 1, 3 and 5 years overall survival, respectively (Figure 5B). Multivariate analysis for the TCGA-LIHC cohort also demonstrated that the WM-Score could serve as an independent prognostic predictor in HCC (Figure 5C). All of these results indicated that the WM-Score model has accurate prognostic value for HCC patients. The analysis of difference in WM-Scores between different TNM grades and clinical grades in the TCGA database indicated that samples with higher clinical grades and TNM stages usually have higher WM-Score values (Figures 5D,E). In addition, considering the EMT-related pathways, the samples with different WM-Score value had different pathway characteristics. For the TCGA database, samples with high WM-Score value were significantly related to cell cycle, DNA damage repair, and DNA replication, while samples with low WM-Score value were related to EMT, WNT target, and cell cycle regulators (Figure 5F).
FIGURE 5. The clinical characteristics and prognosis of HCC correlated with the WM-Score model. (A). Kaplan-Meier overall survival for HCC patients in the high and low WM-Score groups. (B). The predictive accuracy of the WM-Score model in the TCGA-LIHC cohort (AUC: 0.84, 076, and 0.79, corresponding to 365, 1,095, and 1825 days OS, respectively). (C). Multivariate cox regression analysis of factors, which included WM-Score, stage-T, stage, stage-M, gender, stage-N, patient age, and grade in the TCGA-LIHC cohorts. (E). WM-Score differences among grade and stage-T of HCC in TCGA-LIHC. (F). Heatmap showing the GSVA score of EMT signaling pathways between different WM-Score groups in TCGA-LIHC.
Value of WM-Score Model in Chemotherapy and Therapy Sensitivity
For several years, sorafenib has been approved a treatment option for advanced HCC patients, while efficacy of sorafenib is limited by drug resistance (Gnoni et al., 2019). Aiming to further investigate whether the WM-Score value affected drug sensitivity, we evaluated the correlation between the scoring model and the drug response of tumor cell lines. Using Spearman’s correlation analysis, 15 significant correlation pairs were identified in the Cancer Drug Sensitivity Genomics (GDSC) database between scoring model and drug reaction (Yang et al., 2013). Among them, eight pairs of drug sensitivity were related to WM-Score value, and seven pairs showed resistance related to WM-Score value (Figure 6A). In addition, we also analyzed the signaling pathways of these drugs to determine target genes. We found that drugs associated with high WM-Score value mainly target KIT, CLAP, and cell cycle signaling pathways. In contrast, drugs related to low WM-Score value mostly target apoptosis regulation and cell cycle signaling pathways (Figure 6B). Taken together, these findings indicate that the WM-Score values are related to drug reaction, and thus might offer a framework to guide the treatment strategy of HCC.
FIGURE 6. Correlation of scoring model with drug reaction and immunotherapeutic response. (A). Spearman’s analysis was used to determine the correlation between score and drug response in GDSC. (B). The association between drugs and targeted signaling pathways. (C). The difference in the score between immunotherapeutic responder in the IMvigor210 cohort. (D). The percentage of patients with different responses (including SD, PD, CR, and PR) to PD-L1 blockade immunotherapy. E-G. Total samples, or Stage I + II samples, or Stage III + IV samples in the IMvigor210 cohort all showed a significant difference in survival between samples with high and low WM-Scores based on survival analysis.
Role of WM-Score Model in Predicting anti-PD-1/L1 Immunotherapy
In recent years, immune checkpoint inhibitors (ICIs) have made breakthroughs in the treatment of advanced HCC, while biomarkers that could effectively predict the efficacy of immunotherapy are still lacking. Herein, we explored whether the WM-Score model could predict therapeutic response to ICI therapy in HCC patients. For the IMvigor210 cohort, the therapeutic efficacy was significantly better in patients with low WM-Score value compared to those with high WM-Score value (Figure 6C). The frequency of response to anti-PD-1/L1 treatment in the low WM-Score value group was higher than that in the high WM-Score value group (Figure 6D). We also analyzed the survival difference of all samples of IMvigor210 and those under different stages. The results showed that total samples (Figure 6E), or Stage I + II samples (Figure 6F), or Stage III + IV samples (Figure 6G) all exhibited a marked difference in survival between samples with high and low WM-Score value. Especially in the prediction of high-stage clinical samples, the WM-Score value demonstrated extremely high power. Collectively, our results proved that the WM-Score model might serve as a potential predictor of response to anti-PD-1/L1 immunotherapy.
Discussion
A growing pool of evidence indicates that RNA modifications play a key role in gene expression, whose disruption impacts the pathogenesis of human disease, including cancer (Frye et al., 2016). Although RNA modifications as genetic or epigenetic alterations of genes are not traditionally considered as cancer drivers, cumulative evidence suggests that abnormal RNA modifications are functionally correlated with many hallmarks of cancer, such as proliferation, invasion, migration, differentiation, self-renewal, and response to therapy (Cui et al., 2017; Weng et al., 2018; Jin et al., 2019).
For instance, N6-methyladenosine (m6A) is an RNA methylation that is the most abundant form of internal mRNA modification. Yang et al. reported the involvement of the m6A modification in the 3′-UTR of oncogene CDCP1 mRNA in bladder cancer cell growth and progression (Yang et al., 2019). Lang et al. indicated that the m6A modification showed an important function in regulating the stability of viral transcripts and EBV-mediated tumorigenesis (Lang et al., 2019). Furthermore, Lan et al. reported that m6A methyltransferase KIAA1429 was high expressed in HCC tissues and knockdown KIAA1429 inhibited cell proliferation and metastasis in vitro and in vivo (Lan et al., 2019). Chen et al. found the writer CPSF1 of APA was significantly increased in HCC tissues and associated with poor survival outcomes (Chen S.-l. et al., 2021). All these studies focused on one or two modification regulators to explore their dysregulation, function, and underlying mechanism in cancer, however, the deposition of RNA modifications is a dynamic process involving multiple modification regulators. In the present study, we comprehensively described the molecular and biological features of different regulators of RNA modifications and identified two distinct RNA modification subtypes based on multiple modification regulators. Importantly, the two subtypes (cluster 1 and cluster 2) are not only associated with clinical survival, but also with the abundance of immune cell infiltration.
Considering the diversity and complexity of TME, the thorough understanding of its implications in cancer is a significant challenge. In recent years, some research groups have documented that RNA modifications were closely associated with TME. Shen et al. attempted to explore the role of m6A regulators in HCC immune cell infiltration and prognosis, and identified three m6A subtypes based on TCGA and GEO database, which were related to three known immune phenotypes (including immune-inflamed phenotype, immune-excluded phenotype, and immune-desert phenotype) (Shen et al., 2021). Chong et al. also discovered three m6A modification patterns among 1,370 colon cancer cases, which were correlated with different outcomes and TME characterization (Chong et al., 2021). Three m6A modification patterns with distinct TME cell-infiltrating characteristics were also determined in gastric cancer (Zhang et al., 2020), lung adenocarcinoma (Li et al., 2020), pancreatic adenocarcinoma (Tang R. et al., 2020), and gliomas (Xu et al., 2020). Similar with our analysis, these studies were based on a large number of samples in the subject database, such as TCGA and GEO, in order to clarify the role of modification in tumor immune regulation and progression. In our study, we further identified two stable transcriptomic subtypes based on the DGEs of the two RNA modification clusters. Especially, the transcriptomic subtypes were significantly associated with the immune cell activation and prognosis of HCC patients. Thus, the systematical evaluation of RNA modification patterns provides novel clues for understanding TME characterization in HCC. Gu et al. found 3 m 5C regulator-mediated methylation modification patterns based on the expression of 13 m 5C regulators which were closely associated with different immune cell infiltration characteristics in HCC (Gu et al., 2021). Shen et al. demonstrated three m6A modification patterns which affect tumor immune infiltrates and prognosis of patients with HCC (Shen et al., 2021). Previous studies mainly centered upon one types of RNA modification to explore their effect on TME. Here, we performed a comprehensive analysis of multiple types of RNA modification and highlights the cross-talk and the roles of RNA modifications in the TME and response to immunotherapy.We developed the WM-Score model to accurately predict the prognostic value of the RNA modification in individual patients. We found that this model could be applied to assessing clinicopathological features, such as clinical grades and TNM grades, and patients with higher clinical grades and TNM grades usually had higher WM-Score value.
In addition, the RNA modification pattern with higher WM-Score value tended to correlate with immune cell suppression in the tumor microenvironment, while the pattern with lower WM-Score value was usually associated with immune activation. In the IMvigor210 cohort, WM-Scores model was found to be linked with immune cell infiltration in TME as well as response to anti-PD-1/L1 immunotherapy, suggesting the application potential of WM-Score model for predicting HCC anti-PD-1/L1 immunotherapy.
Conclusion
In the present work, the RNA modification regulators were comprehensively analyzed, and the correlation was demonstrated between RNA modification patterns and cell-infiltrating characteristics in the TME. The systematic evaluation of individual tumor RNA modification pattern might serve as a useful predictor of prognosis for HCC patients and act as a valuable tool for developing more effective immunotherapy strategies.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material further inquiries can be directed to the corresponding author.
Author Contributions
ZY and JX designed the study, JX, SS, and ZD wrote this manuscript XH, LX, XL, and QL searched the articles and made figures YZ and GC collected and analyzed the data All authors read and approved the final manuscript.
Funding
This work was supported by Henan Medical Science and Technology Joint Building Program (LHGJ20200320 and LHGJ20200387), the Science and Technology Research Project of Henan Province (202102310115), and National Natural Science Foundation of China (82070643 and U1904164).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We thank the patients and investigators who participated in TCGA and other datasets for providing data.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.761391/full#supplementary-material
References
Azambuja, J. H., Ludwig, N., Braganhol, E., and Whiteside, T. L. (2019). Inhibition of the Adenosinergic Pathway in Cancer Rejuvenates Innate and Adaptive Immunity. Ijms 20, 5698. doi:10.3390/ijms20225698
Barbieri, I., and Kouzarides, T. (2020). Role of RNA Modifications in Cancer. Nat. Rev. Cancer 20, 303–322. doi:10.1038/s41568-020-0253-2
Binnewies, M., Roberts, E. W., Kersten, K., Chan, V., Fearon, D. F., Merad, M., et al. (2018). Understanding the Tumor Immune Microenvironment (TIME) for Effective Therapy. Nat. Med. 24, 541–550. doi:10.1038/s41591-018-0014-x
Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., and Jemal, A. (2018). Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA: A Cancer J. Clinicians 68, 394–424. doi:10.3322/caac.21492
Chen, H., Yao, J., Bao, R., Dong, Y., Zhang, T., Du, Y., et al. (2021a). Cross-talk of Four Types of RNA Modification Writers Defines Tumor Microenvironment and Pharmacogenomic Landscape in Colorectal Cancer. Mol. Cancer 20, 29. doi:10.1186/s12943-021-01322-w
Chen, S.-l., Zhu, Z.-x., Yang, X., Liu, L.-l., He, Y.-f., Yang, M.-m., et al. (2021b). Cleavage and Polyadenylation Specific Factor 1 Promotes Tumor Progression via Alternative Polyadenylation and Splicing in Hepatocellular Carcinoma. Front. Cel Dev. Biol. 9, 616835. doi:10.3389/fcell.2021.616835
Chong, W., Shang, L., Liu, J., Fang, Z., Du, F., Wu, H., et al. (2021). m6A Regulator-Based Methylation Modification Patterns Characterized by Distinct Tumor Microenvironment Immune Profiles in colon Cancer. Theranostics 11, 2201–2217. doi:10.7150/thno.52717
Cui, Q., Shi, H., Ye, P., Li, L., Qu, Q., Sun, G., et al. (2017). m 6 A RNA Methylation Regulates the Self-Renewal and Tumorigenesis of Glioblastoma Stem Cells. Cel Rep. 18, 2622–2634. doi:10.1016/j.celrep.2017.02.059
Davalos, V., Blanco, S., and Esteller, M. (2018). SnapShot: Messenger RNA Modifications. Cell 174, 498. e491. doi:10.1016/j.cell.2018.06.046
Frye, M., Jaffrey, S. R., Pan, T., Rechavi, G., and Suzuki, T. (2016). RNA Modifications: what Have We Learned and where Are We Headed. Nat. Rev. Genet. 17, 365–372. doi:10.1038/nrg.2016.47
Fu, Y., Liu, S., Zeng, S., and Shen, H. (2019). From Bench to Bed: the Tumor Immune Microenvironment and Current Immunotherapeutic Strategies for Hepatocellular Carcinoma. J. Exp. Clin. Cancer Res. 38, 396. doi:10.1186/s13046-019-1396-4
Giannini, E. G., Farinati, F., Ciccarese, F., Pecorelli, A., Rapaccini, G. L., Di Marco, M., et al. (2015). Prognosis of Untreated Hepatocellular Carcinoma. Hepatology 61, 184–190. doi:10.1002/hep.27443
Gnoni, A., Licchetta, A., Memeo, R., Argentiero, A., Solimando, A. G., Longo, V., et al. (2019). Role of BRAF in Hepatocellular Carcinoma: A Rationale for Future Targeted Cancer Therapies. Medicina 55, 754. doi:10.3390/medicina55120754
Gu, X., Zhou, H., Chu, Q., Zheng, Q., Wang, J., and Zhu, H. (2021). Uncovering the Association between m5C Regulator-Mediated Methylation Modification Patterns and Tumour Microenvironment Infiltration Characteristics in Hepatocellular Carcinoma. Front. Cel Dev. Biol. 9, 727935. doi:10.3389/fcell.2021.727935
Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7
Jiang, Y., Wan, Y., Gong, M., Zhou, S., Qiu, J., and Cheng, W. (2020). RNA Demethylase ALKBH5 Promotes Ovarian Carcinogenesis in a Simulated Tumour Microenvironment through Stimulating NF‐ κ B Pathway. J. Cel Mol Med 24, 6137–6148. doi:10.1111/jcmm.15228
Jin, D., Guo, J., Wu, Y., Du, J., Yang, L., Wang, X., et al. (2019). m6A mRNA Methylation Initiated by METTL3 Directly Promotes YAP Translation and Increases YAP Activity by Regulating the MALAT1-miR-1914-3p-YAP axis to Induce NSCLC Drug Resistance and Metastasis. J. Hematol. Oncol. 12, 135. doi:10.1186/s13045-019-0830-6
Jonkhout, N., Tran, J., Smith, M. A., Schonrock, N., Mattick, J. S., and Novoa, E. M. (2017). The RNA Modification Landscape in Human Disease. Rna 23, 1754–1769. doi:10.1261/rna.063503.117
Lan, T., Li, H., Zhang, D., Xu, L., Liu, H., Hao, X., et al. (2019). KIAA1429 Contributes to Liver Cancer Progression through N6-methyladenosine-dependent post-transcriptional Modification of GATA3. Mol. Cancer 18, 186. doi:10.1186/s12943-019-1106-z
Lang, F., Singh, R. K., Pei, Y., Zhang, S., Sun, K., and Robertson, E. S. (2019). EBV Epitranscriptome Reprogramming by METTL14 Is Critical for Viral-Associated Tumorigenesis. Plos Pathog. 15, e1007796. doi:10.1371/journal.ppat.1007796
Li, X., Xiong, X., Wang, K., Wang, L., Shu, X., Ma, S., et al. (2016). Transcriptome-wide Mapping Reveals Reversible and Dynamic N1-Methyladenosine Methylome. Nat. Chem. Biol. 12, 311–316. doi:10.1038/nchembio.2040
Li, X., Xiong, X., Zhang, M., Wang, K., Chen, Y., Zhou, J., et al. (2017). Base-Resolution Mapping Reveals Distinct m1A Methylome in Nuclear- and Mitochondrial-Encoded Transcripts. Mol. Cel 68, 993–1005. e1009. doi:10.1016/j.molcel.2017.10.019
Li, Y., Gu, J., Xu, F., Zhu, Q., Chen, Y., Ge, D., et al. (2020). Molecular Characterization, Biological Function, Tumor Microenvironment Association and Clinical Significance of m6A Regulators in Lung Adenocarcinoma. Brief Bioinform 22, bbaa225. doi:10.1093/bib/bbaa225
Llovet, J. M., Kelley, R. K., Villanueva, A., Singal, A. G., Pikarsky, E., Roayaie, S., et al. (2021). Hepatocellular Carcinoma. Nat. Rev. Dis. Primers 7, 6. doi:10.1038/s41572-020-00240-3
Lu, C., Rong, D., Zhang, B., Zheng, W., Wang, X., Chen, Z., et al. (2019). Current Perspectives on the Immunosuppressive Tumor Microenvironment in Hepatocellular Carcinoma: Challenges and Opportunities. Mol. Cancer 18, 130. doi:10.1186/s12943-019-1047-6
Nachtergaele, S., and He, C. (2017). The Emerging Biology of RNA post-transcriptional Modifications. RNA Biol. 14, 156–163. doi:10.1080/15476286.2016.1267096
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust Enumeration of Cell Subsets from Tissue Expression Profiles. Nat. Methods 12, 453–457. doi:10.1038/nmeth.3337
Nishikura, K. (2010). Functions and Regulation of RNA Editing by ADAR Deaminases. Annu. Rev. Biochem. 79, 321–349. doi:10.1146/annurev-biochem-060208-105251
Nombela, P., Miguel-López, B., and Blanco, S. (2021). The Role of m6A, m5C and Ψ RNA Modifications in Cancer: Novel Therapeutic Opportunities. Mol. Cancer 20, 18. doi:10.1186/s12943-020-01263-w
Ontiveros, R. J., Stoute, J., and Liu, K. F. (2019). The Chemical Diversity of RNA Modifications. Biochem. J. 476, 1227–1245. doi:10.1042/bcj20180445
Ou, Q., Yu, Y., Li, A., Chen, J., Yu, T., Xu, X., et al. (2020). Association of Survival and Genomic Mutation Signature with Immunotherapy in Patients with Hepatocellular Carcinoma. Ann. Transl Med. 8, 230. doi:10.21037/atm.2020.01.32
Peifer, C., Sharma, S., Watzinger, P., Lamberth, S., Kötter, P., and Entian, K.-D. (2013). Yeast Rrp8p, a Novel Methyltransferase Responsible for m1A 645 Base Modification of 25S rRNA. Nucleic Acids Res. 41, 1151–1163. doi:10.1093/nar/gks1102
RajBhandary, U. L., Stuart, A., Faulkner, R. D., Chang, S. H., and Khorana, H. G. (1966). Nucleotide Sequence Studies on Yeast Phenylalanine sRNA. Cold Spring Harbor Symposia Quantitative Biol. 31, 425–434. doi:10.1101/sqb.1966.031.01.055
Roundtree, I. A., Evans, M. E., Pan, T., and He, C. (2017). Dynamic RNA Modifications in Gene Expression Regulation. Cell 169, 1187–1200. doi:10.1016/j.cell.2017.05.045
Safra, M., Sas-Chen, A., Nir, R., Winkler, R., Nachshon, A., Bar-Yaacov, D., et al. (2017). The m1A Landscape on Cytosolic and Mitochondrial mRNA at Single-Base Resolution. Nature 551, 251–255. doi:10.1038/nature24456
Shen, S., Yan, J., Zhang, Y., Dong, Z., Xing, J., and He, Y. (2021). N6-methyladenosine (m6A)-Mediated Messenger RNA Signatures and the Tumor Immune Microenvironment Can Predict the Prognosis of Hepatocellular Carcinoma. Ann. Transl Med. 9, 59. doi:10.21037/atm-20-7396
Sotiriou, C., Wirapati, P., Loi, S., Harris, A., Fox, S., Smeds, J., et al. (2006). Gene Expression Profiling in Breast Cancer: Understanding the Molecular Basis of Histologic Grade to Improve Prognosis. J. Natl. Cancer Inst. 98, 262–272. doi:10.1093/jnci/djj052
Sun, T., Wu, R., and Ming, L. (2019). The Role of m6A RNA Methylation in Cancer. Biomed. Pharmacother. 112, 108613. doi:10.1016/j.biopha.2019.108613
Tang, R., Zhang, Y., Liang, C., Xu, J., Meng, Q., Hua, J., et al. (2020a). The Role of m6A-Related Genes in the Prognosis and Immune Microenvironment of Pancreatic Adenocarcinoma. PeerJ 8, e9602. doi:10.7717/peerj.9602
Tang, S. J., Shen, H., An, O., Hong, H., Li, J., Song, Y., et al. (2020b). Cis- and Trans-regulations of Pre-mRNA Splicing by RNA Editing Enzymes Influence Cancer Development. Nat. Commun. 11, 799. doi:10.1038/s41467-020-14621-5
Tian, B., and Manley, J. L. (2017). Alternative Polyadenylation of mRNA Precursors. Nat. Rev. Mol. Cel Biol 18, 18–30. doi:10.1038/nrm.2016.116
Weng, H., Huang, H., Wu, H., Qin, X., Zhao, B. S., Dong, L., et al. (2018). METTL14 Inhibits Hematopoietic Stem/Progenitor Differentiation and Promotes Leukemogenesis via mRNA m6A Modification. Cell Stem Cell 22, 191–205. e199. doi:10.1016/j.stem.2017.11.016
Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a Class Discovery Tool with Confidence Assessments and Item Tracking. Bioinformatics 26, 1572–1573. doi:10.1093/bioinformatics/btq170
Xu, S., Tang, L., Dai, G., Luo, C., and Liu, Z. (2020). Expression of m6A Regulators Correlated with Immune Microenvironment Predicts Therapeutic Efficacy and Prognosis in Gliomas. Front. Cel Dev. Biol. 8, 594112. doi:10.3389/fcell.2020.594112
Xue, C., Zhao, Y., and Li, L. (2020). Advances in RNA Cytosine-5 Methylation: Detection, Regulatory Mechanisms, Biological Functions and Links to Cancer. Biomark Res. 8, 43. doi:10.1186/s40364-020-00225-0
Yang, F., Jin, H., Que, B., Chao, Y., Zhang, H., Ying, X., et al. (2019). Dynamic m6A mRNA Methylation Reveals the Role of METTL3-m6A-CDCP1 Signaling axis in Chemical Carcinogenesis. Oncogene 38, 4755–4772. doi:10.1038/s41388-019-0755-0
Yang, W., Soares, J., Greninger, P., Edelman, E. J., Lightfoot, H., Forbes, S., et al. (2013). Genomics of Drug Sensitivity in Cancer (GDSC): a Resource for Therapeutic Biomarker Discovery in Cancer Cells. Nucleic Acids Res. 41, D955–D961. doi:10.1093/nar/gks1111
Zeng, D., Li, M., Zhou, R., Zhang, J., Sun, H., Shi, M., et al. (2019). Tumor Microenvironment Characterization in Gastric Cancer Identifies Prognostic and Immunotherapeutically Relevant Gene Signatures. Cancer Immunol. Res. 7, 737–750. doi:10.1158/2326-6066.Cir-18-0436
Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y. L. (2020). m6A Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric Cancer. Mol. Cancer 19, 53. doi:10.1186/s12943-020-01170-0
Zhu, S., Wu, Q., Zhang, B., Wei, H., Li, B., Shi, W., et al. (2020). Autophagy-related Gene Expression Classification Defines Three Molecular Subtypes with Distinct Clinical and Microenvironment Cell Infiltration Characteristics in colon Cancer. Int. Immunopharmacology 87, 106757. doi:10.1016/j.intimp.2020.106757
Keywords: hepatocellular carcinoma, RNA methylation modification, TME, immunotherapy, prognosis
Citation: Xing J, Shen S, Dong Z, Hu X, Xu L, Liu X, Li Q, Zhang Y, Cui G and Yu Z (2021) Analysis of Multi-Layer RNA Modification Patterns for the Characterization of Tumor Immune Microenvironment in Hepatocellular Carcinoma. Front. Cell Dev. Biol. 9:761391. doi: 10.3389/fcell.2021.761391
Received: 19 August 2021; Accepted: 28 October 2021;
Published: 10 November 2021.
Edited by:
Jiang Chen, Zhejiang University, ChinaReviewed by:
Shaohua Xu, Tongji University, ChinaAntonio Giovanni Solimando, University of Bari Aldo Moro, Italy
Gongbo Fu, Nanjing General Hospital of Nanjing Military Command, China
Jinhui Liu, Nanjing Medical University, China
Mingdi Liu, Jilin University, China
Copyright © 2021 Xing, Shen, Dong, Hu, Xu, Liu, Li, Zhang, Cui and Yu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zujiang Yu, johnyuem@zzu.edu.cn