High Expression of RRM1 Mediated by ncRNAs Correlates with Poor Prognosis and Tumor Immune Infiltration of Hepatocellular Carcinoma

Introduction

HCC ranks as the fourth leading cause of cancer death worldwide.¹ Globally, the cancer ranks second in cancer-related mortality among men.² In the United States, the incidence of HCC has been increasing and is expected to be one of the top three causes of cancer-related deaths by 2030.³ Currently, there is no effective treatment for liver cancer. Despite these encouraging advances, treatment options for patients with liver cancer are still very limited and further development of new treatment options is needed. Therefore, understanding the cellular mechanisms underlying the pathogenesis of HCC is essential for the development of effective anti-cancer therapies.

Ribonucleotide reductase M1 (RRM1), the large and catalytic subunit of ribonucleotide reductase, is an enzyme required for the conversion of ribonucleotides to deoxyribonucleotides,^4,5 and deoxyribonucleoside diphosphate is further phosphorylated to deoxyribonucleoside triphosphate, which is ultimately used in DNA synthesis and repair processes.⁶ RRM1 is overexpressed in a variety of cancers, including HCC, and correlates with tumor progression and metastasis.⁷ It is also an independent prognostic factor for HCC. A growing body of evidence suggests that RRM1 is involved in the development of HCC and may be a promising therapeutic target. For example, RRM1 is a molecular target of gemcitabine.⁸ Anti-RRM1 shRNA had significant antitumor effects on RRM1-high expressing tumors.⁶ However, comprehensive studies on the expression, prognosis and mechanisms of RRM1 are still lacking in HCC. In addition, the relevance of RRM1 to tumor immune infiltration is uncertain in HCC.

In our study, first, we carried out RRM1 expression analysis and survival analysis for a variety of human cancers. Then, regulation of ncRNA-associated RRM1 involving small RNAs (miRNAs) and long non-coding RNAs (lncRNAs) has been explored in HCC. Finally, we clarified the relationship between RRM1 expression and immune cell infiltration, immune cell biomarkers or immune checkpoints in HCC. Overall, these results show that high expression of RRM1 mediated by ncRNAs correlates with poor prognosis and tumor immune infiltration of HCC.

Materials and Methods

Human Cancer Data

All available data of 21 cancer types were retrieved from the TCGA data (RNAseq data in HTSeq-FPKM format in TCGA-ALL project), and then the data were transformed into transcripts per million (TPM). Finally, Mann Whitney U-test was performed on the data using R package ggplot2. p < 0.05 was set as statistically significant. This study complies with the publication guidelines and access rules of TCGA.

GEPIA

GEPIA is a powerful database containing high-throughput RNA sequencing data,⁹ which is often used for gene expression correlation analysis. We applied it to determine the expression of RRM1 and lncRNAs in a variety of human cancers. Results were considered statistically significant when p < 0.05. We used GEPIA to analyze the survival of RRM1 in 8 different cancers and the prognostic value of possible lncRNAs in HCC. A value of Log rank p < 0.05 was considered statistically significant. Furthermore, the correlation between RRM1 expression and immune checkpoints in HCC was also assessed by GEPIA. The following conditions were considered statistically significant: |R| > 0.1 and p value < 0.05.

Prediction of Candidate miRNA

The prediction of upstream possible miRNAs is performed by multiple target gene prediction programs, including PITA, RNA22, miRmap, microT, miRanda, PicTar and TargetScan. Only miRNAs predicted by two or more programs were selected for later analysis. Thus, we obtained possible miRNAs of RRM1.

starBase

starBase is a powerful database for pan-cancer analysis of RNA-RNA and RBP-RNA interactions, as well as for survival and differential expression analysis of miRNA, lncRNA, pseudogenes and mRNA. We used Starbase v3.0 to predict candidate lncRNAs that may bind to miR-22-3p and miR-101-3p. In addition, starBase was used to analyze the expression correlation for miRNA-RRM1, MIR4435-2HG or SNHG6 in HCC.

Kaplan-Meier Survival Curves

The Kaplan-Meier Plotter is a tool that used to predict survival of tumor patients, and finally draw survival maps.¹⁰ In this study, we performed survival analysis of ncRNAs (lncRNAs and miRNAs) in patients with HCC. All results (the best cut-off value, log-rank P value, HR, and 95% CIs) were automatically analyzed and presented by the website. Results were considered statistically significant when p < 0.05.

TIMER

TIMER database is a web-based tool for comprehensive analysis and visualization of immune cell infiltration in tumors. We used it to generate scatter plots to visualize the correlation of RRM1 expression with tumor-infiltrating or immune checkpoint expression in HCC. Results were considered statistically significant when p < 0.05.

Results

Expression Analysis of RRM1 in Human Pan-Cancer

To determine whether RRM1 is associated with tumorigenesis, we first performed its expression analysis in 21 human tumors using the TCGA database (Figure 1A). Our results indicated that RRM1 expression was significantly increased in 14 cancer types, including BRCA, CESC, CHOL, COAD, ESCA, GBM, HNSC, KIRC, LIHC, LUAD, LUSC, STAD, THCA and UCEC, but significantly decreased only in KICH. At the same time, no significant differences in the expression of RRM1 in BLCA, KIRP, PAAD, PCPG, PRAD or READ were found. Next, we further carried out RRM1 expression analysis using GEPIA database, and found that RRM1 expression in BRCA, CESC, COAD, CHOL, GBM, LUSC, LIHC, PAAD, READ and STAD was statistically increased when compared with normal tissues (Figure 1B–K). In addition, no significant differences in the expression of RRM1 in remaining cancers were found. Altogether, the expression of RRM1 was increased in eight tumors: BRCA, CESC, COAD, CHOL, GBM, LUSC, LIHC and STAD, suggesting that RRM1 may function as a key promoter in the development of these 8 cancers.

Figure 1 Pan-cancer analysis of RRM1 expression across 21 cancer types using TCGA database (A). (B–K) RRM1 expression analysis in BRCA (B), CESC (C), CHOL (D), COAD (E), GBM (F), LIHC (G), LUSC (H), PAAD (I), READ (J) and STAD (K) compared with TCGA and GTEx normal tissues. *p value < 0.05; **p value < 0.01; ***p value < 0.001. Abbreviation: ns, no significance.

Prognostic Value of RRM1 in Human Cancers

Here, we used GEPIA database to perform survival analysis for RRM1 in various tumors, including BRCA, CESC, COAD, CHOL, GBM, LUSC, LIHC and STAD. For these 8 cancers, we used two indicators, disease-free survival (RFS) and overall survival (OS), to evaluate prognostic analysis. In OS analysis, high RRM1 expression in COAD predicted a better prognosis, however LIHC patients with high expression of RRM1 demonstrated an unfavorable prognosis (Figure 2). In RFS analysis, only high expression of RRM1 in LIHC was related to poor prognosis (Figure 3). RRM1 was not found to be statistically relevant in predicting the prognosis of patients with other cancers. Combined with the prognostic analysis of OS and RFS, RRM1 can be used as a biomarker for HCC patients.

Figure 2 The OS analysis of RRM1 in a variety of human cancers evaluated by GEPIA database. (A–H) The OS curves of RRM1 in BRCA (A), CESC (B), CHOL (C), COAD (D), GBM (E), LIHC (F), LUSC (G), and STAD (H).

Figure 3 The disease-free survival (RFS) analysis of RRM1 in a variety of human cancers evaluated by GEPIA. (A–H) The RFS curves of RRM1 in BRCA (A), CESC (B), CHOL (C), COAD (D), GBM (E), LIHC (F), LUSC (G), and STAD (H).

Prediction and Analysis of Upstream miRNAs of RRM1

It has been repeatedly demonstrated that many ncRNAs regulate gene expression at the transcriptional and post transcriptional levels. We first predicted the upstream miRNAs that might bind to RRM1 to determine whether RRM1 is regulated by some upstream ncRNAs, and finally found 13 candidate miRNAs. We further used Cytoscape software to visualize the miRNA-RRM1 regulatory network (Figure 4A). Then, we conducted correlation analysis between RRM1 and candidate miRNAs. According to the negative regulation mechanism of miRNA, RRM1 should be negatively correlated with miRNAs. As shown in Figure 4B, RRM1 was significantly negatively related to miR-22-3p and miR-101-3p, and positively related to miR-138-5p, miR-216b-5p, miR-330-3p, miR-1323, miR-212-3p, miR-132-3p, miR-577 and miR-340-5p in HCC. There was not statistically relationship between RRM1 and the other three candidate miRNAs. Furthermore, we determined the expression and prognostic value of miR-22-3p and miR-101-3p in HCC. As shown in Figure 4C and D, miR-22-3p was significantly down-regulated in HCC with statistical significance, and its up-regulated expression was positively correlated with the prognosis of HCC patients. Similar to the results of miR-22-3p, miR-101-3p was significantly down-regulated in HCC, and its up-regulated expression positively correlated with the prognosis of HCC patients (Figure 4E and F). Taken together, all these results indicate that miR-22-3p and miR-101-3p may be the most possible upstream regulators of RRM1 in HCC.

Figure 4 Prediction and identification of miR-22-3p and miR-101-3p as potential upstream regulatory miRNAs of RRM1 in HCC. (A) The miRNA-RRM1 prediction network produced by Cytoscape. (B) The correlation between the candidate miRNAs and RRM1 in HCC analyzed by starBase. (C and D) The expression and prognostic value of miR-22-3p HCC assessed by starBase and Kaplan Meier plotter. (E and F) The expression and prognostic value of miR-101-3p HCC assessed by starBase and Kaplan Meier plotter.

Upstream lncRNAs Prediction of miR-22-3p and miR-101-3p

In this section, we first used starbase to predict the upstream lncRNAs of miRNAs. For miR-22-3p, we found 54 possible lncRNAs. Further, we used Cytoscape software to construct the visualization of miR-22-3p regulatory network (Figure S1). For miR-101-3p, we found 28 possible lncRNAs. Similarly, we used Cytoscape software to construct the visualization of miR-101-3p regulatory network (Figure S2). Then, we used GEPIA to evaluate the expression profile of these lncRNAs in HCC, and for the subsequent entry of candidate lncRNAs, we further evaluated their prognostic value in HCC. As shown in Figure 5, for miR-22-3p, compared with the corresponding control, only the expression of MIR4435-2HG was significantly up-regulated among all the 54 lncRNAs candidates (Figure 5A). Subsequently, we further found that the high expression of MIR4435-2HG indicates bad OS (Figure 5C), but not RFS (Figure 5E). Similarly, for miR-101-3p, we found that only SNHG6 was significantly up-regulated in the 28 candidate lncRNAs compared with the control (Figure 5B). Subsequently, it was further found that the high expression of SNHG6 indicated poor OS (Figure 5D) rather than RFS (Figure 5F).

Figure 5 Expression and survival analysis of upstream lncRNAs of miR-22-3p or miR-101-3p in HCC. (A and B) The expression of MIR4435-2HG (A) and SNHG6 (B) in TCGA HCC compared with “TCGA and (or) GTEx normal” data. (C and D) The OS for MIR4435-2HG (C) and SNHG6 (D) in HCC. (E and F) The RFS for MIR4435-2HG (E) and SNHG6 (F) in HCC. ***p value < 0.001.

The competitive endogenous RNA (ceRNA) hypothesis holds that lncRNA can increase mRNA expression by competitively binding shared miRNAs for the purpose of regulating gene expression.¹¹ It follows that lncRNA expression should negatively correlate with miRNA expression. Therefore, the correlation between MIR4435-2HG and miR-22-3p in HCC was evaluated by starBase. Similarly, the correlation between SNHG6 and miR-101-3p in HCC was also included in the analysis. We found that the relationship between the two groups is negatively correlated with statistically significant (Table 1), in line with the ceRNA hypothesis. Combined with all the above analysis, we believed that MIR4435-2HG may be the most potential upstream lncRNA of miR-22-3p/RRM1 axis in HCC, while SNHG6 may be the most potential upstream lncRNA of miR-101-3p/RRM1 axis in HCC. Together, they form the upstream regulatory network of RRM1.

Table 1 Association Between lncRNA and miRNA in HCC Analyzed by starBase

Expression Correlation of RRM1 and Immune Cell Infiltration in HCC

As an important component of the tumor microenvironment, the tumor-infiltrating immune cells promote tumor growth, invasion and metastasis by changing the immune status of tumor cells. We investigated the relationship between RRM1 and immune cell infiltration in HCC, and found that only the copy number of high amplification in CD4+ T cell was significantly overexpressed in HCC (Figure 6A). The subsequent evaluation indicated that the expression of RRM1 significantly positively correlates with all the immune cells analyzed, including B cells, CD8+T cells, CD4+T cells, macrophages, neutrophils and dendritic cells in HCC (Figure 6B–G).

Figure 6 The correlation between immune cell infiltration and RRM1 expression in HCC. (A) Effect of different copy numbers of RRM1 on the infiltration level of various immune cells in HCC. (B–G) Correlation between expression level of RRM1 and infiltration level of B cells (B), CD4 + T cells (C), CD8 + T cells (D), macrophages (E), neutrophils (F) and dendritic cells (G) in HCC. ***p value < 0.001.

Relationship Between RRM1 and Biomarkers of Immune Cells in HCC

We further used GEPIA to explore expression correlation between RRM1 and immune cell biomarkers in HCC. As shown in Table 2, it further showed that RRM1 expression positively correlated with biomarkers of B cell (CD19), CD8+ cell (CD8A and CD8B), CD4+ cell (CD4), Neutrophil (CEACAM8, ITGAM and CCR7), M1 macrophage (NOS2, IRF5, and PTGS2), M2 macrophage (VSIG4 and MS4A4A), or dendritic cell (HLA-DPB1, HLA-DRA, HLA-DPA1, NRP1, CD1C, and ITGAX) in HCC. Together, our results further supported that RRM1 is positively associated with immune cell infiltration in HCC.

Table 2 Association Between RRM1 and Biomarkers of Immune Cells in HCC Analyzed by GEPIA

Correlation of RRM1 with Immune Checkpoints in HCC

The well-known immune checkpoints, PD1/PD1 and CTLA-4, are used by cancer cells to help evade immune surveillance. To further evaluate the possible carcinogenic effect of RRM1 in HCC, we used two different databases to evaluate the relationship between RRM1 and PD1, PD-L1 or CTLA-4. For TIMER, the expression of RRM1 in HCC was significantly positively correlated with PD1, PD-L1 and CTLA-4 (Figure 7A–C). We observed the same positive correlation in GEPIA (Figure 7D–F). These results further indicated the carcinogenic effect of RRM1 in HCC, and tumor immune escape might be involved in the occurrence of this process.

Figure 7 Correlation between RRM1 expression and PD-1, PD-L1 or CTLA-4 expression in HCC. (A–C) The correlation of RRM1 with expression of PD-1 (A), PD-L1 (B), or CTLA-4 (C) in HCC determined by TIMER. (D–F) The correlation of RRM1 with expression of PD-1 (D), PD-L1 (E), or CTLA-4 (F) in HCC determined by GEPIA.

Discussion

HCC is the second leading cause of cancer mortality worldwide with the 5-year survival rate below 20%. Despite some signal regulators related to hepatocarcinogenesis have been clarified, the underlying signaling network remains largely unknown. The more molecular mechanism of liver cancer is elucidated, the more it will help to develop targeted drugs or discover new and ideal prognostic clinical markers. More and more evidences show that RRM1 plays a role as a prognostic marker in a variety of human tumors, including HCC. However, so far, our understanding of RRM1 in HCC is still insufficient and needs further research.

In our study, we first evaluated RRM1 expression profile in 21 human tumors using TCGA data, and then further confirmed the expression of RRM1 using GEPIA. The survival analysis of RRM1 indicated that RRM1 expression was elevated in most cancer types (including HCC), indicating a poor prognosis. Wang et al demonstrated that RRM1 was overexpressed in HCC tissues and was one of the independent prognostic factors in HCC.¹² Kong et al further suggested that RRM1 may exert its importance in HCC progression and prognosis through cell cycle-related pathways.¹³ All of the above reports together with the results of our analysis exhibited the carcinogenic effect of RRM1 in HCC.

In the past few decades, ncRNAs, including small ncRNAs such as miRNAs, lncRNAs and circular RNAs have attracted widespread attention for their extensive involvement in protein synthesis,^14–16 control of transcription,^17–20 and gene expression.^21–24 This has increased the complexity of the mechanisms regulating gene expression. To search for the possible upstream miRNAs of RRM1, seven prediction programs were applied to predict, and finally 13 candidate miRNAs that may bind to RRM1 were obtained. Most of these miRNAs exerted tumor suppressive effects in HCC. The study by Li et al revealed that EIF4G2 promotes hepatocarcinogenesis through activation of the ERK pathway and negatively regulated by the tumor suppressor miR-144-3p.²⁵ miR-144-3p acts on downstream SEPT2 and SGK3 and is regulated by upstream lncRNA NORAD to exert an inhibitory effect on HCC progression.^26,27 miR-338-3p binds to NONO, HN1 and smoothened and is regulated by upstream ST8SIA6-AS1 and hsa_circ_0000092 to inhibit hepatocarcinogenesis and infiltration.^28–30 miR-138-5p acts as an inhibitory activity and predicts negative prognosis by targeting the FOXC1 and PTK2/AKT pathways and is regulated by upstream CircC16orf62 in HCC.^31,32 Several studies have reported that miR-22-3p and miR-101-3p play a role in inhibiting the proliferation and metastasis of HCC.^33–37 Together with correlation analysis, expression analysis and survival analysis, miR-22-3p and miR-101-3p were identified as the most promising upstream miRNAs for RRM1.

According to the ceRNA hypothesis, we used starBase program to predict potential upstream lncRNAs for miR-22-3p and miR-101-3p. For miR-22-3p, we predicted a total of 54 candidate lncRNAs. According to ceRNA theory, lncRNAs should be negatively correlated with miRNAs and highly expressed in HCC. After expression analysis and survival analysis of all 54 candidate lncRNAs, only MIR4435-2HG was screened out. Similarly, for miR-101-3p, we obtained 28 candidate lncRNAs. After expression analysis and survival analysis of all lncRNAs, only SNHG6 was highly expressed in HCC and predicted a poor prognosis. Ultimately, we identified that MIR4435-2HG may be the most promising upstream lncRNA in the miR-22-3p/RRM1 axis in HCC, while SNHG6 may be the most promising upstream lncRNA in the miR-101-3p/RRM1 axis in HCC. Together, they constituted the regulatory network of RRM1.

Tumor-infiltrating immune cells are thought to be associated with tumor suppression, immunosuppression, immunotherapeutic response, and tumor-related proliferation and metastasis during treatment.^38–41 Our study shows that RRM1 expression in HCC is significantly positively correlated with immune cells including B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils and dendritic cells. In addition, RRM1 also correlated significantly and positively with biomarkers of the infiltrating immune cells. Our findings suggest that tumor immune infiltration may be involved in RRM1-mediated hepatocarcinogenesis.

The high level of immune infiltration of tumor tissue drives tumor cells to defensively upregulate various immune checkpoints in order to evade the immune system.⁴² The induction of immune cells into the tumor by immune checkpoint blockade is the key to tumor immunotherapy.⁴³ Thus, the effectiveness of immunotherapy depends primarily on sufficient immune cell infiltration and adequate expression of immune checkpoints in the tumor microenvironment.⁴⁴ Therefore, we evaluated the correlation of RRM1 with immune checkpoints. Our findings showed that high expression of RRM1 in HCC was closely associated with PD1, PD-L1 and CTLA-4, suggesting that targeting RRM1 may improve the efficacy of immunotherapy in HCC.

Conclusion

In conclusion, we clarified that RRM1 is highly expressed in a variety of human cancers (including HCC) and predicts poor prognosis. We further identified an upstream regulatory mechanism of RRM1 in HCC, naming the regulatory network of RRM1 co-regulated by MIR4435-2HG/miR-22-3p and SNHG6/miR-101-3p in hepatocarcinogenesis (Figure 8). Furthermore, our findings suggest that RRM1 may exert its oncogenic effects by increasing tumor immune cell infiltration and expression of immune checkpoints. However, our comprehensive analysis findings need more experiments to verify.

Figure 8 The network mode of RRM1 co-regulated by MIR4435-2HG/miR-22-3p and SNHG6/miR-101-3p in hepatocarcinogenesis.