Development of an Oxidative Phosphorylation-Related and Immune Microenvironment Prognostic Signature in Uterine Corpus Endometrial Carcinoma

Background: As the fourth most common malignant tumors in women, uterine corpus endometrial carcinoma (UCEC) requires novel and reliable biomarkers for prognosis prediction to improve the overall survival. Oxidative phosphorylation (OXPHOS) is found to be strongly correlated with the progression of tumor. Here, we aimed to construct an OXPHOS-related and immune microenvironment prognostic signature to stratify UCEC patients for optimization of treatment strategies. Method: Prognosis-associated OXPHOS-related differentially expressed genes were identified by multivariable Cox regression from TCGA–UCEC cohort. Based on the candidate genes, an OXPHOS-related prognostic signature was constructed by the train set data and verified by the entire set. When integrated with relevant clinical characteristics, a nomogram was also created for clinical application. Through comparison of tumor microenvironment between different risk groups, the underlying mechanism of the model and the inner correlation between immune microenvironment and energy metabolism were further investigated. Results: An OXPHOS-related signature containing ATP5IF1, COX6B1, FOXP3, and NDUFB11 was constructed and had better predictive ability compared with other recently published signatures in UCEC. Patients with lower risk score showed higher immune cell infiltration, higher ESTIMATE score (p = 2.808E−18), lower tumor purity (p = 2.808E−18), higher immunophenoscores (IPSs) (p < 0.05), lower expression of mismatch repair (MMR) proteins (p < 0.05), higher microsatellite instability (MSI), lower expression of markers of N6-methyladenosine (m6A) mRNA methylation regulators, higher tumor mutation burden (TMB) (p = 1.278E−9), and more sensitivity to immune checkpoint blockade (ICB) (p < 0.001) and chemotherapy drugs, thus, possessing improved prognosis. Conclusion: An OXPHOS-related and immune microenvironment prognostic signature classifying EC patients into different risk subsets was constructed in our study, which could be used to predict the prognosis of patients and help to select a specific subset of patients who might benefit from immunotherapy and chemotherapy, thus, improving the overall survival rate of UCEC. These findings may contribute to the discovery of novel and robust biomarkers or target therapy in UCEC and give new insights into the molecular mechanism of tumorigenesis and progression of UCEC.


INTRODUCTION
Uterine corpus endometrial carcinoma (UCEC) is the fourth most common malignancy among women in the United States, with a trend of increasing morbidity and mortality worldwide (Cancer Genome Atlas Research et al., 2013). According to the GLOBOCAN, there would be an estimated 382,069 (2.1% of the total cancer cases) new cases, with 89,929 (0.9% of total cancer deaths) deaths in 2018 (Bray et al., 2018). Due to the advance in techniques for early diagnosis, the patients diagnosed at an early stage account for approximately 75% of UCEC, and most of the patients in stage I can be nearly cured by surgery, with a 5years overall survival rate of 47-69%. In comparison, advanced-stage patients (stage III or IV) have a poor prognosis due to high risk for recurrence and limited therapeutic strategies, with a 5-years overall survival rate of 15-17% (Lee et al., 2017). Therefore, attention should be paid to discovering novel and reliable biomarkers for prognosis prediction, which can also work as a sensitive classifier for a specific subset of UCEC patients who will benefit from immunotherapy and have improved survival.
Recently, increasing attention has been paid into the exploration of the relationship between cancer cell metabolic plasticity and migration and metastasis (Mosier et al., 2021). Cancer cells with a more invasive and distal metastasis phenotype are found to be strongly correlated with the upregulated expression of peroxisome proliferator-associated receptor gamma and coactivator 1-alpha, which are hallmarks of active mitochondrial biogenesis and oxidative phosphorylation (OXPHOS) (LeBleu et al., 2014). OXPHOS is supposed to be a potential biomarker of tumor progression (Commander et al., 2020). In addition, emerging evidence shows that OXPHOS inhibition can be a valuable target in cancer treatment (Ashton et al., 2018;Molina et al., 2018;Xu et al., 2020a;Cardenas et al., 2020). OXPHOS provides the required energy for cancer cells to strive, and cancer stem cells with primary or acquired resistance against chemotherapy or tyrosine kinase inhibitors are characterized with markedly enhanced OXPHOS dependency (Sica et al., 2020). Metformin and thiazolidinediones can inhibit the mitochondrial electron transport chain, which provides energy for ovarian cancer growth. Thus, interfering with the process of OXPHOS may be considered as a new target for cancer therapy (Nayak et al., 2018). With the development of molecular medicine and next-generation sequencing technology, OXPHOS has been employed in the cancer risk prediction model in lung adenocarcinoma (Xu et al., 2020b). OXPHOS has also attached great importance to UCEC. As a regulator participating in oxidative phosphorylation and glycolysis, PKM2 can function as a biomarker for malignant and premalignant endometrial lesions, and the presence of PKM2 high tumor cells in UCEC tissue also indicates a poor prognosis (Lai et al., 2019). To the best of our knowledge, OXPHOS has not been included in the UCEC prognostic model in previous studies, which deserves to be explored for subgrouping of the patients with distinct prognosis for personal treatment and identifying potential novel therapeutic targets.
Here, we identified prognosis-associated OXPHOS-related differentially expressed genes (DEGs) from TCGA-UCEC cohort (Tomczak et al., 2015). Based on the candidate genes, we constructed an OXPHOS-related prognostic signature model using the train dataset and verified it using the entire dataset. To further investigate the underlying mechanism of the model and its correlation with immune microenvironment and energy metabolism, then between high-and low-risk groups, we compared tumor-infiltrating immune cells (TIICs), ESTIMATE score, tumor purity, immunophenoscores (IPSs), the expression of mismatch repair (MMR) proteins, mRNA expression-based stemness index (mRNAsi), microsatellite instability expression (MSI), the expression of markers of N 6 -methyladenosine (m 6 A) mRNA methylation regulators, tumor mutation burden (TMB) and the response to immune checkpoint blockade (ICB) and chemotherapy agents, which may provide evidence for clinical application of the signature for immunotherapy and chemotherapy option ( Figure 1). The findings in this study not only provided a useful predictive model for UCEC but also offered evidence for mutual impacts of OXPHOS, immune microenvironment, and tumor malignant biology.

Data Collection
The gene expression profiles were collected from TCGA-UCEC database (https://tcga-data.nci.nih.gov/tcga/), which contained 552 tumor and 23 normal tissues. Meanwhile, the corresponding clinical and pathological information was also obtained from TCGA database, including patient age, tumor grade, tumor stage, and histological type of each sample. After removing the samples with unknown survival time, <30 days, and no survival status, we finally integrated the transcriptome and complete clinical data, enrolling a total of 511 patients in this study (Table 1).

Oxidative Phosphorylation-Related Survival-Related Differentially Expressed Genes and Functional Annotation
We used the R package limma to define the differentially expressed genes among 289 OXPHOS-related genes acquired from the National Center for Biotechnology Information-GENE and the Metabolic Atlas database (Xu et al., 2020b). The genes with p < 0.05 and |log2 (fold change)| >1 were regarded as DEGs. We used the survival package to show the association between OXPHOS-related DEGs and survival time. We defined the DEGs with p less than 0.05 as survival-related DEGs (sDEGs). To explore the potential biological functions of the OXPHOS-related sDEGs, we implemented the ClusterProfiler package of R software for Gene Ontology (GO) term enrichment analysis (Ashburner et al., 2000), which included three categories: biological process (BP), molecular functions (MF), cellular components (CC), as well as the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (Kanehisa and Goto, 2000).

Construction of an Oxidative Phosphorylation-Related Prognostic Signature
We randomly assigned 256 patients as the train set and took the entire set for validation. First, we utilized the univariate Cox regression in the train set to identify candidate sDEGs with a significant level of p < 0.05. Then, we performed the least absolute shrinkage and selection operator (LASSO) Cox regression analysis by the glmnet package (Qu et al., 2018). Finally, we employed multivariate Cox proportional hazard regression analysis to construct an OXPHOS-related prognostic signature. The formula for calculating risk score was as follows: risk score n i 1 the expression level of n* the regression coefficient calculated by multivariate Cox proportional hazard regression.
We treated the median of risk score in the train set as the cutoff value to divide the validation data into the high-and low-risk groups and plotted the Kaplan-Meier curves, as well as the ROC curves by survminer, survivalROC package (Heagerty et al., 2000). To confirm the independence of the signature, a conjoint univariate and multivariate Cox analyses were conducted. In

Gene Set Enrichment Analysis and Immunity Analyses
Gene set enrichment analysis (GSEA) was utilized to elucidate the molecular mechanisms of the OXPHOS-related sDEGs. We separated the samples in the entire set into the high-and lowrisk groups based on scores calculated by the signature and then compared the enriched BP between the two groups. The value of p < 0.05 was regarded as the cutoff criterion (Zhu et al., 2016). Furthermore, we employed the ESTIMATE algorithm by the ESTIMATE package to calculate the stromal, immune, and estimate scores of each sample and evaluate the relationship between the OXPHOS-related prognostic signature and immunization (https://bioinformatics.mdanderson.org/publicsoftware/estimate/) (Yoshihara et al., 2013).

Analysis of Tumor-Infiltrating Immune Cells
We employed CIBERSORT (https://cibersort.stanford.edu/) to explore the immune cell infiltration of each sample based on the RNA-seq data (Newman et al., 2015). To further analyze the association between risk score and tumor-infiltrating immune cells (TIICs), single-sample gene-set enrichment analysis (ssGSEA) was also performed to quantify the immune activity in different risk groups by exploring 29 immune-related genes .

Analysis of Immune Status Between Different Risk Groups
The immunogenicity of tumor is determined by the four dominant components containing immunomodulators, immunosuppressive cells, effector cells as well as MHC molecules. Immunophenoscores (IPSs) (including IPS, IPS-CTLA4, IPS-PD1-PD-L1-PD-L2, IPS-PD1-PD-L1-PD-L2-CTLA4 scores) are calculated based on the representative celltype gene expression z-scores to evaluate and compare the potential response to immune checkpoint inhibitor between different risk groups (Charoentong et al., 2017). IPS data of each UCEC patient in the entire set was downloaded from The Cancer Immunome Atlas (TCIA) (https://tcia.at/home). Besides, the immunotherapy response could be predicted with high accuracy by an Immune Cell Abundance Identifier (ImmuCellAI) result-based model, which precisely estimates the abundance of 24 immune cell types including 18 T-cell subsets from gene expression data (http://bioinfo.life.hust.edu. cn/ImmuCellAI) (Wang et al., 2020a;Toh et al., 2021). Microsatellite instability (MSI) (Toh et al., 2021) and the expression of markers of N 6 -methyladenosine (m 6 A) mRNA methylation regulators (Wang et al., 2020a) of each sample were also collected from TCGA, which can predict potential response to immune checkpoint inhibitors (ICIs) and help explain the association between the signature and immunogenicity. Tumor cells sharing more similarities with stem cells represents cancer progression. Malta et al. developed an innovative mRNAsi algorithm to calculate the degree of similarities for UCEC patients in TCGA database (Malta et al., 2018). We also downloaded the mRNAsi information from the research for analyses.

Analysis of Tumor Mutation Burden
The mutation data collected from TCGA were analyzed by R package maftools (Robinson et al., 2017). The tumor mutation burden (TMB) was calculated by the formula: TMB (total mutation/total covered bases) * 10 6 .

Statistical Analysis
All statistical analyses were applied by R version 3.6.1 (Package: limma, survival, ClusterProfiler, glmnet, survminer, survivalROC, rms, ESTIMATE, pRRophetic). Categorical variables were presented by counts and percentages. Continuous variables in normal distribution were analyzed using Student's t-test and presented as mean ± standard deviation, while in abnormal distribution, they were presented as median (range). Multiple groups of continuous variables were analyzed by one-way ANOVA. The hazard ratio and 95% confidence interval were calculated to identify genes associated with overall survival. Unless with special explanation, p < 0.05 was considered statistically significant.

Identification of Oxidative Phosphorylation-Related Differentially Expressed Genes
We conducted the differential expression analysis of 289 OXPHOS-related genes between tumor and normal endometrium (Xu et al., 2020b). Sixty-six OXPHOS-related DEGs were identified, including 54 upregulated genes and 12 downregulated genes (Supplementary Figures S1A,B and Supplementary Table S3). The GO term enrichment analysis presented that they were related to ATP metabolic process (p < 0.0005) and purine ribonucleoside triphosphate metabolic process (p < 0.0005) in the BP group, mitochondrial inner membrane (p < 0.0005), and organelle inner membrane (p < 0.0005) in the CC group, oxidoreductase activity, acting on NAD(P)H and electron transfer activity (p < 0.0005) in the MF group (Supplementary Figures S1C); the KEGG results manifested that they mainly took part in pathways of oxidative phosphorylation (p < 0.05), endocannabinoid signaling (p < 0.05), and HIF-1 signaling (p < 0.05) in cancer (Supplementary Figures S1D).

Establishment of Oxidative Phosphorylation-Related Prognostic Signature
Using univariate Cox regression from the train data, we defined seven OXPHOS-related sDEGs, including ATP5IF1, MRPL12, FOXP3, NDUFA13, ATP5F1E, and NDUFB11 (Supplementary Table S4). ATP5IF1 and FOXP3 were protective elements of UCEC, while the others were dangerous ones. Then we conducted LASSO Cox regression analysis and multivariate Cox regression analysis (Supplementary Figures S2A,B) to construct a reliable OXPHOS-related prognostic signature. Consequently, ATP5IF1, COX6B1, FOXP3, and NDUFB11 were included in the signature as key prognosis-associated ones (Supplementary Figures S2C). Furthermore, the results of the expression analyses of the four genes showed that all of them were highly expressed in tumor compared with normal tissue, but only ATP5IF1 and FOXP3 were highly associated with overall survival (OS) of UCEC patients (Supplementary Figures S3A-F). The qRT-PCR results of our own samples also demonstrated that the expression level of ATP5IF1 (p 0.0311), COX6B1 (p 0.0212), FOXP3 (p 0.1082), and NDUFB11 (p 0.0450) were significantly higher in UCEC tissues than in normal endometrial tissues (Supplementary Figures S3G-J). The formula to calculate the risk score of patients is shown as follows: risk score (−0.011678 * expression value of ATP5IF1) + (0.0014649 * expression value of COX6B1) + (−0.309569 * expression value of FOXP3) + (0.0029791 * expression value of NDUFB11).

Validation of Prognostic Signature in Uterine Corpus Endometrial Carcinoma
We used the validation data to evaluate the performance of the risk model in survival prediction. Using the preceding equation, we calculated the risk score and divided the samples of the train and entire sets separately into high-and low-risk groups by taking the median risk score in the train set as the cutoff line. The distribution of risk score (Figure 2A), survival status of each sample ( Figure 2B), and the key gene expression profiles in each set ( Figure 2C) are exhibited in Figure 2. By comparison with the low-risk group, Kaplan-Meier survival curves revealed that the OS of the high-risk group was significantly worse (log rank test, p 1.35E−4) ( Figure 2D). Besides, a survival ROC curve analysis was conducted, and the area under the curve (AUC) at 1 year in the train set was 0.748 ( Figure 2E), which verified the predictive value of the signature, and the AUC value in the entire set also suggested a similar potential for prediction of survival with the above results ( Figure 2E). The PCA plots indicated that patients in different risk groups tended to distribute differently ( Figure 2F).

Clinical Utility of the Oxidative Phosphorylation-Related Prognostic Signature
The conjoint univariate and multivariate analyses in both sets confirmed the independence of the OXPHOS-related prognostic signature ( Figure 3). Integrated with clinical factors, we performed subgroup survival analysis in stages I and II, stages III and IV, grades 1 and 2, grades 3 and 4, endometrial, mixed, and serous, age ≤6 0, age >60, and we found that a higher risk score was related to poorer prognosis in all subgroups (Supplementary Figures S4A-H). The correlation between risk core and clinical characteristics is shown in Supplementary  Figures S4I-L and a heatmap (Supplementary Figure S5), indicating that patients with lower stage, lower grade, endometrial histological type, and lower age tended to have a lower risk score. Aimed at creating a feasible method for clinical application, we constructed a nomogram cooperating risk core, patient age, histological type, tumor grade, and tumor stage ( Figure 4A). By comparing the AUC values at 1, 3, and 5 years plotted by the signature and clinical factors, we discovered that the gene-based signature was superior to others in the OS prediction. What is more, when the gene-based signature was combined with the clinical characteristics, it had better predictive ability than the signature alone ( Figure 4B). The calibration curves also demonstrated that the results predicted by the nomogram were in high consistency with actual survival, which suggested that the nomogram had a high clinical value ( Figure 4C).

Gene Set Enrichment Analyses
To ascertain involved biological processes, we employed GSEA to analyze the transcript message of patients in different risk groups. In the high-risk group, representative KEGG pathways were axon guidance, glycosaminoglycan biosynthesis ghondroitin sulfate, glycosaminoglycan biosynthesis heparan sulfate, RNA polymerase, and tight junction ( Figure 5A), while in the lowrisk group, representative KEGG pathways were cytokine-cytokine receptor interaction, chemokine signaling pathway, and so on ( Figure 5B). Through comparison, we discovered that in the low-risk group, the representative biological processes were mainly related to immunization. Furthermore, we implemented the ESTIMATE algorithm and found that the higher the risk score, the lower the immune score, the stromal score, and the ESTIMATE score, and the higher the tumor purity ( Figures 5C-F), which further proved that the signature was significantly related to immunization.

Difference of Tumor-Infiltrating Immune Cells in Different Risk Groups
To further explore the relationship between the OXPHOSrelated prognostic signature and immunization, we found the difference of TIICs between different risk groups. M0 macrophages, M2 macrophages, and activated mast cells were significantly elevated in the high-risk group (p < 0.05). In contrast, plasma cells, CD8 + T cells, activated memory CD4 + T cells, and regulatory T cells (Tregs) were evidently elevated in the low-risk group (p < 0.05) ( Figure 6A). We analyzed the correlation of risk score and TIICs; as we can see in Figures  6B,C, M0 macrophages, M2 macrophages, plasma cells, activated memory CD4 + T cells, CD8 + T cells, gamma delta T cells, and regulatory T cells (Tregs) were the seven most relevant types of immune cells with risk core. There was a significant difference in the enrichment scores of diverse immune cell subpopulations and related functions or pathways with ssGSEA between different risk groups, such as TIL and HLA, which showed coherence with the GO and KEGG analyses results (Figures 6D,E). , and IPS-PD1-PD-L1-PD-L2-CTLA4 (p 3.776E−11) in the highrisk group were clearly higher than that in the low-risk group ( Figures 7A-D). The expression of immune modulators including CD27, CTLA4, ICOS, PD-L2, TIGIT, PD-1, LAG3, TIM-3, CD86, PD-L1, CD70, and CD270 (p < 0.001) in the lowrisk group was also remarkably higher than in the high-risk group ( Figure 7E). We further discovered that the genes related to immune checkpoint inhibitor (ICI), including TIGIT, ICOS, CTLA4, CD27, CD58, PD-1, CD86, TIM-3, IDO1, and PD-L1 were negatively related to risk score ( Figures 7F-Q), which suggested that the patients in the low-risk group promised to respond to immunotherapy.
By comparing the reaction to immune checkpoint blockade (ICB), we found that the percentage of patients with response to ICB in the high-risk group (53.18%, 142/267) was lower than that in the low-risk group (75.41%, 184/244) ( Figure 8A). The risk scores of the responders were also lower than those of the nonresponders ( Figure 8B). These results could demonstrate that the OXPHOS-related prognostic signature had the potential to predict the possible effect of ICB in UCEC patients and work as a classifier for therapy selection.
Moreover, we compared the MSI status of different risk groups and found, as shown in Figures 8C-G, the expression level of dominant mismatch repair (MMR) proteins, such as MLH1, MSH2, MSH6, PMS2, and EPCAM, in the low-risk group were lower, and risk score was negatively correlated with the MSI expression ( Figure 8I). The patients in the low-risk group had a lower score of mRNAsi ( Figure 8H), which can partly explain the poor prognosis in the high-risk group. In addition, significantly decreased gene expression of KIAA1429, FMR1, HINRNPC, ZC3H13, YTHDF1, and YTHDF3, which are key factors of m 6 A mRNA methylation could be found in Figure 8J. These findings all provided evidence for the efficiency of immunotherapy for patients in the low-risk group and helped to elucidate the mechanisms of different outcomes. Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 753004 8

Relationship Between Tumor Mutational Burden and Risk Score
To further clarify the mechanism of tumorigenesis and development of UCEC patients, we compared the tumor mutational burden (TMB) between different risk groups based on somatic mutation data. We discovered that the low-risk group has remarkably higher TMB than the high-risk group (p 1.278E−8) ( Figure 9A). The correlation between risk score and TMB was negative (Cor −0.269, p 6.178E−11) ( Figure 9B). We also found that only the mutation of TP53 was found in the high-risk group, whereas the others were found in the low-risk group ( Figure 9C). Next, we divided UCEC patients into H-TMB and L-TMB groups based on the median TMB; as shown in Figure 9D, the OS of the H-TMB group was better than that of the L-TMB group (log-rank test, p < 0.001). Stratified survival analysis was conducted to assess the synergistic effect of these factors, which corroborated that TMB would not influence the effect of the signature in prognosis prediction. Significant survival differences were observed in both TMBstratified subgroups [log-rank test, H-TMB and high risk score (HH) versus H-TMB and low risk score (HL), p < 0.001; L-TMB and high risk score (LH) versus L-TMB and low risk score (LL), p < 0.001; Figure 9E], which indicated that the signature could function as an independent predictive indicator and effectively evaluate the potential response to immunotherapy. The most frequent somatic mutations in the high-risk group were as follows: PTEN > TP53 > PIK3CA > ARID1A > TTN > PIK3R1 > CTNNB1 > KMT2D > CHD4 > CSMD3 ( Figure 9F). While the most frequent somatic mutations in the low-risk group followed the order PTEN > ARID1A > Figure 9G). Other mutation details, including variant classification, mutation type, SNV class, mutation load, and top 10 mutated genes in UCEC between different risk groups, are shown in Supplementary Figure S6.

DISCUSSION
Though the overall survival state of UCEC patients has improved a lot due to the development of early diagnostic techniques (Braun et al., 2016), there are still 15% of patients with UCEC diagnosed at an advanced stage, thus, suffering from tumor recurrence owing to a limited response to unsuitable therapies (Li and Wan, 2020). Therefore, developing valid and reliable biomarkers for survival prediction and treatment selection to avoid overtreatment in patients who will not relapse and propose individual adjuvant treatment to patients who will recur counts (Lee et al., 2017), and efforts have been put into developing a robust prognostic signature based on DEGs for identification of biologic subsets to guide treatment strategies for UCEC patients, such as a cell cycle-related prognostic signature (Liu et al., 2020a), an autophagy-related prognostic signature , and so on. However, there is no consensus on predictive biomarkers for scientifically validated therapy. In recent years, emerging evidence has proved that OXPHOS could function as a new target in cancer therapy (Ashton et al., 2018). However, OXPHOS has not been involved in developing a prognostic signature for UCEC, thus far, and the potential underlying mechanism and possible clinical application of OXPHOS-related signature in cancer are worth exploring.
Based on the 289 OXPHOS-related gene matrix obtained from the research of Zihao Xu, we used differential expression analysis and identified 66 OXPHOS-related DEGs, including 54 upregulated genes and 12 downregulated genes; most of them have been reported to be associated with cancer. For instance, EPAS1 could enhance the effect of paclitaxel on breast cancer cells by inhibiting growth and promoting apoptosis of MCF-7/TAX cells (Song et al., 2020). Furthermore, functional annotation of the OXPHOS-related DEGs also indicated that the candidate genes were significantly related to ATP metabolic process, purine ribonucleoside triphosphate metabolic process, and the essential pathways reported in many cancers (Matsuda et al., 2017;Fallah and Rini, 2019), including retrograde endocannabinoid signaling and HIF-1 signaling pathways. Then, we randomly selected samples from the entire set in half into the train set and constructed the model based on it. After univariate Cox analysis and LASSO Cox regression analysis, finally, a prognostic model composed of four OXPHOS-related genes (ATP5IF1, COX6B1, FOXP3, and NDUFB11) was built. Except for ATP5IF1, other candidate genes included in the prognostic signature have been proven to play a vital role in cancer development. COX6B1, as one of the core proteins involved in the molecular basis of the below-background radiation stress response, might inhibit the proliferation of laryngeal squamous cell carcinoma cells . FOXP3 can activate the Wnt/beta-catenin signaling pathway, thus, inducing epithelial-mesenchymal transition and promoting the progression of non-small cell lung cancer (Yang et al., 2017). In addition, the regulation of the expression of the NDUFB11 gene can regulate the programmed cell death process, so it promises to be a new target for cancer therapy (Panelli et al., 2013).
Furthermore, Kaplan-Meier survival curves in the train and entire set indicated that the patients with high risk score tend to FIGURE 8 | The response to immune checkpoint blockade (ICB) therapy and the microsatellite instability (MSI) status of different risk groups. (A) Comparison of the response to immune checkpoint blockade (ICB) therapy between the low-and the high-risk groups. Patients in the low-risk group showed more sensitive response to ICB therapy than those in the high-risk group. (B-G) Patients in the low-risk group had lower expression of the dominant mismatch repair (MMR) proteins including MLH1, MSH2, MSH6, PMS2, and EPCAM than patients in the high-risk group. (H) The mRNAsi of the low-risk group is lower than that of the high-risk group. (I) The correlation between the MSI expression and the risk score. (J) The expression of the markers of m 6 A mRNA methylation regulators in the high-and low-risk groups. m 6 A, N 6 -methyladenosine. In the high-risk group, N 267; in the low-risk group, N 244.
Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 753004 have a worse prognosis. Besides, a survival ROC curve analysis demonstrated the predictive value of the OXPHOS-related prognostic signature. In addition, by conjoint univariate and multivariate analyses, the independence of the OXPHOSrelated signature for prognosis prediction was verified. By combining the clinical features, we performed subgroup survival analysis in different risk groups and discovered that high risk was related to poor prognosis in all subgroups. Primarily, we observed that the risk score was highly correlated to tumor stage, tumor grade, histological type, and patient age. What is more, multi-ROC curve analyses and the evaluation of the nomogram provided evidence for the potential clinical application of the signature in prognosis prediction. The results of GSEA indicated that the representative pathways in the low-risk group were mainly related to immunization. What is more, the immune score, stromal score, and ESTIMATE score of the low-risk group were higher, and the tumor purity was lower, which implied that patients with a low risk score appeared to be more immunogenic and might preferentially benefit from immunotherapy. To penetrate the immune mechanism of the signature, we compared TIICs between the two different risk groups and analyzed the relationship between the infiltrating immune cells and risk score. We found that plasma cells, CD8 + T cells, activated memory CD4 + T cells, and regulatory T cells (Tregs) were remarkably elevated in the low-risk group, while M0 macrophages, M2 macrophages, and activated mast cells were significantly elevated in the high-risk group. Immune cell profiles in the tumor microenvironment (TME) of UCEC tissues could be utilized to predict the survival of patients with UCEC (Li and Wan, 2020). Under the influence of environmental stimuli, tumor-associated macrophages (TAMs) can differentiate into two types of macrophage with opposite functions, in which M1 macrophages are protective immune cells that can remove pathogens and malignant cells. In contrast, M2 macrophages are dangerous factors, which can promote angiogenesis and the growth of tumor. High M2 macrophage infiltration in most tumor types is correlated with an adverse prognosis (Zhao et al., 2017;Wang et al., 2018), which is consistent with our results. A meta-analysis proved that UCEC patients with a high CD8 + T-cell density had a favorable prognosis (Guo et al., 2020). Increased expression of activated CD4 + memory T cells had been confirmed to be beneficial for prognosis in bladder cancer , and the previous study of our team also provided evidence for the opinion (Liu et al., 2020c). The regulatory T cell (Treg) could suppress the antitumor immune response, thus, becoming a target of immunotherapy (Tanaka and Sakaguchi, 2017;Antomarchi et al., 2019), which is different from our results. However, another study showed that FOXP3 lymphocytic infiltration had no significant impact on the survival state of UCEC patients (Giatromanolaki et al., 2008). Our results suggested that CD8 + T cells, activated memory CD4 + T cells, regulatory T cells (Tregs), M0 macrophages, and M2 macrophages were crucial to the development of UCEC, which had the potential to be targeted for immunotherapies. The immune cell subpopulation enrichment analyses, as well as the ssGSEA-related functions or pathway analyses, showed that the scores of all the items except for aDCs and type I IFN response were higher in the low-risk group, which further proved that UCEC patients with low risk scores appeared to be more immunogenic. Furthermore, patients with low risk scores had significantly higher scores of IPS, IPS-PD1-PD-L1-PD-L2, IPS-CTLA4, and IPS-PD1-PD-L1-PD-L2-CTLA4; thereby, they tend to respond more sensitively to immune checkpoint inhibitors. The Food and Drug Administration (FDA) has approved a series of monoclonal antibodies (mAbs) targeting PD-1, PD-L1, and CTLA-4 in several malignancy therapies (Arora et al., 2018). Besides, several clinical trials demonstrated the clinical value of anti-PD-1 antibody, pembrolizumab, which can inhibit tumor immune escaping in endometrial cancer. It has been approved as an alternative therapy for patients in the advanced stage (Musacchio et al., 2020). Also, the expression of CD27, CTLA4, ICOS, PD-L2, TIGIT, PD-1, LAG3, TIM-3, CD86, PD-L1, CD70, and CD270 were also remarkably higher in the low-risk group, which deserved to be explored as a new target for immunotherapy. TIGIT could enhance T-cell function by inhibiting PVRIG, so TIGIT-PVRIG pathways play a vital role in human cancers (Whelan et al., 2019). LAG-3 has been identified on various tumor cells, including UCEC, as a critical inducer for the FIGURE 10 | Prediction of the response to common chemotherapy drugs between the low-and high-risk groups. Patients in the high-risk group possessed higher estimated IC50 than patients in the low-risk group. In the high-risk group, N 267; in the low-risk group, N 244. malignancy progression (Friedman et al., 2020). There is also increasing evidence for TIM-3 as an emerging target for malignancy immunotherapy (Moore et al., 2019). The difference found in the comparison of the response to immune checkpoint blockade (ICB) therapy between different risk groups showed the same trend with the above findings. We also compared the risk score of the patients showing a different response to ICB, and found that the responders tended to possess a low risk score. These results demonstrated that the OXPHOSrelated prognostic signature had the potential to be the classifier for a specific subset of UCEC patients, who would have a favorable effect on ICB therapy. Compared with UCEC patients in the advanced stage, patients with primary tumor tend to have more somatic loss of mismatch repair (MMR) protein expression; thus, the loss of MMR should be considered as guidance for immunotherapy (Ta et al., 2018). We compared the MSI status of different risk groups and found that the expression of dominant MMR proteins (MLH1, MSH2, MSH6, PMS2, and EPCAM) and mRNAsi in the low-risk group was downregulated, which illustrated that the patients with low risk score were burdened with high-MSI and showed fewer stemness features associated with oncogenic dedifferentiation. Haoya Xu et al. found that cell stemness can predict patient prognosis and constructed an mRNAsi-related prognostic signature in endometrial cancer (Xu et al., 2021). Accumulating evidence revealed that m 6 A mRNA methylation participates in the pathogenesis and progression of cancers through molecular mechanisms, such as inhibiting the antitumor response of CD8 + T cells Sun et al., 2019). Patients with a low risk score had decreased gene expression of the m 6 A mRNA methylation, thus, showing a favorable response to immunotherapy.
To explore the cause of tumorigenesis, we compared the TMB of each group, and we discovered that the low-risk group had remarkably higher TMB, and the risk score was negatively correlated with TMB. It is reported that patients with high TMB had a more sensitive response to immune checkpoint blockade, such as anti-PD-1 agents (Fancello et al., 2019). A multicohort study demonstrated the value of tissue TMB (tTMB) as a predictive biomarker for response to pembrolizumab monotherapy (Marabelle et al., 2020). In our present study, the most frequent somatic mutations in the low-risk group were PTEN, ARID1A, PIK3CA, TTN, PIK3R1, CTCF, MUC16, ZFHX3, KMT2D, and MUC5B. PTEN functions as a tumor inhibitor by decreasing the activity of phosphoinositide 3kinase (PI3K) (Bian et al., 2018), and its expression diminished with the progression of UCEC (Xi et al., 2019). It has been proven that ARID1A alterations indicated improved prognosis of tumor patients for compromising MMR and enhancing the infiltration of lymphocytes (Jiang et al., 2020). UCEC patients with high PIK3CA had shorter survival time because PIK3CA could change the tumor immune microenvironment by altering the fraction of tumorassociated neutrophils (Pan et al., 2019). Eighty-five percent of lung squamous cell carcinoma patients with TTN mutant type had missense variations and favorable overall survival (Cheng et al., 2019). In addition, patients in the high TMB group showed better overall survival than in the low TMB group, and the combined Kaplan-Meier analyses proved that risk score was a reliable predictive indicator independent of TMB. Our results were consistent with previous studies, indicating that UCEC patients with low risk score tend to have improved response to immunotherapy and favorable prognosis.
These findings indicated that patients with low risk scores possessed higher IPS, higher TMB, and responded more sensitively to ICI and chemotherapy, which partly explained the reason for the improved prognosis of the low-risk patients.
However, there are still some deficiencies with our study. First of all, our data were based entirely on TCGA databases and lacked an external database for validation. Second, our molecular results lack the evidence of the experiments in vitro or in vivo. Third, these OXPHOS-related prognostic genes are supposed to be further examined by cell function assay for biological function and pathway analyses.

CONCLUSION
In conclusion, an OXPHOS-related and immune microenvironment prognostic signature classifying EC patients into different risk subsets was constructed in our study. This model could predict the prognosis of patients and help to select a specific subset of patients who might benefit from immunotherapy and chemotherapy, thus, improving the overall survival rate. These results may also contribute to the identification of novel immune biomarkers or target therapy in UCEC and give a new insight into the tumorigenesis and progression of UCEC.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Nanjing Medical University. The patients/ participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
JB, HY, JL, TC, and ZZ conceived the study and participated in the study design, performance, and manuscript writing. MY, SN, and SY conducted the bioinformatics analysis. LC and FS revised the manuscript. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We would like to thank the researchers and study participants for their contributions. Supplementary Figure S6 | Variant classification, mutation type, SNV class, mutation load and top 10 mutated genes in UCEC between low (A) and high risk groups (B). In the high risk group, N 267; in the low risk group, N 244.