Identification of m7G Methylation-Related miRNA Signature Associated with Survival and Immune Microenvironment Regulation in Uterine Corpus Endometrial Carcinoma

Background N7-methylguanosine (m7G) has been implicated in the development of cancer. The role of m7G-related miRNAs in the survival prediction of UCEC patients has not been investigated. Current research was the first to construct an m7G-related miRNA model to accurately predict the survival of patients with uterine corpus endometrial carcinoma (UCEC) and to explore immune cell infiltration and immune activity in the tumor microenvironment. Methods RNA-seq data and clinical information of UCEC patients were derived from The Cancer Genome Atlas (TCGA) database. Using the TargetScan online database, we predicted miRNAs linked to the m7G-related genes and identified miRNAs which were significantly associated with the survival in UCEC patients and constructed a risk scoring model. The TCGA-UCEC cases were scored according to the risk model, and the high- and low-risk groups were divided by the median risk value. Gene enrichment analysis and immune cell infiltration and immune function analysis were performed using “clusterProfiler” and “GSVA” packages in R. Results The survival prediction model consisted of 9 miRNAs, namely, hsa-miR-1301, hsa-miR-940, hsa-miR-592, hsa-miR-3170, hsa-miR-876, hsa-miR-215, hsa-miR-934, hsa-miR-3920, and hsa-miR-216b. Survival of UCEC patients in the high-risk group was worse than that in the low-risk group (p < 0.001). The receiver operating characteristic (ROC) curve showed that the model had good predictive performance, and the area under the curve was 0.800, 0.690, and 0.705 for 1-, 3-, and 5-year survival predictions, respectively. There were differences in the degree of immune cell infiltration and immune activity between the low-risk and high-risk groups. The expression levels of the identified differentially expressed genes correlated with the susceptibility to multiple anticancer drugs. Conclusions The survival prediction model constructed based on 9 m7G-related miRNAs had good predictive performance.


Introduction
Uterine corpus endometrial carcinoma (UCEC) is a common gynecological malignancy that endangers women's health. According to the 2020 global cancer statistics, there were more than 410,000 new UCEC cases and more than 97,000 new deaths worldwide [1]. In recent years, with the identification of an increasing number of markers and signaling pathways, rising novel approaches to the treatment of UCEC have emerged. However, the prognosis prediction of UCEC patients still needs to be explored. N 7 -methylguanosine, or m7G, is a ubiquitous posttranscriptional RNA modification [2]. A methyl group is added to riboguanosine at the N7 position during transcription initiation [3]. A positively charged modification of RNA results from the modification of m7G at the 5′ cap catalyzed by methyltransferase cotranscriptionally [4]. m7G has been shown to be associated with a variety of biological processes and is involved in the development of various diseases, including cancers [5,6].
In this study, we constructed a model of m7G-related microRNAs (miRNAs) for predicting survival in UCEC patients. Moreover, UCEC patients were divided into highand low-risk groups based on risk score results, and differentially expressed genes were analyzed to explore the correlation with immune cell infiltration, immune activity, and anticancer drug sensitivity. The flow chart of our study is demonstrated in Figure 1.

Materials and Methods
2.1. Data Acquisition. RNA-seq expression data and clinical data of UCEC were derived from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). A list of twenty-eight m7G-related genes was obtained from a previous literature [7] (Supplementary Table S1).

Construction of the m7G-Related miRNA Prognosis
Prediction Model in UCEC. The obtained proteins encoded by 28 m7G-related genes play multiple roles in m7G biological processes, which in turn affect cancer progression [8][9][10][11]. It has been demonstrated that miRNAs regulate gene expression posttranscriptionally [12]. miRNAs can affect tumor progression by regulating the expression of downstream m7G-related genes [13,14]. Based on this, we designed and constructed an m7G-related miRNA prognosis model. Using the TargetScan online database (https://www .targetscan.org/vert_72/) [15], we predicted miRNAs that may have regulatory relationships with m7G-related genes (Supplementary Table S2). Next, based on the transcriptome data obtained from the TCGA-UCEC cohort, differentially expressed miRNAs (DEmiRNAs) between tumor samples and normal samples were obtained by screening with the "edgeR" package in R [16] according to the following conditions: FDR < 0:05 and jlog 2 FCj ≥ 1. To explore the relationship between the expression of DEmiRNAs and the survival of UCEC patients, we used the "survival" package in R to perform the Cox univariate and multivariate analyses and screened out DEmiRNAs with p value < 0.01. For the purpose of obtaining a generalized linear model and reducing error, 1000 Cox LASSO regression calculations and ten cross-validation cycles were performed with the R package "glmnet" [17]. The obtained DEmiRNAs significantly associated with survival were named SDEmiRNAs, and a risk prognostic model was constructed. Finally, to assess the sensitivity and specificity of the survival risk assessment model, we plotted ROC curves using the "timeROC" package in R for predicting survival at 1-, 3-, and 5-year survival [18]. Risk score = 〠Coef miRNA × log 2 miRNA expression + 1 ð Þ : We divided the high-and low-risk groups in UCEC patients by the median risk score.
2.4. Screening of Differentially Expressed Genes. A risk score was calculated for each TCGA-UCEC case according to the constructed survival risk model, and the median risk score was used to distinguish high risk from low risk. Following this, we identified differentially expressed genes (DEGs) between the high-risk and low-risk groups by looking at the transcriptome data in the TCGA-UCEC cohort. Genes that were differentially expressed between the high-risk and low-risk groups were found by screening with the "edgeR" package in R under the following conditions: FDR < 0:05 and jlog 2 FCj ≥ 1. We named the obtained genes RDEGs. Next, we collected the immune score data of TCGA-UCEC from the ESTIMATE database (https://bioinformatics .mdanderson.org/estimate/disease.html), differentiated the high score group and the low score group based on median immune scores, and designated those differentially expressed genes that met the conditions of FDR < 0:05 and jlog 2 FCj ≥ 1 as IDEGs. To obtain codifferentially expressed genes (coDEGs), we took the intersection of RDEGs and IDEGs.
2.5. Gene Enrichment Analysis. We performed functional and pathway enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) on coDEGs using the "clusterProfiler" package in R [19].
2.6. Immune Cell Infiltration and Immune Function Analysis. We used the "GSVA" package in R to quantify immune cell infiltration and immune function in the TCGA-UCEC cohort using the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm [20]. Subsequently, we analyzed the correlation between different immune cell infiltrations and different immune functions. We also explored the association between risk scores and immune infiltration, as well as the correlation between coDEGs and immune infiltration.
2.7. Drug Sensitivity Analysis. We used the GSCA online tool (http://bioinfo.life.hust.edu.cn/GSCA/#/drug) to explore the relationship between the expression of coDEGs and the sensitivity of antitumor drugs, including two major modules: The Cancer Therapeutics Response Portal (CTRP) drug sensitivity analysis and the Genomics of Drug Sensitivity in Cancer (GDSC) drug sensitivity analysis [21].

Statistical
Analysis. For statistical studies, R software (version: 4.1.2) was used. To compare the data between two groups, Student's t-test was utilized. One-way ANOVA was employed to analyze multiple groups followed by Tukey's post hoc test. All tests were two-sided and statistical significance was defined as p < 0:05. The "ggplot2" package in R was used for plotting.

Construction of a Prognostic
Model from m7G-Related miRNAs. We obtained 33 specimens of normal tissue and 546 specimens of tumor tissue from the TCGA-UCEC cohort. The differentially expressed miRNA analysis results of the tumor group compared to the normal group showed 91 upregulated DEmiRNAs and 61 downregulated DEmiR-NAs (Supplementary Table S3). In the heatmap (top 20 DEmiRNAs), the differences in the expression of DEmiRNAs between the normal group and the tumor group are shown ( Figure 2). Next, we identified 9 SDEmiRNAs that were        The median risk score was used to differentiate between the high-and low-risk groups (Figure 4(a)). It is clear from Figures 4(b) and 4(c) that the high-risk group had a significantly higher mortality rate than the low-risk group and the low-risk group had better survival. The sensitivity and specificity of the survival risk assessment model were evaluated by the ROC curve, and the areas under the curve of the risk score model were 0.800, 0.690, and 0.705 for the prediction of 1-year, 3-year, and 5-year survival probability, respectively ( Figure 4(d)).

Predictive Value of Risk Scoring Models Combined with
Clinical Parameters. We incorporated clinical parameters (age and FIGO stage) and risk score into univariate and multivariate Cox regression analyses models and found that the risk score was an independent factor affecting the survival of UCEC patients (Figures 5(a) and 5(b)). In addition, we drew a nomogram for risk assessment combined with clinical parameters to predict the 1-, 3-, and 5-year survival rates of UCEC patients ( Figure 5(c)), and the calibration curve showed that the nomogram had good predictive performance ( Figure 5(d)).

Association of SDEmiRNAs with Survival in UCEC
Patients. Survival analysis was used to explore the correlation between the expression of SDEmiRNAs and the survival of UCEC patients. In Figure 6, the Kaplan-Meier curves showed that UCEC patients with high expression of hsa-miR-876, hsa-miR-934, hsa-miR-940, and hsa-miR-1301 had a significantly lower survival rate than those with low expression, and the survival rate of UCEC patients with high expression of hsa-miR-3170 was significantly increased compared with those with low expression.

BioMed Research International
IDEGs and obtained 1516 coDEGs by taking the intersection (Figure 7(a)). We performed enrichment analysis of the obtained coDEGs with KEGG, and the results showed that the top three pathways were neuroactive ligand-receptor interaction, cAMP signaling pathway, and oxytocin signal-ing pathway (Figures 7(b) and 7(c)). An overview of a given gene's function is provided by GO, which includes molecular function (MF), cellular components (CC), and biological process (BP). In the GO enrichment analysis, BP module's results were especially enriched in cell-cell adhesion via Correlation between GDSC grug sensitivity and mRNA expression   (Figures 7(d) and 7(e)).

3.5.
Tumor Microenvironment Immune Infiltration. The ssGSEA results are represented as a heatmap, and it is likely that different immune cells have different degrees of infiltration in UCEC tumor samples, while at the same time, immune function activities are also different (Figure 8(a)). Then, we examined the correlations between 16 subtypes of immune cells, as well as 13 immune functions (Figures 8(b) and 8(c)). Interestingly, in Figure 8(d), the proportions of CD8+ T cells, DCs, iDCs, neutrophils, and Th2 cells in the high-risk score group were significantly lower than those in the low-risk score group. On the other hand, the proportion of aDCs in the highrisk score group was significantly higher than that in the lowrisk group. In addition, among various immune functions, only the degree of parainflammation and type I IFN response was significantly different between the high-risk score group and the low-risk score group (Figure 8(e)). Subsequently, we analyzed the correlation between the expression of coDEGs and 16 different subtypes of immune cells, as well as the correlation between the expression of coDEGs and 13 immune functions. The results showed that the expression of most coDEGs was negatively correlated with the degree of immune cell infiltration and immune functions, except for IGLV3-24 and IGLV9-49 (Figure 8(f)).
3.6. Drug Sensitivity Analysis. In order to comprehensively understand the correlation between the expression of coDEGs and the sensitivity of antitumor drugs, we used CTRP and GDSC drug data for mining analysis and the results are shown in Figures 8(g) and 8(h).

Discussion
Methylation occurs in biological systems via enzymes. The process of methylation can control gene expression, RNA processing, and protein function. In epigenetics, it is considered a key process. The N 7 -methylguanosine (m7G) modification of RNA has recently gained significant attention. m7G modifications affect a variety of RNA molecules, including messenger RNA, ribosomal RNA, microRNA, and transfer RNA, that participate in biological and pathological functions [6]. It is now clear that m7G plays a critical role in the development of human diseases, especially cancers, and aberrant m7G levels are linked to tumorigenesis and progression through the regulation of multiple oncogenes and tumor suppressor genes [22,23]. Currently, the molecular mechanisms underlying m7G modification in UCEC are not well understood. Therefore, we explored m7G-related miRNA signatures in UCEC and constructed a risk score model for survival prediction of UCEC patients. In our study, we identified 9 m7G-related miRNAs associated with the survival of UCEC patients by integrated bio-informatics approach, including: hsa-miR-1301, hsa-miR-940, hsa-miR-592, hsa-miR-3170, hsa-miR-876, hsa-miR-215, hsa-miR-934, hsa-miR-3920, and hsa-miR-216b. These SDEmiRNAs were used to construct a risk scoring model for the survival prediction of UCEC patients at 1, 3, and 5 years. Previous studies have shown that these miRNAs are associated with the phenotype of various malignant tumors. According to Wang et al., miR-1301 inhibits tumor cell migration and invasion by regulating the UBE4B-p53 pathway in multiple human cancer cells [24]. In gastric cancer, the miR-1301-3p/KIF23 axis inhibits gastric cancer cell proliferation, migration, and invasion by knocking down circ_ 0067934 [25]. Guo et al. reported that PVT1 inhibits proliferation and promotes apoptosis in human retinal epithelial cells by binding to microRNA-1301-3p and activating KLF7 [26]. In pancreatic cancer, Zhang et al. reported that miR-1301-3p inhibits epithelial-mesenchymal transition through targeting RhoA [27]. Wang et al. revealed that miR-1301-3p directly binds to METTL3 and regulates hepatocellular carcinoma progression. It has been reported that METTL3 associates with both m7G cap-binding complexes CBP80 and eIF4E [28]. Researchers found that miR-940 increases proliferation and metastasis in endometrial cancer by regulating MRVI1 [29]. Through regulating FOXO3, miR-940 promotes malignant progression of breast cancer [30]. In esophageal squamous cell carcinoma cells, miR-940 inhibits cell proliferation and promotes apoptosis, which may affect the outcome after surgery [31]. By absorbing hsamir-940, hsa_circ_0092339 targets C-MYC indirectly and plays an important role in castration-resistant prostate cancer [32]. It is reported that miR-592 enhances medullary thyroid cancer tumorigenesis through cyclin-dependent kinase 8 [33]. MiR-592 suppresses the malignant phenotypes of thyroid cancer through the regulation of lncRNA NEAT1 and downregulation of NOVA1 [34]. A miR-592-mediated activation of mTOR (mammalian target of rapamycin), ERK1/ERK2 signaling, and neuronal differentiation impart group 4 medulloblastoma characteristics [35]. It has been found that lncRNA MEF2C-AS1 inhibits cervical cancer by targeting RSPO1 by suppressing miR-592 [36]. There is no relevant experimental research literature on miR-3170, but one literature incorporates miR-3170 into the prognostic model through bioinformatics analysis and considers it to be a predictor of UCEC [37]. In a glucose-induced tumor microenvironment, the HOXC-AS2/miR-876-5p/HKDC1 axis regulates endometrial cancer progression [38]. By activating the PI3K/AKT signaling pathway, miR-876-5p targets GNG12 and contributes to glioma progression [39]. Through targeting TMED3, miR-876-3p modulates the resistance of gastric cancer cells to cisplatin and their stemlike properties [40]. A study by Liang et al. reported that hsa_circ_0097922 regulates ACTN4 expression via miR-876-3p, thus promoting tamoxifen resistance in breast cancer cells [41]. miR-215 interferes with cell migration and invasion by targeting stearoyl-CoA desaturase in colorectal cancer [42]. Through targeting RB1 and triggering the Wnt/β-catenin pathway, miR-215 promotes nasopharyngeal carcinoma progression [43]. By targeting Sox9, miRNA-215-5p inhibits aggressiveness in breast cancer cells [44]. A

19
BioMed Research International decrease in miR-215 levels is correlated with an increase in KDM1B levels in enzalutamide-resistant prostate cancer cells that promotes AR-dependent AGR2 transcription and regulates the sensitivity to AR-targeted therapy [45]. In breast cancer, miR-934 regulates PTEN and the epithelialmesenchymal transition [46]. An increase in the expression of miR-934 has been shown to serve as an independent prognostic factor in lung cancer and contributes to proliferation, migration, and invasion in non-small-cell lung cancer cells [47]. miR-934 derived from colorectal cancer tumors can induce macrophage polarization toward M2 to promote liver metastasis [48]. At present, there is no literature report that miR-3920 is related to tumors, and further research is needed in the future. miR-216b overexpression enhances the activation of PI3K/AKT by partially regulating PXN in gastric cancer cells [49]. Through binding to the 3′-UTR of HDAC8, miR-216b-5p inhibits proliferation and progression in breast cancer cells [50]. Inhibition of osteosarcoma cell proliferation, migration, and invasion is mediated by miRNA 216b targeting Forkhead box M1 [51].
The tumor microenvironment (TME) is a key factor in the occurrence and development of tumors. In particular, various immune cell infiltration and changes in immune activity in the TME affect the biological properties of tumors. Immunotherapy has been known as a promising antitumor regimen in recent years. To propose an appropriate antitumor regimen, it is necessary to first understand the specific conditions of immune cell infiltration and immune activity in the TME. Our findings showed that in UCEC, the high-risk group had a lower degree of immune cell infiltration than the low-risk group. Adaptive immunity is mediated by CD8-positive T cells, a subgroup of MHC class I-restricted cells. These include cytotoxic T cells, which can destroy cancerous or virally infected cells, and CD8-positive suppressor T cells, which suppress specific forms of immunity. Cancer is primarily targeted by CD8+ cytotoxic T lymphocytes (CTLs) [52]. TME-associated immune-related tolerance and immunosuppression causes CTL dysfunction and exhaustion during cancer progression, resulting in adaptive immune resistance [52]. Innate immune cells known as dendritic cells (DCs) infiltrate tumors and present tumor-derived antigens to naïve T cells. Hence, DCs serve as a major therapeutic target for cancer immunotherapy since they play a critical role in priming antitumor immunity [53]. As the most dominant immune cell, neutrophils also play a complex and critical role in cancer. The peripheral blood counts of neutrophils have been found to be elevated in a number of cancer studies. Researchers have demonstrated that the neutrophil-to-lymphocyte ratio can serve as an independent prognostic indicator for cancer patients [54]. Th2 cells have been shown to mediate both pro-and antitumor effects. Studies showed that Th2 cells are indeed involved in mediating antitumor immunity, contrary to the traditional view that type 2 immunity inhibits these responses [55][56][57]. In addition, we also found differences in parainflammation and type I IFN response between the high-and low-risk groups, which may suggest a potential role in TME.
In addition, we analyzed the correlation between the expression of coDEGs and immune cell infiltration or immune func-tion, and the results showed that the negative correlation was more than the positive correlation. The results of drug sensitivity analysis showed that the expression of coDEGs was significantly correlated with a variety of antitumor drugs, which provided suggestions for the selection of anticancer drugs.
Currently, several studies have reported prognostic prediction models for UCEC. Ni et al. constructed a prognostic model composed of 12 miRNAs for UCEC, which can accurately predict the survival of UCEC [58]. Using redoxrelated genes, Geng et al. developed a clinical outcome prognostic model [59]. Compared with previous studies, our study started with m7G-related miRNAs and included only 9 miRNAs, which had better accuracy in 1-year prognosis prediction. And the downstream genes regulated by miR-NAs included in the model were related to immune cell infiltration in the tumor microenvironment, providing a basis for antitumor drug selection.

Conclusion
In conclusion, we constructed a model of nine m7G-related miRNAs for predicting survival in UCEC patients. The TCGA-UCEC samples were divided into high-and lowrisk groups based on the risk model, and the immune cell infiltration and immune function of each group were explored. Moreover, we performed functional analysis of coDEGs and investigated the relationship between their expression and immune activity, as well as their correlation with anticancer drug susceptibility.

Data Availability
The data included in the present study are available from the corresponding author upon reasonable request.

Disclosure
The results shown here are in part based upon the data generated by the TCGA Research Network: https://www.cancer .gov/tcga.