Identification of a Novel Immune Landscape Signature for Predicting Prognosis and Response of Endometrial Carcinoma to Immunotherapy and Chemotherapy

Introduction

Uterine Corpus Endometrial Carcinoma (UCEC) is the most common gynecological cancer. In 2018, 382,069 new cases and 89,929 deaths were reported worldwide (Bray et al., 2018). Despite the emergence of targeted therapy and immunological therapy, the incidence and mortality of UCEC have shown a consistent increase (Lortet-Tieulent et al., 2018). The overall 5-year survival rate can reach 75–86% (Gottwald et al., 2010), however, the survival time of patients with cancer metastasis or recurrence after treatment may drop below 16 weeks (Chaudhry and Asselin, 2009). Besides, the therapeutic regimens such as immunotherapy and chemotherapy are mainly designed according to the clinical stage of patients and regardless of the patients’ varying responses. Therefore, it is an urgent need of the scientific community to build a new prognostic model to identify patients that are at a high risk and suitable for certain regimens.

Surgery is the most preferred route to treat UCEC, commonly supported by radiotherapy and chemotherapy that are designed according to histopathologic parameters of the patients. Surprisingly, chemotherapy may exert different or even opposite effects on patients with identical pathological grade. Furthermore, there is limited evidence regarding the type of patients who can draw benefit from chemotherapy. To further complicate this, immunotherapy can trigger strong response in patients with DNA polymerase ϵ (POLE) mutation, microsatellite instability and high-tumor mutational burden (TMB), however, the difficulty of assessing these factors makes them unsuitable as prognostic markers.

Recently, the tumor immune microenvironment and infiltration of immune cells have been found to be associated with cancer development, prognosis and therapeutic response (Galon et al., 2013; Jain et al., 2016; Pages et al., 2018; Yu et al., 2018). Immune and stromal cells play critical roles in cancer biology. Immune related genes may regulate the infiltration of immune cells, a process that has close correlation with immunotherapeutic response (Binnewies et al., 2018). Therefore, we hypothesized that immune- related genes may be utilized to predict the prognosis and therapeutic response of UCEC patients. In this study, we identified two immune-related genes, their different expression levels have significant prognostic value, and developed a model for predicting the survival and therapeutic response of UCEC patients.

Materials and Methods

Data Sources and Clustering

Downloaded from the Cancer Genome Atlas (TCGA) database were data about mRNA expression data of 547 UCEC patients and clinical characteristics, including age, tumor grade, histological type and clinical stage from TCGA¹ on Dec 1, 2019. All the mRNA expression data were derived from 552 tumor cases and 23 normal cases. 32 patients without well-annotated clinical information and survival time less than 30 days were excluded. After that, 515 patients were obtained. The tumor purity, infiltration level and stromal content were calculated through the ESTIMATE method (Yoshihara et al., 2013). The single-sample gene-set enrichment analysis (ssGSEA) scores were implemented via invoking the R package“GSEAbase”scores. to obtain the enrichment level of 29 immune signatures in each UCEC tissue by evaluating the mRNA expression level of UCEC patients and perform hierarchical clustering of UCEC using R package “hclust” (Barbie et al., 2009; Hanzelmann et al., 2013).

A total of 15 UCEC specimens and 15 adjacent tissues were obtained from patients at the Wuxi Maternal and Child Health Hospital, the Affiliated Hospital to Nanjing Medical University from 2018 to 2019 and routine written informed consent was obtained from all patients. These tissues were used to validate the mRNA and protein expression of KLRC1 and CCL13 in an external set.

Differentially Expressed Genes and Immune-Related Genes

To identity the differentially expressed genes (DEGs) among all of the three groups, we first compared the DEGs between Immunity_L and Immunity_H, Immunity_L and Immunity_M, and Immunity_M and Immunity_H. After that, we imported the immune gene set from Immport database². Then, the overlapping genes were obtained by Venn analysis.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Enrichment Analyses

We performed functional enrichment analyses to investigate the potential mechanisms of different hierarchical clustering based on 29 immune signatures. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were utilized to reveal the enriched biological process, cellular component, molecular function and signaling pathway. Terms with a false discovery rate (FDR) < 0.05 were listed using R package “ClusterProfiler”.

Establishment of the Immune-Related Prognostic Signature

We divided all the cases into a training set and a testing set with the ratio of 1:1. We used the training set to identify the prognostic immune-related genes and to establish the IRPS. The testing set and entire set were used to validate its prognostic capability. First, a univariate Cox regression analysis was used to identify prognosis-related genes in the training set. The inclusion criterion was set at P < 0.05 and least absolute shrinkage and selection operator (LASSO) regression was utilized to minimize the overfitting. We then utilized multivariate Cox model to verify the correlation and developed an immune risk score model using the coefficients of multivariate Cox analysis. The risk score for patients in training set, testing set and total set was calculated using the following equation:

Patients were then divided into high-risk and low-risk groups according the risk score.

Validation of the IRPS

The receiver operating characteristic (ROC) curve was plotted to validate the prognostic value of IRPS. The area under the curve (AUC) was calculated using R package “survivalROC”. The survival analyses were conducted using Kaplan-Meier survival curves and “survival” R package. We also used the decision curve analysis (DCA) curve to obtain the predictive power of the IRPS and other clinical characteristics.

Construction and Validation of a Predictive Nomogram

To fully expand the predictive power of a prognostic model, a nomogram was constructed based on the clinical characteristics of UCEC, including age, stage, grade and histological type. Validation of the nomogram was evaluated using calibration plot.

Gene Set Enrichment Analysis

To identify potential biological mechanism related IRPS, we performed GSEA and GO analysis. KEGG terms with FDR ≤ 0.05 were considered enriched. Based on IRPS, patients were divided into different groups, the different expression genes with a fold change (FC) > 1 and an adjusted P-value < 0.05 were identified using R package “limma”. The GO analysis was then performed using the “clusterProfiler” R package.

Estimate of Tumor-Infiltrating Immune Cells

We used the CIBERSORT tool to quantify 22 types of immunocyte fractions based on TCGA RNA-sequencing data. P < 0.05 was set as the threshold. P < 0.05 was set as the threshold.

TIMER Database Analysis

TIMER is a comprehensive resource for systematical evaluations of the clinical impact of different immune cells on diverse cancer types³. We analyzed the expressions of KLRC1 and CCL13 in UCEC and evaluated their correlation with the infiltration of immune cells. Besides, correlations of KLRC1 and CCL13 expression with markers of several immune cells were also statistically evaluated using Spearman’s correlation and represented via scatterplots.

TISIDB Database Analysis

The TISIDB online platform was used to analyze the correlation of KLRC1 and CCL13 expression with 28 immune infiltrating cells⁴.

Mutation Analysis

We downloaded the mutation data of UCEC patients from the TCGA database⁵ and utilized the maftools to analyze the mutation data. The tumor mutational burden (TMB) score was calculated using following formula:

IPS Analysis

IPS can be generated in an unbiased manner using machine learning based on four major gene categories that determine immunogenicity. The IPS was calculated with z-scores of representative genes associated with immunogenicity. The IPSs of patients were obtained from the Cancer Immunome Atlas (TCIA)⁶.

Immunotherapy Response Prediction

The response to immunotherapy was predicted using an online tool called Immune Cell Abundance Identifier (ImmuCellAI) (Miao et al., 2020), which can estimate the abundance of 24 immune cells from gene expression datasets, including RNA-Seq and microarray data, and predict the patient’s response to an existing immune checkpoint blockade therapy.

Verification of Gene Correlation in GEPIA

To further verify the correlation of KLRC1 and CCL13 expression with immune cells markers, the Gene Expression Profiling Interactive Analysis (GEPIA)⁷ database was employed. Statistical analysis was performed using Spearman’s correlation.

Chemotherapy Response and Candidate Small Molecule Drugs Prediction

The response of chemotherapy in UCEC patients was determined using a public database called Genomics of Drug Sensitivity in Cancer (GDSC⁸). The half-maximal inhibitory concentration (IC50) was estimated which represented the drug response. The potential small molecule drugs for UCEC were predicted using Connectivity map (CMap)⁹. This database comprises of genome-wide transcriptional expression data from small molecule drugs, and can discover the connections between drugs, genes and diseases through the variation of gene-expression profiles. These small molecule drugs were predicted based on 382 DEGs between high-risk and low-risk group with | log2fold change (FC) | > 1 and FDR < 0.05. The 3D structures of the three most significant drugs were obtained from Pubchem¹⁰.

Quantitative Real-Time RT-PCR

Total RNA from 15 UCEC samples and 15 adjacent tissues was extracted using TRIzol reagent (Invitrogen) and the total RNA integrity were checked by RNA 6000 Nano kit. Before reverse transcription to cDNA, 4 × gDNA wiper Mix (Vazyme R323-01), DEPC and total RNA (1 μg) were resuspended and reacted at 42°C for 2 min to remove the residual genomic DNA from total RNA. PrimeScript^®, RT reagent kit was used to synthesize the complementary RNA. The SYBR^®, Premix Ex Taq^TM Kit (TaKaRa DRR041) was utilized to perform real-time quantification. The relative expression levels of target genes were normalized by GAPDH and estimated using the 2^{–△ △ Ct} method. The PCR primers are listed in Supplementary Table 2.

Immunohistochemical Staining

The protein expression levels of CCL13 and KLRC1 were estimated via immunohistochemical (IHC) staining. Briefly, the tissues slides were deparaffinized, rehydrated and treated with 3% H₂O₂ for 15 min to eliminate endogenous peroxidase. Then, antigen retrieval was performed by heating the slides in sodium citrate buffer for 3 min. Next, the slides were incubated with rabbit anti-CCL13 or anti-KLRC1 primary antibodies (Affinity, Biosciences, 1:200) at 4°C overnight. The slides were washed and incubated with HRP-conjugated donkey anti-rabbit secondary antibodies (Abcam) for 15 min. The staining was visualized using DAB solution and samples were counterstained with hematoxylin.

Immunostaining of CCL13 and KLRC1 were analyzed by two pathologists who were blinded to the same information. The staining intensity score was defined on a scale of 0 to 3 in which 0 means no staining, 1 means mild staining, 2 means medium staining and 3 means intense staining. The percentage score of stained cells were also calculated on a scale of 1 to 4 in which 1 represents (0–25%), 2 = (26–50%), 3 = (51–75%) and 4 = (76–100%). In order to obtain the final score, the intensity score and percentage score were multiplied to reach the final score ranging from 0 to 12.

Statistical Analysis

We adopted the R project (version 3.6.2; R Foundation) for all analysis¹¹. The following R packages was adopted in this study (“pheatmap”, “rms” “ggplot2”, “forest plot”, “limma”, “glmnet”, “preprocessCore”, “GSVA”, “survminer”, “survival ROC”, “beeswarm”, “ggstatsplot”). Two-side statistical analyses were performed and samples with P-value < 0.05 were considered statistically significant.

Results

Construction of UCEC Subgrouping

The total workflow is as shown in the following figure (Figure 1). With the help of the ssGSEA scores of 29 immune signatures and R package “hclust”, we divided the patients into three clusters according to immune infiltration: Immunity High (Immunity_H), Immunity Medium (Immunity_M), and Immunity Low (Immunity_L). The three distinct clusters, Immunity_H, Immunity_M, and Immunity_L, showed different immune activities. The hierarchical clustering map was shown in Supplementary Figure 1. We found that the patients in the Immunity_H group had higher ESTIMATE Score, Immune Score and Stromal Score and lower Tumor Purity (Figures 2, 3A-C) than other groups. Besides, the expression levels of most HLA genes were significantly higher in Immunity_H group than that in Immunity_L group (Figure 3D). The type of immune cells was different among three groups (Figure 3E). We also compared the expression of several immune regulators, including PD-1, PD-L1, TIM-3, LAG-3, and CTLA4. The expression levels of these immune regulators in Immunity_H group were all higher than those in Immunity_L group (Figures 3F–J). We then conducted Kaplan-Meier survival analysis which highlighted that patients in three groups had distinct clinical outcomes (P = 0.027, Figure 3K).

Differentially Expressed Genes and Immune-Related Genes

To identity the differentially expressed genes (DEGs) among all of the three groups, we first compared the DEGs between Immunity_L and Immunity_H, Immunity_L and Immunity_M, and Immunity_M and Immunity_H and identified 2314, 411 and 1378 DEGs (Supplementary Figures 2A–C). Then, according to the gene set from Immport database (see text footnote 2), we obtained 1811 immune genes (Supplementary Material B). Then, the overlapping genes were obtained by Venn analysis. Finally, 89 overlapping genes were found differentially expressed in all four subgroups, which suggested to play a crucial role in immune status of UCEC (Supplementary Material C). Thus, the 89 overlapped DEGs were selected as key immune related DEGs for further analysis (Supplementary Figure 2D).

Identification of Potential Biological Function-Related Genes

Go and KEGG analyses were performed which highlighted 89 biological function-related key genes in UCEC. We found that biological functions like such as T cell activation, leukocyte cell-cell adhesion, positive regulation of leukocyte activation, etc. were associated with the identified 89 genes. Furthermore, these genes participated in the KEGG pathways including Cytokine-cytokine receptor interaction, Natural killer cell mediated cytotoxicity, and Viral protein interaction with cytokine and cytokine receptor etc. (Supplementary Figures 3A–D).

Development and Validation of the IRPS

To construct the IRPS based on 89 identified overlapping genes, the univariate Cox regression analysis was utilized to identify prognosis-related genes (Supplementary Table 1). 3 overlapping genes were identified. After that, we used LASSO Cox analysis to decrease overfitting (Supplementary Figures 4A,B). After analysis, three genes were all reserved, including KLRC1, CCL13 and LTA. All these genes were associated with the overall survival of UCEC patients (Supplementary Figure 5). We then performed multivariate Cox proportional hazards regression analysis to build the IRPS (Table 1). Two hub genes were reserved. The mRNA and protein expression of these two gens were presented in Figure 4. The mRNA expression of CCL13 and KLRC1 in tumor tissues were significantly lower than that in the normal tissues (Figures 4A,B). Similarly, the protein expression of CCL13 and KLRC1 was consistent with their mRNA expression (Figures 4C,D). The risk score was obtained according to the corresponding coefficients and the expression levels of hub genes. Risk score was calculated using the following formula:

A cutoff Risk score value of 1.40 was selected based on the median value of the risk score in the training set and used to divide the patients into low-risk (Risk score ≤ 1.4028) and high-risk (Risk score > 1.4028) groups.

In the training set, we found the IRPS can distinguish risk score, survival status and expression of 2 hub genes as displayed in Figures 5A–C. Kaplan-Meier survival analysis showed statistical difference between two groups (Figure 5D), and the areas under the ROC curves (AUC) were 0.599, 0.649, and 0.661 for 1-, 3-, and 5-year survival, respectively (Figure 5E). The testing set cohort and the entire cohort were used to validate the prognostic power of the IRPS model. The distribution of risk score, survival status and expression of two hub genes in the testing and entire sets are presented in Figures 5F–H,K–M. Patients in high-risk group showed worse prognosis than low-risk group in both testing and entire sets (Figures 5I,N). ROC analysis revealed the prognostic accuracy of the IRPS in testing and entire sets (Figures 5J,O).

Furthermore, we analyzed the prognostic power of the IRPS with different clinical features in the entire set. We first represented the data in a heatmap to obtain the general distribution of risk score, hub genes and other clinical features (Supplementary Figure 5A) and then utilized subgroup Kaplan-Meier analysis to evaluate the prognostic value of IRPS in some specific conditions (Supplementary Figures 6B–G). We found the IRPS reached satisfactory prognostic discrimination in patients with age ≤ 60 (Supplementary Figure 6B), grade G1&G2 (Supplementary Figure 6C), grade G3&G4 (Supplementary Figure 6D), histological type endometrial (Supplementary Figure 6E), stage I&II (Supplementary Figure 6F) and stage III&IV (Supplementary Figure 6G).

Construction and Validation of a Prediction Nomogram

We conducted the univariate and multivariate Cox regression analyses to determine whether the IRPS was an independent prognostic indicator for UCEC. According to univariate Cox regression analysis, the hazard ratio (HR) of risk score and 95% confidence interval (CI) were 2.717 (1.407–5.248), 1.783 (1.013–3.140), and 2.191 (1.425–3.370) in training, testing, and entire sets, respectively (Table 2). When turns to multivariate Cox regression analysis, the HR and 95% CI were 2.224 (1.137–4.350) and 1.949 (1.273–2.983) in training and entire set, respectively. However, in testing sets, the HR of risk score and 95% CI were 1.653 (0.956–2.859) (Table 2). According to the univariate and multivariate Cox regression analyses, the age, stage, histological type, grade and risk score are significant prognostic factors and should be involved in the construction of the prognostic models.

To expand the prognostic power of the IRPS and other clinical characteristics, we constructed a nomogram that integrated age, clinical stage, grade, histological type, and risk score. Each parameter was assigned with a score and their total score was calculated. To validate the performance of the nomogram (Figure 6A), 1, 3, and 5-year calibration curves were constructed (Figures 6B–D), which revealed a close association between the predicted and actual curves. We further compared the AUC of IRPS and other clinical characteristics for 1-, 3-, and 5-year survival (Figures 7A–C) and found that the results were not suitable for clinical usage. However, when these factors were combined, the AUC reached 0.736, 0.746 and 0.796 for 1-, 3-, and 5-year survival, respectively (Figures 7D–F), suggesting the combination of IRPS and other clinical characteristics was highly reliable. This methodology was further confirmed by the decision curve analysis (DCA) (Figures 7G,H).

Potential Biological Pathways Associated With IRPS

According to the GSEA analysis. We identified that pathways, such as axon guidance, basal cell carcinoma, glycosaminoglycan biosynthesis chondroitin sulfate, were enriched in the high-risk group, whereas autoimmune thyroid disease, B cell receptor signaling pathway, and chemokine signaling pathway were enriched in the low-risk group (Supplementary Figures 7A,B). Besides, we also identified several immune-related GO terms such as T cell activation, regulation of leukocyte activation, regulation of lymphocyte activation, regulation of T cell activation and leukocyte cell-cell adhesion (Supplementary Figures 6C,D).

Correlation Between IRPS and Immune Cell Infiltration

We used CIBERSORT to obtain the proportion of the 22 immune cells (Figure 8A) and found that the proportions of several types of immune cells, including plasma cells, CD8+ T cells, CD4+ memory T cell, follicular helper T cells and M1 macrophages, were higher in the low-risk group, but those of immune cells like, CD4+ memory T cells, M0 macrophages and mast cells were lower in the high-risk group. Besides, we also investigated the correlation between IRPS and different types of immune cells. The IRPS showed positive correlation with memory B cells, activated dendritic cells, M0 macrophages, mast cells, CD4+ memory T cells, and negative correlation with M1 macrophages, NK cells, CD4+ memory T cells, CD8+T cells and follicular helper T cell (Figures 8B,C).

Correlation Analysis Between 2 Hub Genes and Immune Infiltration Level

To investigate the relationship between the two hub genes, tumor purity and the immune filtrating cells, we used the TIMER database to obtain their relationship. We first analyzed the correlations between CCL13 expression, tumor purity and immune infiltration level of 6 immune cells. The results showed that CCL13 expression level had negative correlation with tumor purity. Besides, CCL13 expression level had significant positive correlations with infiltrating levels of B cell, CD8+ T cell, CD4+ T cell, Macrophage, Neutrophil and Dendritic cell (Supplementary Figure 8A). Similarly, the KLRC1 expression level had negative correlation with tumor purity and positive correlation with infiltrating levels of B cell, CD8+ T cell, CD4+ T cell, Macrophage, Neutrophil and Dendritic cell (Supplementary Figure 8B). The above results were also validated using the TISIDB dataset (Supplementary Figures 9A–C). These results suggest that CCL13 and KLRC1 play a specific role in immune infiltration in UCEC.

We further revealed the relationship between the two hub genes and several immune markers, including CD8+ T cells, T cells (general), B cells, monocytes, TAMs, M1 and M2 macrophages, neutrophils, NK cells, DCs and several functional T cells. After adjustment by purity, the results showed that CCL13 expression level was associated with most immune marker sets of immune cells and different T cells except for the M1 Macrophage and Dendritic cell (Table 3). In addition, the KLRC expression level was associated with nearly all of the immune cells. These results were also verified using the GEPIA database (Table 4). In general, according to the above results, CCL13 and KLRC1 play a remarkable role in immune regulation in UCEC.

In addition, we analyzed the relationships of the mutants of these 2 hub genes with immune infiltrates in UCEC. Compared with the immune infiltration levels in samples with wild type signatures, diverse forms of mutation in two hub genes could inhibit the immune infiltration levels of several immune cells, including CD8+ T cell, macrophages and dendritic cells (Supplementary Figures 8C,D).

Correlation Between IRPS and Mutation Profile

The relationship between IRPS and mutation profile was evaluated in UCEC patients using somatic mutation data. The top 10 mutated genes in high-risk and low-risk group are shown in the Figures 9A,B. And the most frequently mutated genes in high-risk and low-risk group are presented in Figure 9C. The results revealed that somatic mutation was more frequently observed in the low-risk group. And the TMB scores in low-risk group were significantly higher than that in high-risk group (P < 0.05, Figure 9D). Further results revealed that TMB score had negative correlation with IROS (P = 4.015e-09, Figure 9E).

Correlation Between IRPS and Two Therapeutic Regimens

We also analyzed the expression of four immune checkpoint molecules in high-risk and low-risk groups. The results revealed that IRPS was negatively corelated with the listed four immune checkpoint molecules (Figure 10A). Besides, we performed IPS analysis to acquire immunogenicity. The results showed that four molecules displayed higher scores in the low-risk group (Figure 10B). Besides, according to the ImmuCellAI, patients in the low-risk group showed higher immunotherapy response rate compared with patients in the high-risk group (Figures 10C,D), which implied that patients in the low-risk group would benefit from immunotherapy.

Chemotherapy is the most common way to treat UCEC cancer, in this research, we used GDSC database to predict the likelihood of response to several common chemotherapy drugs. We estimated the IC50 of each sample and observed a significant difference of IC50 between high-risk and low-risk groups among eight chemo drugs. Patients in the low-risk group were more sensitive to commonly administered chemodrugs (P = 1.467e-05 for cisplatin, P = 4.412e-06 for gemcitabine, P = 0.039 for paclitaxel, P = 0.002 for bleomycin, P = 1.458e-06 for vinblastine, P = 0.048 for vinorelbine, P = 4.620e-05 for vorinostat, P = 0.005 for methotrexate) (Figure 11). In contrast, the chemotherapeutic response of Docetaxel and Doxorubicin was not significantly different between both groups.

Potential Small Molecular Drugs for UCEC

In order to explore new therapeutic regimens for UCEC, the Cmap database was employed (Lamb et al., 2006). We found eight associated small molecule drugs that are listed in the Table 5. Among these small molecule drugs, the 3D chemical structures of three most significant small molecule drugs were obtained from PubChem (Figure 12).

Discussion

UCEC is the most common tumor affecting female reproductive system, with a 5-year survival rate of 16% in patients with distant metastasis (Siegel et al., 2019). To date, the therapeutic regimens, such as immunological therapy and chemotherapy, are mainly designed according to the clinical stages of the tumor. Due to physiological differences, not all the patients can benefit from the current therapeutic regimens (Brooks et al., 2019). To overcome this challenge, in this research, we developed a model for predicting the survival and therapeutic response of UCEC patients using two immune-related genes.

We first estimated the relative levels of 24 immune cells based on the training set online data and used hierarchical clustering analysis to profile the infiltration of immune cells. The results revealed that the infiltration of immune cells varied much among UCEC patients. We found different tumor purities, immune scores, stromal scores, fractions of different immune cells, expression of several HLA genes and expression of five immune checkpoint molecules (CTLA4, CD274, HAVCR2, LAG3 and TIGHT) among UCEC subtypes. These results strongly suggest that tumor immune microenvironment has different landscapes in UECE patients. Emerging evidence demonstrates that tumor immune microenvironment is closely associated with the prognosis of several cancers¹³ (Liu et al., 2020). In this research, we found the overall survival of three UCEC subtypes differed significantly. The hierarchical clustering analysis is capable of distinguishing patients with different prognosis. However, this method is complicated and produces irrelevant information, making it less clinically applicable.

To overcome the shortcomings, we filtered out two hub genes (CCL13 and KLRC1) closely related to the prognosis of UCEC patients. The mRNA and protein expression levels of both genes were verified via qPCR and IHC. The survival analyses confirmed the ability of IRPS in distinguishing patients with different prognosis. In order to enhance the predictive power of the IRPS, we added other clinical characteristics and built a nomogram model. According to the ROC and DCA curves, this nomogram model exhibited remarkable ability in predicting the prognosis of patients.

We further investigated the potential biological function of IRPS. The GSEA analysis revealed that several immune-related pathways were significantly enriched in the low-risk group. In contrast, the same pathways in the high-risk group were scattered. It is widely acknowledged that innate and adaptive immune cells play a major role in regulating the cancer growth. Increasing evidence shows that some immune cells (like Neutrophils, Macrophage M2, T regulatory cell) can stimulate, while some (like Macrophage M1, CD8+ T cell and Th1 CD4+ T cell) can inhibit cancer growth (Coussens and Werb, 2002; Disis, 2010). In the present study, low-risk and high-risk groups differed in the proportion of immune cells in UCEC, including plasma cells, CD8+ T cells, CD4+ memory T cells, CD4+ memory T cells, follicular helper T cells, M0 macrophages, M1 macrophages and activated mast cells. To be specific, M1 Macrophages, CD4+ and CD8+ T cells and plasma cells were activated in the low-risk group, suggesting that they can inhibit cancer growth and improve the prognosis of UCEC patients.

As listed in the methods section, the IRPS was established using the expression profiles of CCL13 and KLRC1. CCL13 is a gene located on chromosome 17q11.2 that encodes monocyte chemoattractant protein 4 (MCP-4), a Cys-Cys (CC) type cytokine characterized by two adjacent cysteines. In the immunoregulatory and inflammatory processes, CCL13 demonstrates chemotaxis to monocytes, lymphocytes, basophils and eosinophils, but not neutrophils, and plays a role in the accumulation of leukocytes during inflammation. Increasing evidence has confirmed that chemokines and their receptors can facilitate the entry of specific immune cells into tumors, thus enhancing anti-tumor response and improving patient prognosis (Rusakiewicz et al., 2013; Jacquelot et al., 2018).

KLRC1 (Killer Cell Lectin Like Receptor C1), also known as NKG2A, is a protein-coding gene associated with Natural killer (NK) cells. NK cells can mediate the lysis of certain tumor cells and virus-infected cells, and specific humoral and cell-mediated immunity. The protein encoded by KLRC1 belongs to the killer cell lectin-like receptor family, also called NKG2 family, which is a group of transmembrane proteins preferentially expressed in NK cells. KLRC1 can form a complex with KLRD1/CD94 and participate in the recognition of the MHC class I HLA-E molecules in NK cells. Researcher has also proved that crystal structure of CD94-NKG2A in complex with HLA-E bound to a peptide derived from the leader sequence of HLA-G. A previous study found KLRC1 expression changed with CD8+ T cell infiltration in 34 types of human cancers (Chen et al., 2019).

It is well known that tumors can escape the immune system via several mechanisms, including expanding T regulatory cells, inducing the production of certain inhibitory cytokines, altering the function of antigen presenting cells (APCs) (Disis, 2010). In this research, we found that the expression of CCL13 and KLRC1 had a positive correlation with the activation of several types of immune cells. Mutations in these two genes can inhibit the infiltration of some immune cells, especially in CD8+ T cells. Thus, IRPS based on both genes can distinguish cellular immunoactivation and immunosuppression.

We also investigated whether IRPS can provide valuable information about the host response to immunotherapy and chemotherapy. Immune checkpoint molecules are traditional biomarkers for evaluating the therapeutic benefit of immunotherapy. In this research, we found that the expression levels of four immune checkpoint molecules (PD-1, PD-L1, PD-L2 and CTLA4) were significantly low in the high-risk group, suggesting that the patients in the high-risk group might not benefit from immunotherapy based on immune checkpoint inhibitors. Apart from immune checkpoint molecules, tumor mutational burden (TMB) has emerged as a promising predictive biomarker for immunotherapy based on immune checkpoint inhibitors in several tumor types (Chen et al., 2019). High TMB and high neoantigen load have positive correlation with sensitivity to immunotherapy (Hollern et al., 2019). Multiple studies have proved that TMB may be a surrogate for overall neoantigen load²⁰ (Rizvi et al., 2015; Rooney et al., 2015). During the cancer onset, somatic cells mutate and express neoantigens (Gubin et al., 2015). These neoantigens can sometimes induce T-cell-dependent immune responses by activating CD8+ T cells that can recognize those neoantigens and initiate tumor cell lysis (Chen and Flies, 2013).

In this research, when it turns to the evidence regarding immune checkpoint molecules, According to TMB and ImmuCellAI tool, patients in the low-risk group may benefit more from immunotherapy with immune checkpoint inhibitors, but those in the high-risk group may not. Similarly, the patients in the low-risk group were sensitive to chemodrugs such as cisplatin, gemcitabine, paclitaxel, bleomycin, vinblastine, vinorelbine, vorinostat and methotrexate. However, patients in the high-risk group were resistant to these drugs, which may explain their poor prognosis.

Fortunately, we found that several small molecule drugs, such as carbenoxolone, emetine, lovastatin and MG-262, could provide potential benefits for patients in the high-risk group. There is limited research regarding the effects of these drugs on tumor. Carbenoxolone is widely used as an antiulcer drug, but with unknown effect on tumor. Carbenoxolone can also act as the inhibitor of Pannexin 1 (Panx-1) and suppress the migration and invasion of testicular cancer cells to counter cancer progression and metastasis (Penuela et al., 2012; Furlow et al., 2015; Jankowski et al., 2018; Liu et al., 2019). Emetine, a potent anti-protozoal and emetic drug, recent evidence has verified its anti-malarial, anti-bacterial and anti-amoebic effects (Matthews et al., 2013; Hudson et al., 2016; Khandelwal et al., 2017; Yang et al., 2018). Over the past decades, emetine has been reported to have anti-tumor effects on leukemia, ovarian cancer, bladder cancer and lung cancer by inhibiting tumor growth by regulating multiple mechanisms such as apoptosis and autophagy (Moller and Wink, 2007; Moller et al., 2007; Kim et al., 2015; Sun et al., 2015). Lovastatin, an HMG-CoA reductase inhibitor, can decrease cholesterol biosynthesis and is an ideal medicine for treating coronary heart disease. In 2004, lovastatin was found to be a useful adjuvant drug for breast cancer (Shibata et al., 2004). Besides, lovastatin can reduce cancer-related deaths. MG-262, also known as Z-Leu-Leu-Leu-B(OH)2, is a proteasome inhibitor that can reversibly and selectively inhibit chymotryptic activity of the proteasome. As we know, proteasome inhibition has emerged as a novel approach to treat cancer. In some studies, MG-262 has exhibited obvious inhibitory effect on the growth of malignant cells.

The above drugs are untraditional anti-tumor drugs and there is limited evidence of their effects on tumors especially UCEC. However, for patients in the high-risk group, all drugs with potential benefits should be tried. For patients who may not benefit from traditional drugs, adjuvant agents should be tried.

Nevertheless, there are still some limitations in this research. First, this study only includes the immune-related genes and did not take other biomarkers into consideration. Another important issue, there are lots of immune cells varies significantly among individuals, it is hard to distinguish if the gene expression level mainly depends on the varieties of immune cells. The variable number of NK cells might be a confounder for KLRC1 expression. We used Timer2.0 to explore the correlation of KLRC1 and NK cells, the result confirmed that the expression of KLRC1 is positively correlated with NK cells. However, according to the results of CIBERSORT analysis, the proportion of NK cells showed no significant difference in high- and low- risk group. Besides, by analyzing the correlation between the RiskScore and different types of immune cells, we found the RiskScore was negatively correlated with NK cells. Thus, we speculated that the expression of KLRC1 might not mainly depend on the number of NK cells varieties among individuals. It was an independent prognosis related factor. Additionally, this research is based on the online data, and large-sample clinical studies are still needed to validate the predictive value of our IRPS model.

In summary, our study identified two immune-related genes, CCL13 and KLRC1 in the development of UCEC. The IRPS of both genes can predict the prognosis and immune status of UCEC patients and evaluate their therapeutic response to immunotherapy and chemotherapy.

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 Clinical Research Ethics Committee, Wuxi Maternal and Child Health Hospital, The Affiliated Hospital to Nanjing Medical University. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

YZ conceived the study. JL, YW, and JM participated in the design, analysis, and draft of the study. JL and YW plotted all figures in this manuscript. JM helped in data analysis. All authors approved the final version of this manuscript and agreed to be accountable for all aspects of the work.

Funding

This study was financially supported by the Wuxi Science and Technology Bureau Project (No. CSE31N1720), the Jiangsu Province Youth Medical Talent Fund Project (No. QNRC2016166), and the Jiangsu Provincial Six Talent Peaks Project (No. YY-124).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Cao Yongke, Ph.D., from Nanjing KEY Translation Studio for editing the English text of a draft of this manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.671736/full#supplementary-material

Abbreviations

UCEC, Uterine Corpus Endometrial Carcinoma; LASSO, least absolute shrinkage and selection operator; IRPS, immune-related prognostic signature; GEPIA, Gene Expression Profiling Interactive Analysis; GDSC, Genomics of Drug Sensitivity in Cancer; TMB, tumor mutational burden; TIME, tumor immune microenvironment; TCGA, the Cancer Genome Atlas; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; FDR, false discovery rate; ROC, receiver operating characteristic; AUC, area under the curve; DCA, decision curve analysis; GSEA, gene set enrichment analysis; TIICs, tumor-infiltrating immune cells; KLRC1, Killer Cell Lectin Like Receptor C1.

Footnotes

1. ^ http://cancergenome.nih.gov/
2. ^ https://www.immport.org
3. ^ https://cistrome.shinyapps.io/timer
4. ^ http://cis.hku.hk/TISIDB/index.php
5. ^ https://portal.gdc.cancer.gov
6. ^ https://tcia.at/home
7. ^ http://gepia.cancer-pku.cn/index.html
8. ^ https://www.cancerrxgene.org
9. ^ https://www.broadinstitute.org/connectivity-map-cmap
10. ^ https://pubchem.ncbi.nlm.nih.gov/
11. ^ http://www.r-project.org/

References