A pan-cancer analysis of the oncogenic role of ribonucleotide reductase subunit M2 in human tumors

Gene and protein expression analysis of RRM2 in pan-cancer

The UCSC Xena platform was used to retrieve RNA-seq data and clinical information for the TCGA transcriptome (Goldman et al., 2020). Using log2 (TPM), the expression data were transformed to eliminate duplicate and missing values (Vivian et al., 2017). In the Single Gene Analysis module of GEPIA2, the correlation between differential RRM2 mRNA expression analysis and pathological stages analysis of TCGA tumors was determined (Tang et al., 2017). The violin plots were created using log2 (TPM (transcripts per million) + 1) transformed expression data (Huo et al., 2021). The UALCAN portal was used to analyze data from CPTAC dataset. The RRM2 total protein or phosphoprotein level was compared (Chen et al., 2019; Qiu et al., 2021). The RRM2 protein expression was validated using IHC analysis from HPA database, the results were compared to TCGA RRM2 gene expression data (Uhlén et al., 2015).

Clinical value analysis of RRM2 in pan-cancer

To begin, Cox regression analyses were performed on TCGA datasets to determine the effect of RRM2 on pan-cancer patient survival with the R packages “survival”,“survminer” “rms” “Hmisc” and “ggplot2”. We assessed OS, progression-free survival (PFI), and disease-specific survival (DSS). We visualized the univariate Cox regression analysis using the forestplot package and created a nomogram based on the optimal multivariate Cox regression analysis. GEPIA2 was utilized to obtain OS and DFS significant level map data as well as survival plots for RRM2 across tumors (Tang et al., 2019). To determine the hypothesis, the log-rank test was utilized. Then, the R packages “pROC” and “ggplot2” were employed to generate the receiver operating characteristic (ROC) curve using the RRM2 expression profile data from the TCGA databases. The area under the curve (AUC) was computed to quantify the ability to discriminate between the tumor and normal groups (Chen et al., 2021; Xie et al., 2021a; Xie et al., 2021b; Qiu et al., 2022).

Genetic alteration analysis of RRM2 in pan-cancer

We collected data on RRM2 alteration frequency, mutation type, mutated site information, CNA (Copy number alteration), DNA methylation alterations, and 3D structure of the protein throughout TCGA malignancies using the cBioPortal for Cancer Genomics platform (Cerami et al., 2012; Gao et al., 2013).

RRM2-related gene enrichment analysis

For the subsequent analysis of the protein-protein interaction (PPI) network, we employed the Search Tool for the Retrieval of Interacting Genes (STRING) (Szklarczyk et al., 2019). Using Cytoscape (version 3.8.2), the PPI network was constructed, visualized, and analyzed (Otasek et al., 2019). Using GEPIA2, the top 100 RRM2-associated target genes were identified, and Pearson correlation coefficients were calculated between RRM2 and selected genes. Meanwhile, using RNA-seq data from the TCGA database, the top ten correlation results in GEPIA2 correlation analysis were verified using the adjusted Spearman’s rank correlation test.

A heatmap representation of the selected genes’ expression profiles. Then, we used the “clusterProfiler” and “org.Hs.eg.db” and “ggplot2” R packages to merge and filter the top 100 RRM2-correlated genes data in order to conduct GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (Yu et al., 2012; Cui et al., 2020). In addition, GSEA analysis was utilized to explore the possible biological processes involving RRM2 in pan-cancer using “clusterprofiler” R package. The sets of pathway genes were curated using data from the Molecular Signatures Database (MSigDB). Enrichment significance was defined as gene sets with P adjust < 0.05 and FDR < 0.25 (Subramanian et al., 2015).

Association analysis of RRM2 expression with tumor immune microenvironment

To investigate the function of RRM2 in TIME, we determined the relationship between RRM2 and marker genes for immune activating genes, immunosuppressive genes, MMR-related genes, and immune checkpoint genes using the R packages “stat” and “ggplot2” as well as Spearman’s correlation analysis. TIMER is a tool for analyzing the relationship between RRM2 expression and immune infiltrates in all TCGA tumors (Li et al., 2020; Li et al., 2017).

Clinical samples

From 154 patients undergoing hepatic resection at Huashan Hospital, LIHC tumor tissues and paracancerous tissues were obtained. Then, all of the samples were examined using histology. Following that, histology was used to verify all of the samples. The Huashan Hospital Human Ethics Committee approved the experimental protocols, and all patients sign an informed consent form (HIRB: 2021-783).

Cell culture and transfection

Two LIHC cell lines (HepG2, Huh-7) and one immortalized hepatocyte cell line (LO2) were cultured in Dulbecco’s Modified Eagle’s medium (DMEM/F12, BI) supplemented with 10% FBS (Gibco, Carlsbad, CA, USA). These cell lines were obtained from the cell bank of the Chinese Academy of Sciences in Shanghai. All cells were grown in a humidified atmosphere incubator at 37 °C and 5% CO₂. Small interfering RNA negative control 5′-UUCUCCGAACGUGUCACGUTT-3′ (sense), 5′-ACGUGACACGUUCGGAGAATT-3′ (antisense); hs-RRM2 si-1 5′-GGAGUGAUGUCAAGUCCAATT-3′ (sense), 5′-UUGGACUUGACAUCACUCCTT-3′ (antisense); hs-RRM2 si-2 5′-GGAGGAGAGAGUAAGAGAATT-3′(sense), 5′-UUCUCUUACUCUCUCCUCCTT-3′ (antisense) and hs-RRM2 si-3 5′-GCAAGCGAUGGCAUAGUAATT-3′ (sense), 5′-UUACUAUGCCAUCGCUUGCTT-3′(antisense) targeting the coding region of each mRNA, were transfected into HepG2 or Huh-7 cells with Lipo3000 transfection reagent according to the manufacturer’s instructions (HANBIO, China). After 48 h, transfected cells were used for experiments. To generate stable HBx-overexpression cell lines (LO2-HBx), lentiviral vectors pLenti-HBx-HA-Zeo (GeneChem, DaeJeon, South Korea) were constructed and used for corresponding cells by lentivirus-mediated transfection. Before the subsequent experiments, G418 (GeneChem, DaeJeon, South Korea) was used for 14 days to select stable transfections. The efficiency of over-expression and knockdown were validated by qRT-PCR.

RNA extraction and qRT-PCR

Total RNA was extracted from samples using TRIzol reagent (Pufei, Shanghai, China) and reversely transcribed into cDNA using the PrimeScript RT reagent Kit (Takara Bio, China), following by PCR amplification using SYBRs Premix Ex Tap TM II (Takara Bio, Shiga, Japan) in a real-time PCR system (Bio-Rad, Hercules, CA) (Xie et al., 2021a; Xie et al., 2021b). GraphPad Prism was used to conduct statistical analyses. The primers sequences as below: RRM2 forward-CATTGTGAGGTACAGGCGGAAG, reverse-GAAATGGTCTGAGCTGGCAGAAG, HBx forward-CCCGTCTGTGCCTTCTCATC, reverse-CCCAACTCCTCCCAGTCTTTA, ACTB forward-TGGCACCCAGCACAATGAA, reverse-CTAAGTCATAGTCCGCCTAGAAGCA. GAPDH forward-GACATCAAGAAGGTGGTGAAGCAG,reverse-GTCAAAGGTGGAGGAGTGGGT.

Cell viability and cell growth assay

The CCK-8 and the 5-ethynyl-20-deoxyuridine (Edu) assay were used to determine cell viability and growth. According to manufacturer’s protocol, HepG2 and Huh-7 cells were seeded at a density of 2000 cells/well in 96-well plates and treated with negative control and si-RRM2. Cell viability of transfected cells at 48 h was detected using CCK-8 assay (Dojindo, Rockville, MD, USA). Furthermore, the cells were incubated with Edu (Ribobio, Guangzhao, China). Data were obtained from three separate experiments.

Statistical analysis

R language version 3.6.2 was utilized for statistical analysis. This study displayed all data as means ± standard deviation. Student’s t test and GraphPad Prism 6 were applied to determine the significance of the differences.

RRM2 is upregulated in pan-cancer

Using GEPIA2, we analyzed the mRNA expression status of RRM2 across numerous cancer types. RRM2 expression was found to be elevated in the majority of human tumors, according to the findings (Fig. S1). RNA-seq data from the TCGA and GTex databases were used to further confirm the expression of RRM2 in pan-cancer. As depicted in Fig. 1A, RRM2 expression significantly overexpressed among 28 types of cancer. In contrast, it showed little expression in acute myeloid leukemia (LAML). In paired tumor and normal tissues from the TCGA pan-cancer database, RRM2 was highly expressed in bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), thyroid carcinoma (THCA), uterine corpus (UCEC), esophageal carcinoma (ESCA), head and neck cancer (HNSC), kidney renal papillary cell carcinoma (KIRP), rectum adenocarcinoma (READ), stomach adenocarcinoma (STAD), LIHC, lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC) and prostate adenocarcinoma (PRAD) (Fig. 1B).

We investigated the relationship between RRM2 expression and pathological stages of cancer. Remarkably, we discovered clear associations in the majority of instances, including adrenocortical carcinoma (ACC), BRCA, COAD, kidney chromophobe (KICH), skin cutaneous melanoma (SKCM), KIRP, LUAD, kidney renal clear cell carcinoma (KIRC), LUSC, ovarian serous (OV), LIHC and THCA (Fig. 1C, all P < 0.05). In addition to evaluating RRM2 at the transcriptional level, we evaluated it at the protein level using the National Cancer Institute’s CPTAC dataset, which contains large-scale proteome data. The total protein expression of RRM2 was strongly elevated in BRCA, KIRC, COAD, HNSC, LIHC, LUAD, OV, pancreatic adenocarcinoma (PAAD), and UCEC patients (Fig. 2A, P < 0.05, P < 0.001). Meanwhile, we examined IHC information from the HPA database to validate RRM2 protein expression in CPTAC dataset and compared the results to TCGA data on RRM2 gene expression. The outcomes of the analyses of data from the HAP and TCGA database were congruent. Normal liver, thyroid gland, stomach, and colon were negative for RRM2 immunohistochemical staining, whereas tumor tissues were moderately stained (Fig. 2B). Regarding the protein phosphoprotein of RRM2, we analyzed three types of pan-cancer tumors in greater detail using the CPTAC dataset. Figure 2C summarizes RRM2 phosphorylation sites and significant differences. The S20 locus of RRM2 is more highly phosphorylated in BRCA than in normal tissues, but less so in LUAD.

Multifaceted clinical value analysis of RRM2 in Cancers

We downloaded the TCGA mRNA sequencing and clinical data of 33 cancer types from the UCSC Xena platform and used the Cox proportional hazards model to examine the prognostic assessment value of RRM2 in pan-cancer (Table S1). Figure S2 shows the forest plots indicating RRM2 as a high-risk factor for KIRP, KIRC, ACC, LIHC, LUAD, PAAD, brain lower grade glioma (LGG) and PRAD of OS prognosis assessment in pan-cancer, and RRM2 expression was clearly correlated with worse DSS in LIHC, ACC, BLCA, PAAD, PRAD, KIRP, LGG, LUAD, pheochromocytoma and paraganglioma (PCPG) and KIRC. A nomogram integrating RRM2 expression and pathological stage, an independent clinical risk factor, was developed (Fig. 3). The C-index in ACC, KIRC, and KICH indicated the prediction results match the observation results perfectly, as shown in Table S2. Using the Mantel-Cox test and the GEPIA2 database, we analyzed the survival contribution of RRM2 in various cancer types; the survival maps supplemented by OS and disease-free survival (DFS) curves are depicted in Fig. S3. High RRM2 transcription levels were associated with a poor OS prognosis for cancers, including cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC, logrank P = 0.029), ACC (logrank P = 0.00033), KIRC (logrank P = 0.021), KIRP (logrank P = 0.00028), LGG (logrank P = 1.8e−06), LIHC (logrank P = 0.00058), LUAD (logrank P = 1.5e−05), PAAD (logrank P = 0.0058). The analysis of DFS data revealed that patients in ACC (logrank P = 0.00042), KIRC (logrank P = 0.047), KIRP (logrank P = 0.00017), LGG (logrank P = 0.0014), LIHC (logrank P = 0.00032), LUAD (logrank P = 0.019), PAAD (logrank P = 0.00038), PRAD (logrank P = 0.0093), SKCM (logrank P = 0.027), THCA (logrank P = 0.0048). ROC curve analysis indicates that RRM2 in LIHC (AUC = 0.961, 95% CI [0.030–0.983]), KIRC (AUC = 0.929, 95% CI [0.882–0.976]), BRCA (AUC = 0.961, 95% CI [0.944–0.977]), CESC (AUC = 0.988, 95% CI [0.965–1.000]), LUAD (AUC = 0.968, 95% CI [0.952–0.984]), ESCA (AUC = 0.918, 95% CI [0.820–1.000]), glioblastoma multiforme (GBM, AUC = 0.999, 95% CI [0.998–1.000]), COAD (AUC = 0.941, 95% CI [0.907–0.975]), KIRP (AUC = 0.926, 95% CI [0.882–0.970]), UCEC (AUC = 0.970, 95% CI [0.932–1.000]), LGG (AUC = 0.962, 95% CI [0.951–0.972]), LUSC (AUC = 0.995, 95% CI [0.990–0.999]), OV (AUC = 0.996, 95% CI [0.991–1.000]), PAAD (AUC = 0.980, 95% CI [0.967–0.993]), CHOL (AUC = 1.000, 95%CI [1.000–1.000]), uterine carcinosarcoma (UCS, AUC = 1.000, 95% CI [1.000–1.000]) performed well in the diagnostic accuracy (Fig. 4).

Mutation landscape of RRM2 in cancers

According to the results of our analysis, RRM2 mutations include missense, inframe, splice site, truncating, amplification, and deep deletions. Melanoma patients exhibited the highest frequency of RRM2 alterations, with “mutation” being the predominant type, whereas the “amplification” type of CNA and copy number “deep deletion” were the predominant types in LIHC and Mature B-cell Neoplasms, respectively (Fig. 5A). Miss-mutation, amplification, and deep deletion are the most common types of RRM2 mutations (Fig. 5B). Figure 5C depicts RRM2 gene modification types, sites, and case counts. Between amino acids 0 and 389, we detected a total of 92 mutation sites. The RRM2 missense mutation was the most prevalent type alteration, whereas the R358W/Q mutation was identified in one case of GBM and two cases of UCEC. These mutations were dispersed throughout the entire protein sequence and 3D structure (Fig. 5D). Spearman’s correlation was utilized to examine the relationship between promoter DNA methylation and RRM2 expression levels (Fig. 5E). RRM2 expression was significantly inversely correlated with DNA methylation in patients with OV (cg00506866), diffuse large B-cell lymphoma (DLBC, cg09903779, cg00793280, cg05515713), KIRP (cg18623836), LIHC (cg00506866, cg18623836), STAD (cg18623836) and UCS (cg18623836) (−1 < Pearson r ≤  − 4.0). Amplification, gain function, and diploidy were the most common putative copy-number alterations of RRM2 (Fig. 5F). IGHJ2, GDF5-AS1, SMPD4P1, TRAV16, TRAV4, STAG3L4, FGF2, ARHGAP44-AS1, CDRT15P1, and COX10-DT gene alterations were more prevalent in the altered groups (Fig. 5G). These findings may merit additional in-depth investigation.

Enrichment analysis of RRM2-related partners

Using the STRING tool, we obtained experimentally detected RRM2-binding proteins, and the PPI network of proteins was shown in Fig. 6A. The GEPIA2 tool was then used to identify the top 100 genes that correlate with RRM2 expression in combining TCGA pan-cancer expression. As illustrated in Fig. 6B, RRM2 expression was positively associated with that of BUB1 (BUB1 Mitotic Checkpoint Serine/Threonine Kinase, R = 0.72), CCNA2 (Cyclin A2, R = 0.76), CENPI (Centromere Protein I, R = 0.72), CKAP2L (Cytoskeleton Associated Protein 2 Like, R = 0.74), PLK1 (Polo Like Kinase 1, R = 0.74), KIF4A (Kinesin Family Member 4A, R = 0.72), KIF1 (Kinesin Family Member 1, R = 0.75), MKI67 (Marker of Proliferation Ki-67, R = 0.75), NCAPH (Non-SMC Condensin I Complex Subunit H, R = 0.72) genes (all P < 0.001). The majority of tumors exhibited a positive relationship between RRM2 and the nine genes listed above, as indicated by the results of a heat map (Fig. 6C). As shown in Fig. 6D, GO analysis determined that RRM2 was primarily involved in cell cycle checkpoint, ATPase activity, and condensed chromosome. Several pathways, including “Hepatitis B”, “p53 signaling pathway”, and “cell cycle” were identified as the KEGG pathways most significantly enriched, suggesting that RRM2 played a crucial part in the development and progression of cancers, particularly HBV-related LIHC (Fig. 6E). Meanwhile, GSEA was used to identify functional RRM2 enrichment (Fig. S4). High expression of RRM2 is implicated in numerous carcinogenic and immune-related processes, including PLK1 pathways, Forkhead box M1 (FOXM1) pathways, P53 signaling pathways, immunoregulatory interactions between a lymphoid and a nonlymphoid cell, and DNA repair and methylation.

Correlation of RRM2 expression with tumor immune microenvironment

In CESC, ESCA, GBM, LUSC, and thymoma (THYM), RRM2 had a negative correlation with marker genes of immune activating genes, such as VSTR, TNFSF14, TNFRSF14, and STING1 (Fig. 7A). RRM2 was associated with immunosuppressive genes, including PDCD1LG2, IDO1, HAVCR2 and CD274, in COAD, HNSC, DLBC, BLCA, THCA, OV, KIRP, LGG, LIHC, SKCM, and KIRC (Fig. 7C). Additionally, RRM2 expression was strongly correlated with over 30 immune checkpoint markers, including TNFRSF9, TNFSF9, CD70, LGALS9, PDCD1, CD276, HAVCR2, CD28, CD48, CTLA4, ICOS, LAG3 and LAIR, in BRCA, KIRC, LGG, LIHC, THCA, and THYM (Fig. 7B). RRM2 expression was positively correlated with the level of MMR-related genes in the majority of tumors, particularly LIHC, HNSC, GBM, LUSC, PRAD, and UCEC (Fig. 7D).

As for immune infiltration analysis, RRM2 expression level was linked with the estimated infiltration value of Tregs for PRAD, BRCA, LIHC, and THCA (Fig. 8A). In BLCA, COAD, ESCA, KIRP, PAAD, and STAD, RRM2 expression correlated positively with M1-like macrophages and negatively with M2-like macrophages (Fig. 8B). RRM2 expression was significantly correlated with the immune-associated cells T cells in four cancer types, cytotoxic cells in four cancer types, B cells in two cancer types, CD8+T cells in three cancer types, Treg cells in four cancer types, dendritic cells in one cancer type, macrophages in three cancer types, NK cells in six cancer types, and neutrophils in three cancer types. In particular, we discovered a significant correlation between the RRM2 expression level and seven types of immune-associated cells in THYM (Fig. S5).

The validation of RRM2 expression and its oncogenicity

IHC was used to confirm the expression of the RRM2 protein in LIHC patient cohort from Huashan Hospital. We verified the upregulated protein expression of RRM2 in tumor tissue (n = 154) compared to adjacent normal tissue by IHC staining (Figs. 9A and 9B). To demonstrate the oncogene function of RRM2, we conducted experiments with liver cancer cells. We successfully suppressed RRM2 expression in HepG2 and Huh-7 cells using three RRM2 siRNA (Figs. 9C and 9D). Based on the efficiency of knockout, RRM2 siRNA-1 was chosen for subsequent experiments. CCK-8 (Fig. 9E) and Edu assays (Figs. 9F–9H) showed that RRM2 knockdown significantly inhibited cell proliferation in HepG2 and Huh-7, confirming our conclusion. In addition, HBx plasmid was transfected into LO2 cells to explore the regulatory relationship of HBx on RRM2 (Fig. 9I). We found that RRM2 mRNA expression were significantly upregulated in LO2-HBx cells (Fig. 9J).

Validating the associations of RRM2 with immunological makers of Tregs in LIHC patient cohort

Tumors tended to have a high abundance of Tregs cell infiltration, which was associated with protumor immunity. IHC was used to examine the expression of transforming growth factor β1 (TGF-β1) and activator of transcription (STAT)-5B protein from patients tissues. We demonstrated that both TGF-β1 and STAT5B expression was higher in the RRM2-high group than that in the RRM2-low group (Figs. 10A and 10B). TGF-β1 was significantly more strongly positive in LIHC tissues than in adjacent paracancerous tissues. STAT5B expression followed a similar pattern (Fig. 10C). A significant correlation between STAT5B and RRM2 protein levels was observed in LIHC and paracanerous tissues, as well as mixed tissues from LIHC patients (Fig. 10D). A positive correlation between TGF-β1 and RRM2 was also observed in mixed tissues from LIHC patients (Fig. 10E). All these findings suggested that RRM2 may play a significant role in immune infiltration in LIHC, at least in part.