Reclassification of endometrial cancer and identification of key genes based on neural-related genes

Endometrial cancer (EC) is the most common gynecologic malignancy, and its incidence has been increasing every year. Nerve signaling is part of the tumor microenvironment and plays an active role in tumor progression and invasion. However, the relationship between the expression of neural-related genes (NRGs) and prognosis in endometrial cancer remains unknown. In this study, we obtained RNA sequencing data of EC from The Cancer Genome Atlas (TCGA). Endometrial cancer was classified into two subtypes based on the expression of neural-associated genes (NRGs), with statistical differences in clinical stage, pathological grading, and prognosis. A prognostic prediction model was established by LASSO-Cox analysis, and the results showed that high expression of NRGs was associated with poor survival prognosis. Further, CHRM2, GRIN1, L1CAM, and SEMA4F were found to be significantly associated with clinical stage, immune infiltration, immune response, and important signaling pathways in endometrial cancer. The reclassification of endometrial cancer based on NRG expression would be beneficial for future clinical practice. The genes CHRM2, GRIN1, L1CAM, and SEMA4F might serve as potential biomarkers of EC prognosis.

Introduction

In 2020, endometrial cancer has been the sixth most frequent cancer in women worldwide, with 417,000 new cases and 97,000 deaths (1, 2). Although the incidence has leveled off in recent years, it still has been increasing at a rate of 1% per year, making it one of the few human cancers with a rising fatality rate (3). The choice of surgery, radiation, hormonal and/or chemotherapy, immunotherapy, and targeted therapy depend on the stage of the disease.

The tumor microenvironment is composed of neuronal cells, tumor cells, fibroblasts, and immune cells. Cancer cells can generate electroactive tissue by connecting with neural synapses, which drives cancer cells to migrate and develop (4, 5). Active crosstalk between nerves and tumor cells was first observed in prostate and gastric cancers (6, 7), but its role in endometrial cancer remains largely unknown. Perineural infiltration is a new metastatic pathway in endometrial cancer. Endometrial cancer cells have been found to migrate along the neuropil in vitro, which is associated with DRG and also a risk factor for perineural infiltration (8). Through cytokinesis, sympathetic nerve endings in the uterus release norepinephrine, adenosine triphosphate (ATP), and other molecules with oxytocic and contractile properties. Uterine parasympathetic fibers primarily release acetylcholine to regulate myometrium activity (9). Furthermore, estrogen and progesterone play important roles in reshaping uterine innervation in response to cyclical changes from puberty to menopause.

Using comprehensive genomic analysis of TCGA, Talhouk et al. classified endometrial cancers into four different subgroups: POLE, microsatellite instability, low copy number, and high copy number (10). Clinically, TCGA molecular typing is practical, useful, and beneficial in predicting the prognosis of patients. Here, we reclassified endometrial carcinoma based on the expression of NRGs. The research on the nerve and endometrial cancer crosstalk can assist the identification of the treatment targets for endometrial cancer. As a result, we found two subgroups related to prognosis by clustering endometrial cancer patients in TCGA data based on NRGs. To find prognosis-associated genes, a prognostic model of neural-associated genes was created. In this study, we investigated the relationship between NRGs and endometrial cancer prognosis, clinical staging, pathological grading, signaling pathways, immune infiltration, and immune response in the hopes of assisting future research.

Materials and methods

Source and processing of data sets

The Cancer Genome Atlas (TCGA) database (https://Portal.gdc.cancer.gov/repository) was used to obtain raw data from 552 endometrial cancer patients (11). In addition, clinical endometrial cancer data was retrieved, including survival time, survival status, age, grade, and stage information. 42 NRGs were identified from a previous comprehensive review (12).

Consensus clustering

The chi-square test and the R language package were used to examine the correlations between clustering and clinical features (13). The “ConsensusClusterPlus” package was used to separate endometrial cancer cases into two subgroups (14). The packages “survival” and “survminer” provided survival analysis between subtypes. The “ggplot2” package was used to identify genetic differences between typings (15), while the “pheatmap” tool was used to create heatmaps. The prcomp function in the statistics package was used to perform principal component analysis (PCA) (16).

Differential expression analysis

The R package Limma (v3.40.2) was used to study mRNA differential expression (17). Adjusted p-values (FDR) were analyzed in TCGA to correct for false-positive results. The screening conditions for differentially expressed mRNA were | log2FC| ≥ 1 and FDR < 0.05. Determining the cutoff value by the median is the most commonly used method for determining the cutoff value. Similarly, the cut-off value is determined by the interquartile range. High GRIN1/L1CAM expression (top 25%) and low GRIN1/L1CAM expression (bottom 25%) were defined. Since the total expression of CHRM2/SEMA4F was relatively low, 50% was used as the cutoff for high and low expression of both genes. The LASSO regression algorithm was used for feature selection, and 10-fold cross validation was used (18).

Enrichment analysis and ssGSEA analysis

GO enrichment analysis (19) and KEGG enrichment analysis (20) were done by using R packages ClusterProfiler (21). The Cox regression analysis was performed to identify prognostic genes significantly associated with overall survival (OS) in patients with endometrial cancer (p<0.01) (22). Survival curves were constructed using the R packages “survival” and “survminer”. ROC curves were made using the R packages “survivalROC” and “timeROC”. We collected some functional pathways and calculated the functional pathway scores according to the ssGSEA algorithm.

Analysis of immune infiltrates

The CIBERSORT (23) and EPIC (24) in the R package “immunedeconv” (https://grst.github.io/immunedeconv) were used to analyze immune infiltrates of different subtypes. It was also visualized using the R package (v4.0.3) ggplot2 and pheatmap (15).

Analyses of immune checkpoint genes

The correlation of neural-related gene expression with 8 commonly used immune checkpoint genes (CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LC2, SIGLEC15 and TIGIT) was analyzed and visualized using the R package (v4.0.3) ggplot2 (15) and pheatmap.

Algorithm for Predicting Immune Responses

Treatment response to immune checkpoint inhibitors can be predicted using the TIDE algorithm (25).

Stemness analysis

The stemless of mRNA was evaluated using the OCLR method (26).

Statistical analysis

All statistical analyses were performed in R software (v4.0.3) and were statistically significant at P<0.05.

Results

Identification of EC subtypes based on neural-related genes

We obtained EC data from the TCGA database to explore the relationship between 42 neural-related genes and endometrial cancer. To identify the subtypes, we used the R package ConsensusCluster Plus and two clusters showed up in the result: cluster I (C1) and cluster II (C2) (Figures 1A-C). C1 and C2 were found to be well split into two subgroups using principal component analysis (Figure 1D). From the retrospective 10-year clinical follow-up study, the overall survival (OS) rates of the two groups were statistically different (P<0.05), with C1 having a greater OS than C2 (Figure 1E). In addition, there was a statistically significant difference between C1 and C2 in terms of clinical staging and pathological grading (P<0.05). (Figure 2) (Table 1).

FIGURE 1

Figure 1 Classification of EC subtypes based on neural-related genes (NRGs). (A) Cumulative distribution function (CDF) curve and delta area curve of consensus clustering. (B) Heatmap of consensus clustering. Rows and columns represent samples, and different colors represent different categories. (C) Heatmap of neural-related gene expression in different subtypes of EC. High expression is represented by red, whereas low expression is represented by blue. (D) Principal component analysis (PCA) is a method of analyzing data. (E) Survival curves based on the Kaplan-Meier method. Different subgroups’ overall survival curves.

FIGURE 2

Figure 2 Clinical characteristics of C1 and C2. Demonstration of the proportion of different clinical features in different subgroups. (A) Ethnicity. (B) Stage. (C) Grade. (D) Primary, Recurrence and Metastasis. (E) Chemotherapy and Hormone Therapy. The * means that the difference in clinical features between the two groups is statistically significant (p<0.05).

TABLE 1

Table 1 Clinical characteristics of C1 and C2.

Differential expression analysis and enrichment analysis of C1 and C2

We used the R package Limma to search differentially expressed genes (DEGs) in C1 and C2, with FC>2 and P<0.05 as screening criteria. When compared to C2, C1 showed 209 up-regulated DEGs and 325 down-regulated DEGs. In C1, 24 DEGs were down-regulated (e.g. SEMA4F), 8 DEGs were up-regulated (e.g. ADRB2), and 10 DEGs were unregulated (e.g. GDNF) compared to C2 (Figures 3A, B, S1). To investigate the activated or suppressed signaling pathways in C1 and C2, KEGG and GO analyses were used (Figures 3C, D). KEGG analysis showed that compared to C2, C1 displayed the activation of tumor suppressor pathways like estrogen signaling pathway, ferroptosis, IL-17 signaling pathway, and amino acid metabolic pathway, with the suppression of cancer-associated signaling pathways such as gastric cancer, basal cell carcinoma, Wnt signaling pathway, and neurotransmission. Also, C1 demonstrated the inhibitition of multisystem diseases such as cardiomyopathy, hepatitis C, cushing’s syndrome, human papillomavirus infection, and pathogenic escherichia coli infection. In addition, C1 inhibited lipolysis in adipocytes and oxytocin signaling pathway associated with endometrial carcinogenesis (Figure 3C). GO analysis suggested that compared to C2, C1 induced activation of the processes including humoral immune response and antibacterial humoral response, with the inhibition of synapse organization, neuron projection guidance, modulation of chemical synaptic transmission, extracellular matrix organization, cell−cell adhesion, which are generally recognized to be key processes that promote cancer growth and metastasis (Figure 3D).

FIGURE 3

Figure 3 Differential expression and enrichment analysis of C1 and C2. (A) Volcano plots of C1 and C2 differentially expressed genes. Blue represents genes with high expression, red represents genes with low expression, and gray represents genes without differential expression. (B) Heat map of differential expression. (C) C1 activates or suppresses the KEGG pathway when compared to C2. (D) C1 activates or suppresses the GO pathway in comparison to C2.

Immune status analysis of C1 and C2

We used the immunedeconv R package to assess immune infiltration of C1 and C2. the CIBERSORT showed that there were statistically significant differences between C1 and C2 in B cell naive (P<0.05), CD8+ T cells (P<0.001), T cell CD4+ memory activated (P<0.01), T cell regulation (Tregs) (P<0.001), myeloid dendritic cell resting(P<0.001), myeloid dendritic cell activated (P<0.001), and neutrophil (P<0.01), indicating that C2 exhibited stronger immunosuppression compared to C1 (Figure 4A). The EPIC further confirmed that C1 and C2 showed stronger immunosuppression in terms of T cell CD4+ (P<0.01) and T cell CD8+ (P<0.001) with statistically significant differences (Figure 4B). We also used the ggplot2 and pheatmap R packages to analyze the ICG of both EC subtypes, and the study showed that CTLA4, HAVCR2, PDCD1, and TIGIT expression were increased in C1 (P<0.001) compared to C2 (Figure 4C), implying that immunotherapy may be more effective. Furthermore, we utilized the TIDE algorithm to predict cancer immune response, and the findings revealed that the C2 group had a higher TIDE score than the C1 group, with a significant difference (P=7e-05), indicating that the C1 group may benefit more from immune checkpoint inhibitor therapy (ICBs) (Figure 4D). The OCLR algorithm revealed that C2 had a greater stemness index than C1, with a statistical difference (P=0.0017), indicating that the C2 group had a higher degree of cancer progression, which could help identify novel targets for anticancer drugs (Figure 4E).

FIGURE 4

Figure 4 Analysis of C1 and C2 immune infiltration, immunological response, and stemness. (A, B) CIBERSORT and EPIC scores reveal a difference in immune cell infiltration between C1 and C2. (C) Immune checkpoint-associated genes are expressed differently in C1 and C2. (D) TIDE scores for C1 and C2 groups were compared using the TIDE algorithm for predicting cancer immune response. (E) The difference in stemness index between C1 and C2 as calculated by the OCLR method. (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001).

Prognostic analysis of neural-related gene expression in EC

Using the LASSO regression technique, we attempted to determine a link between NRGs and EC prognosis. High expression of neurologically linked genes was shown to be associated with a poor prognosis (P<0.05) (Figures 5A-D). The expression of 42 NRGs may be a prognostic biomarker for EC patients: the area under the curve (AUC) for 1 year, 3 years, and 5 years was 0.689, 0.693, and 0.653 respectively (Figure 5E). Individual prognostic analysis revealed that four of the 42 NRGs were statistically different and linked with EC prognosis (P<0.01) (Figure 5F). CHRM2 and GRIN1 were shown to be positively related to EC prognosis, while L1CAM and SEMA4F were found to be adversely related (Figure 5F). These results suggest that the expression of neural-related genes including CHRM2, GRIN1, L1CAM, and SEMA4F may be potential biomarkers of EC prognosis.

FIGURE 5

Figure 5 Analysis of neural-associated gene expression in EC for prognosis. (A) Coefficients of 42 neural-associated genes represented by λ parameter. (B) LASSO COX regression model was used to draw the partial likelihood deviance versus log(λ). (C) The correlation between risk and survival. A curve, scatter plot, and heatmap are all used to represent the data. (D) Curves of overall survival for high and low risk groups. (E) The receiver operating characteristic (ROC) analysis of risk scores. At 1, 3, and 5 years, the area under the curve (AUC) was 0.900, 0.919, and 0.939, respectively. (F) A univariate Cox analysis was used to look for genes linked to endometrial cancer prognosis.

Correlation of CHRM2/GRIN1/L1CAM/SEMA4F expression with clinical characteristics

Based on RNA sequencing (RNA seq) and clinical data from the TCGA database, we divided the expression of CHRM2, GRIN1, L1CAM, and SEMA4F into high and low expression groups, and defined the cutoff points for high expression of GRIN1 and L1CAM (top 25%) and low expression of GRIN1 and L1CAM (bottom 25%). Because the total expression of CHRM2 and SEMA4F was low, the cut-off point for high and low expression of both genes was set at 50%. The results showed that the high expression of CHRM2 was negatively correlated with EC grading and staging (Figure S2A) (Table 2); the high expression of GRIN1 was negatively correlated with EC grading and staging (Figure S2B) (Table 3); the high expression of L1CAM was positively correlated with EC tumor grading and staging (Figure S2C) (Table 4); the high expression of SEMA4F was positively correlated with EC tumor grading and staging (Figure S2D) (Table 5).

TABLE 2

Table 2 Clinical characteristics of CHRM2 high expression group and CHRM2 low expression group.

TABLE 3

Table 3 Clinical characteristics of GRIN1 high expression group and GRIN1 low expression group.

TABLE 4

Table 4 Clinical characteristics of L1CAM high expression group and L1CAM low expression group.

TABLE 5

Table 5 Clinical characteristics of SEMA4F high expression group and SEMA4F low expression group.

The biological significance of CHRM2/GRIN1/L1CAM/SEMA4F in EC

We classified the EC data in the TCGA database into groups based on the expression levels of four neural-related genes: CHRM2, GRIN1, L1CAM, and SEMA4F, and performed GO and KEGG analyses on each group. The definitions of high and low expression were the same as those mentioned previously. We gathered some functional pathways and used the ssGSEA algorithm to calculate functional pathway scores.

In the EC with high CHRM2 expression, 92 genes were upregulated and 24 genes were downregulated (FC>2, P<0.05) (Figures 6A, B). Compared with the CHRM2 low expression group, the CHRM2 highly-expressed group contained upregulation of tumor-promoting pathways including endometrial cancer, Wnt signaling pathway, and estrogen signaling pathway. It also inhibits endometrial carcinogenesis and metastasis by activating the PPAR signaling system, suppressing the cAMP signaling pathway, cell adhesion molecules, vasculogenesis, amino acid synthesis and metabolism, and mucin-type O−glycan biosynthesis, among other things (Figures 6C, D). ssGSEA analysis showed that the CHRM2 gene was negatively associated with cellular response to hypoxia, tumor proliferation signature, DNA repair, G2M checkpoint, MYC targets, IL-10 anti-inflammatory signaling pathway, DNA replication, and positively associated with EMT markers, ECM-related genes, angiogenesis, TGFB, collagen formation (Figure 7).

FIGURE 6

Figure 6 Differential expression and enrichment analysis of CHRM2 high and low expression groups. (A) The volcano plot shows the differential gene expression of CHRM2 high expression group and CHRM2 low expression group was drawn with fold-change values and adjusted P. (B) Differential gene expression showed by heatmap (only 50 genes were displayed because of the large quantity of the genes); (C, D) KEGG and GO analysis showed the upregulated/downregulated pathways of the CHRM2 high expression group compared with the low expression group. When P<0.05 or FDR<0.05 is considered to be enriched to a meaningful pathway.

FIGURE 7

Figure 7 Analysis of the correlation between CHRM2 and 19 pathways using the Spearman, ssGSEA algorithm. (Statistically significant difference at P < 0.05). (A, G, J–L, P, S) The ssGSEA analysis showed that the CHRM2 gene was not statistically correlated with tumor inflammation signature, apoptosis, infammatory response, PI3K-AKT-MTOR pathway, p53 pathway, gene up-regulated by reactive oxigen specis (ROS) and degradation of ECM. (B, C, H, I, M, O, Q) CHRM2 gene was negatively associated with cellular response to hypoxia, tumor proliferation signature, DNA repair, G2M checkpoint, MYC targets, IL-10 anti-inflammatory signaling pathway and DNA replication. (D–F, N, R) CHRM2 gene was positively associated with EMT markers, ECM-related genes, angiogenesis, TGFB and collagen formation.

In the EC with high GRIN1 expression, 465 genes were upregulated and 264 genes were downregulated (FC>2, P<0.05) (Figures S3A, S3B). Compared with the GRIN1 low expression group, the GRIN1 high expression group inhibited the Wnt signaling pathway, signaling pathway regulating pluripotency of stem cells, cell adhesion molecules, hepatocellular carcinoma, and gastric cancer. Activation of cancer-related signaling pathways: p53 signaling pathway, cAMP signaling pathway, AMPK signaling pathway, breast cancer, and prostate cancer (Figures S3C, S3D). ssGSEA analysis showed that GRIN1 was negatively associated with cellular response to hypoxia, tumor proliferation signature, apoptosis, DNA repair, G2M checkpoint, inflammatory response, MYC targets, TGFB, IL-10 anti-inflammatory signaling pathway, DNA replication, collagen formation, and ECM degradation (Figure S4A).

In the EC with high L1CAM expression, 568 genes were upregulated and 493 genes were downregulated (FC>2, P<0.05) (Figures S5A, S5B). Compared with the L1CAM low expression group, the L1CAM highly-expressed group is consisted of upregulation of tumor-promoting pathways including PI3K-Akt signaling pathway, MAPK signaling pathway, cell adhesion molecules, and synapse organization. The p53 signaling pathway, Wnt signaling pathway, estrogen signaling pathway, and endometrial cancer pathway were inhibited (Figures S5C, S5D). ssGSEA analysis showed that L1CAM positively connected with cellular response to hypoxia, tumor proliferation signature, DNA repair, G2M checkpoint, MYC targets, TGFB, IL-10 anti-inflammatory signaling pathway, DNA replication, collagen production, and ECM degradation. reactive oxygen species (ROS) upregulation of genes was found to be negatively linked (Figure S4B).

In the EC with high SEMA4F expression, 11 genes were upregulated and 62 genes were downregulated (FC>2, P<0.05) (Figures S6A, S6B). Compared to the SEMA4F low expression group, The SEMA4F high expression group activated multiple diseases of neurodegeneration, Parkinson’s disease, and Alzheimer’s disease while inhibiting the IL-17 signaling pathway, chemical carcinogenesis, protein-coupled receptor signaling pathway, axoneme assembly, and acute inflammatory response (Figures S6C, S6D). ssGSEA analysis showed that SEMA4F was negatively correlated with tumor inflammation signature, ECM-related gene, angiogenesis, apoptosis, inflammatory response, P53 pathway, IL-10 anti-inflammatory signaling pathway, genes up-regulated by reactive oxygen species (ROS), tumor proliferation signature, DNA Repair, G2M checkpoint, MYC targets, TGFB, and DNA replication (Figure S4C).

The findings imply that the genes CHRM2, GRIN1, L1CAM, and SEMA4F in EC have pro- or oncogenic effects as a result of the combined activity of several signaling pathways influencing tumor growth.

Correlation between CHRM2/GRIN1/L1CAM/SEMA4F expression and immune infiltration, immune response and stemness

We used the immunedeconv R package to obtain immune infiltration data for high/low expression of CHRM2, GRIN1, L1CAM, and SEMA4F in EC. In this research, the CIBERSORT and EPIC algorithms were applied.

The CIBERSORT algorithm showed that high expression of CHRM2 was positively correlated with T cell CD4+ memory resting (P<0.05) and T cell regulatory (Tregs) (P<0.001), myeloid dendritic cell resting (P<0.05), mast cell activated (P<0.001), and mast cell resting (P<0.05), while it was negatively correlated with NK cell activated (P<0.05), and macrophage M1 (P<0.05) were negatively correlated (Figure 8A). The EPIC algorithm showed that high expression of CHRM2 was positively correlated with T cell CD4+ (P<0.001), T cell CD8+ (P<0.001), and endothelial cell (P<0.001), but negatively correlated with macrophage (P<0.05) (Figure 8B).

FIGURE 8

Figure 8 An analysis of immune infiltration, immunological response, and stemness in two groups with high and low CHRM2 expression. (A, B) Comparison of CHRM2 high expression group and CHRM2 low expression group in immune infiltration obtained with CIBERSORT and EPIC algorithm; The horizontal axis represents different immune cells, the vertical axis represents the immune scores (*P<0.05, ***P<0 .001). (C) Comparison immune checkpoint genes expression in CHRM2 high expression group and CHRm2 low expression group; The horizontal axis represents different immune checkpoint genes, the vertical axis represents the expression level (*P<0.05,**p<0.01). (D) Statistical table of immune response and the distribution of immune response scores of the different groups in predict results. (*P<0 .05). (E) Comparison of CHRM2 high expression group and CHRM2 low expression group in stemness was exhibited by mRNAsi score with OCLR algorithm (****p<0.0001).

The CIBERSORT algorithm showed that high GRIN1 expression was positively correlated with T cell CD8+ (P<0.01), T cell regulatory (Tregs) (P<0.001), and NK cell resting (P<0.01) compared to low GRIN1 expression (P<0.001), while it was negatively correlated with B cell plasma (P<0.05), NK cell activated (P<0.05), macrophage M2 (P<0.05), and myeloid dendritic cell activated (P<0.05) (Figure S7A). The EPIC algorithm showed that high expression of GRIN1 was positively correlated with T cell CD4+ (P<0.001), T cell CD8+ (P<0.001), endothelial cell (P<0.05), but negatively correlated with macrophage (P<0.01) (Figure S7B).

The CIBERSORT algorithm showed that high L1CAM expression was positively correlated with T cell follicular helper (P<0.01), NK cell activated (P<0.01), macrophage M1 (P<0.001), and myeloid dendritic cell activated (P<0.01), while it was negatively correlated with T cell CD8+ (P<0.05), T cell CD4+ memory resting (P<0.05), T cell regulatory (Tregs) (P<0.001), NK cell resting (P<0.01), macrophage M2(P<0.05), myeloid dendritic cell resting (P<0.001) and neutrophil (P<0.05) (Figure S8A). The EPIC algorithm showed that high expression of L1CAM was positively correlated with B cell (P<0.05) and with T cell CD4+ (P<0.001) and T cell CD8+ (P<0.001), but endothelial cell (P<0.01) negatively (Figure S8B).

The CIBERSORT algorithm showed that high SEMA4F expression was positively correlated with B cell naive(P<0.01), myeloid dendritic cell resting(P<0.01), myeloid dendritic cell activated(P<0.001), and mast cell activated(P<0.05), while it was negatively correlated with T cell CD8+(P<0.001), T cell CD4+ memory activated(P<0.01), T cell regulatory (Tregs) (P<0.001) and neutrophil(P<0.05) (Figure S9A). The EPIC algorithm showed that high expression of SEMA4F was positively correlated with T cell CD4+ (P<0.01), but macrophage (P<0.01) and NK cell(P<0.001) negatively (Figure S9B).

In addition, we analyzed the correlation between ICG and the expression of CHRM2, GRIN1, L1CAM, and SEMA4F. CHRM2 expression was positively correlated with CTLA4 (P<0.05) and SIGLEC15 (P<0.01) (Figure 8C); GRIN1 expression was positively correlated with CTLA4 (P<0.05), while negatively correlated with PDCD1LG2 (P<0.01) (Figure S7C); L1CAM expression was positively correlated with CD274 (P<0.05), LAG3 (P<0.001) and PDCD1LG2 (P<0.05), while negatively correlated with CTLA4 (P<0.001) (Figure S8C). SEMA4F was positively correlated with HAVCR2 (P<0.01), LAG3 (P<0.001), PDCD1 (P<0.001) and TIGIT (P<0.01) (Figure S8C), while negatively correlated with CTLA4 (P<0.001), HAVCR2 (P<0.01), LAG3 (P<0.05), PDCD1 (P<0.001), TIGIT (P<0.001) and SIGLEC15 (P <0.05) (Figure S9C). The TIDE algorithm revealed that high expression of CHRM2, and L1CAM was linked to poor immune response (Figures 8D, S8D), whereas high expression of GRIN1 was linked to a positive immunological response (Figure S7D). Stem cell scores were lower in the high expression group of CHRM2 (Figure 8E) and GRIN1 (Figure S7E) than in the low expression group, according to Spearman correlation analysis of OCLR scores, while the inverse was true for L1CAM (Figure S8E) and SEMA4F (Figure S9E).

Gene landscape of CHRM2/GRIN1/L1CAM/SEMA4F

We obtained mutational, transcriptomic, and clinical data of EC patients from the TCGA database and found the highest rate of PETN mutations in EC (57%), with varying degrees of mutations in neuro-oncology-related genes CHRM2, GRIN1, L1CAM, and SEMA4F: L1CAM (9%), SEMA4F (6%), CHRM2 (5%), and GRIN1 (3%) (Figure 9). The differences between GRIN1 and tumor mutational load (TMB) and microsatellite instability (MSI) were statistically significant, while L1CAM was negatively connected with TMB and SEMA4F was positively correlated with MSI. The findings imply that GRIN1, L1CAM, and SEMA4F are closely linked to immunotherapy and can respond to immunotherapy characteristics (Figure 10).

FIGURE 9

Figure 9 Mutational landscape of CHRM2, GRIN1, L1CAM and SEMA4F. (A) CHRM2. (B) GRIN1. (C) L1CAM. (D) SEMA4F.

FIGURE 10

Figure 10 The correlation of CHRM2, GRIN1, L1CAM, and SEMA4F with TMB and MSI using Spearman’s technique. (Statistically significant difference at P < 0.05).

Discussion

In the tumor microenvironment, neurons are critical biological components. Denervation and regulation of neurotransmitters for tumor treatment have become hot topics of research in recent years (27). Tumors select neuronal programs to promote their development and progression. The frequency of endometrial cancer has been increasing each year as is the number of patients with endometrial cancer brain metastases (28). Although nerve-cancer crosstalk influences tumor growth, the etiology is yet unknown (7).

The identification of cancer subgroups based on gene expression has proven useful in clinical settings, such as endometrial cancer molecular staging (29). Based on NRGs, we classified endometrial cancer into two subtypes: C1 and C2. Prognosis, clinical-stage, pathological grading, and immunological status were all statistically different between the two subtypes. In comparison to C2, C1 had a lower clinical stage and pathological grade, a better prognosis, better immune activation, stronger immune checkpoint gene expression, and was more suited to immunotherapy. In addition, there were statistical differences in enriched pathways and biological processes between C1 and C2. These data support the link between neural-related genes and endometrial cancer, and they suggest that using NRGs to classify EC subtypes could be clinically effective.

A prognostic model was constructed by LASSO-Cox, and CHRM2, GRIN1, L1CAM, and SEMA4F were identified as EC prognostic-related genes. CHRM2 is a gene encoding muscarinic receptor (mAChR) on neuronal cell membranes, which affects cholinergic activity by influencing the transcription level, mRNA stability, and affinity of the receptor (30). Previous studies have found that CHRM2 is enriched in the PI3K-Akt signaling pathway and its methylation rate rises as a progression of gastric cancer (31). CHRM2 inhibits the invasion and migration of non-small cell lung cancer through the M2R/ERK/Akt/NF-κB axis (32). In the central nervous system (CNS), glutamate receptor subunit 1 (GRIN1) is essential for synaptic transmission and plasticity (33). GRIN1 mutations are linked to schizophrenia, neurodevelopmental delay, epilepsy, and glioma, but no other tumors are linked to them (34–36). For the first time, our findings reveal that CHRM2 and GRIN1 play significant roles in endometrial cancer and are positively related to endometrial cancer prognosis. L1 cell adhesion molecule (L1CAM) is a membrane glycoprotein of the immunoglobulin family (37). Consistent with previous findings, L1CAM plays a key role in EC cancer cell migration and adhesion (38). Furthermore, we discovered that L1CAM was linked to immune cell infiltration and ICG (CD274, LAG3, PDCD1LG2, CTLA4), and we hypothesized that L1CAM may have a regulatory role in the tumor microenvironment, influencing tumor growth and metastasis. SEMA4F is a membrane-bound glycoprotein of the signaling element receptor family that has been linked to cancer in prior research, including being associated to breast cancer development (39), axonogenesis and neurogenesis in prostate cancer (40), and glioma prognosis (41). We propose that SEMA4F is a key regulator of tumor growth, angiogenesis, migration, and apoptosis and that it plays a role in endometrial cancer.

In summary, our study reveals the relevance of neural-related genes to endometrial cancer. Our findings suggest that EC reclassification based on neural-related genes is expected to be translated into clinical applications. The genes CHRM2, GRIN1, L1CAM, and SEMA4F, which are prognostically associated with endometrial cancer, play important roles in immune cell infiltration, immune response and stem cell relevance, clinical features, enriched pathways, and immunotherapy, and are potential biomarkers for EC with significant clinical translational potential. Further investigation can considered to quantify the indicators through tissue specimens and animal experiments to validate them for greater application in the treatment of tumors.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Author contributions

FC and YZ did the analysis. FC wrote the paper. YD and PL did the data sorting and charting. TQ and WJ conceived the paper. All authors contributed to the article and approved the submitted version.

Funding

This work was in part supported by the Innovation and Entrepreneurship Talent Project of Lanzhou (2020-RC-52) and the “Innovation Fund” for graduate students of the First Clinical Medical College of Gansu University of Chinese Medicine (LCCX2021011).

Acknowledgments

We apologize to those colleagues whose important work could not be cited due to space constraints.

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.

Supplementary material

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

Supplementary Figure 1 | Comparison of differential expression of 42 neural-related genes between C1 and C2.

Supplementary Figure 2 | Comparison of the clinical characteristics of the two groups with high and low expression of CHRM2, GRIN1, L1CAM and SEMA4F. (A) CHRM2. (B) GRIN1. (C) L1CAM. (D) SEMA4F.

Supplementary Figure 3 | Differential expression and enrichment analysis of GRIN1 high expression group and GRIN1 low expression group.

Supplementary Figure 4 | The connection of GRIN1, L1CAM, and SEMA4F with 19 pathways was investigated using the Spearman, ssGSEA algorithm. (S4A) GRIN1.(A) GRIN1. (B) L1CAM. (C) SEMA4F.

Supplementary Figure 5 | Differential expression and enrichment analysis of L1CAM high expression group and L1CAM low expression group.

Supplementary Figure 6 | Differential expression and enrichment analysis of SEMA4F high expression group and SEMA4F low expression group.

Supplementary Figure 7 | Analysis of immune infiltration, immunological response, and stemness in two groups with high and low GRIN1 expression.

Supplementary Figure 8 | Analysis of immune infiltration, immunological response, and stemness in two groups with high and low L1CAM expression.

Supplementary Figure 9 | Analysis of immune infiltration, immunological response, and stemness in two groups with high and low SEMA4F expression.

References

1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2021. CA: Cancer J Clin (2021) 71:7–33. doi: 10.3322/caac.21654

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: Cancer J Clin (2021) 71:209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Xia C, Dong X, Li H, Cao M, Sun D, He S, et al. Cancer statistics in China and united states, 2022: profiles, trends, and determinants. Chin Med J (2022) 135:584–90. doi: 10.1097/CM9.0000000000002108

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Faulkner S, Jobling P, March B, Jiang CC, Hondermarck H. Tumor neurobiology and the war of nerves in cancer. Cancer Discovery (2019) 9:702–10. doi: 10.1158/2159-8290

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Silverman DA, Martinez VK, Dougherty PM, Myers JN, Calin GA, Amit M. Cancer-associated neurogenesis and nerve-cancer cross-talk. Cancer Res (2021) 81:1431–40. doi: 10.1158/0008-5472.CAN-20-2793

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Magnon C, Hall SJ, Lin J, Xue X, Gerber L, Freedland SJ, et al. Autonomic nerve development contributes to prostate cancer progression. Sci (New York N.Y.) (2013) 341:1236361. doi: 10.1126/science.1236361

CrossRef Full Text | Google Scholar

7. Zhao C-M, Hayakawa Y, Kodama Y, Muthupalani S, Westphalen CB, Andersen GT, et al. Denervation suppresses gastric tumorigenesis. Sci Trans Med (2014) 6:250ra115. doi: 10.1126/scitranslmed.3009569

CrossRef Full Text | Google Scholar

8. Ni T, Huang T, Gu S-L, Wang J, Liu Y, Sun X, et al. DRG neurons promote perineural invasion of endometrial cancer via GluR2. J Cancer (2020) 11:2518–28. doi: 10.7150/jca.40055

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Español P, Luna R, Soler C, Caruana P, Altés-Arranz A, Rodríguez F, et al. Neural plasticity of the uterus: New targets for endometrial cancer? Women's Health (London England) (2022) 18:17455057221095537. doi: 10.1177/17455057221095537

CrossRef Full Text | Google Scholar

10. Talhouk A, McConechy MK, Leung S, Li-Chang HH, Kwon JS, Melnyk N, et al. A clinically applicable molecular-based classification for endometrial cancers. Br J Cancer (2015) 113:299–310. doi: 10.1038/bjc.2015.190

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Liu J, Lichtenberg T, Hoadley KA, Poisson LM, Lazar AJ, Cherniack AD, et al. An integrated TCGA pan-cancer clinical data resource to drive high-quality survival outcome analytics. Cell (2018) 173:400–16.e11. doi: 10.1016/j.cell.2018.02.052

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zhang Y-G, Jin M-Z, Zhu X-R, Jin W-L. Reclassification of hepatocellular cancer with neural-related genes. Front In Oncol (2022) 12:877657. doi: 10.3389/fonc.2022.877657

CrossRef Full Text | Google Scholar

13. Giorgi FM, Ceraolo C, Mercatelli D. The r language: An engine for bioinformatics and data science. Life (Basel Switzerland) (2022) 12:648. doi: 10.3390/life12050648

CrossRef Full Text | Google Scholar

14. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinf (Oxford England) (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

CrossRef Full Text | Google Scholar

15. Villanueva RAM, Chen ZJ. Ggplot2: Elegant graphics for data analysis (2nd Ed.). Measurement: Interdiscip Res Perspect (2019) 17:160–7. doi: 10.1080/15366367.2019.1565254

CrossRef Full Text | Google Scholar

16. David CC, Jacobs DJ. Principal component analysis: a method for determining the essential dynamics of proteins. Methods In Molecular Biology (Clifton, N.J.) (2014) 1084:193–226. doi: 10.1007/978-1-62703-658-0_11

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Software (2010) 33:1–22. doi: 10.18637/jss.v033.i01

CrossRef Full Text | Google Scholar

19. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Gene Ontology Consortium. Nat Genet (2000) 25:25–9. doi: 10.1038/75556

CrossRef Full Text | Google Scholar

20. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res (2000) 28:27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an r package for comparing biological themes among gene clusters. Omics J Integr Biol (2012) 16:284–7. doi: 10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

22. Abd ElHafeez S, D'Arrigo G, Leonardis D, Fusaro M, Tripepi G, Roumeliotis S. Methods to analyze time-to-Event data: The cox regression analysis. Oxid Med Cell Longevity (2021) 2021:1302811. doi: 10.1155/2021/1302811

CrossRef Full Text | Google Scholar

23. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Racle J, Gfeller D. EPIC: A tool to estimate the proportions of different cell types from bulk gene expression data. Methods In Mol Biol (Clifton N.J.) (2020) 2120:233–48. doi: 10.1007/978-1-0716-0327-7_17

CrossRef Full Text | Google Scholar

25. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24:1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, et al. Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell (2018) 173:338–54. doi: 10.1016/j.cell.2018.03.034

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Sartori R, Hagg A, Zampieri S, Armani A, Winbanks CE, Viana LR, et al. Perturbed BMP signaling and denervation promote muscle wasting in cancercachexia. Science Translational Medicine (2021) 13:eaay9592. doi: 10.1126/scitranslmed.aay9592

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Bhambhvani HP, Zhou O, Cattle C, Taiwo R, Diver E, Hayden Gephart M. Brain metastases from endometrial cancer: Clinical characteristics, outcomes, and review of the literature. World Neurosurg (2021) 147:e32–9. doi: 10.1016/j.wneu.2020.11.087

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Vermij L, Smit V, Nout R, Bosse T. Incorporation of molecular characteristics into endometrial cancer management. Histopathology (2020) 76:52–63. doi: 10.1111/his.14015

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Luo X, Kranzler HR, Zuo L, Zhang H, Wang S, Gelernter J. CHRM2 variation predisposes to personality traits of agreeableness and conscientiousness. Hum Mol Genet (2007) 16:1557–68. doi: 10.1093/hmg/ddm104

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Wang J, Ding Y, Wu Y, Wang X. Identification of the complex regulatory relationships related to gastric cancer from lncRNA-miRNA-mRNA network. J Cell Biochem (2020) 121:876–87. doi: 10.1002/jcb.29332

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Zhao Q, Yue J, Zhang C, Gu X, Chen H, Xu L. Inactivation of M2 AChR/NF-κB signaling axis reverses epithelial-mesenchymal transition (EMT) and suppresses migration and invasion in non-small cell lung cancer (NSCLC). Oncotarget (2015) 6:29335–46. doi: 10.18632/oncotarget.5004

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Nishimura N, Kumaki T, Murakami H, Enomoto Y, Katsumata K, Toyoshima K, et al. Arthrogryposis multiplex congenita with polymicrogyria and infantile encephalopathy caused by a novel variant. Hum Genome Variation (2020) 7:29. doi: 10.1038/s41439-020-00116-8

CrossRef Full Text | Google Scholar

34. Krzystanek M, Asman M, Witecka J, Pałasz A, Wiaderkiewicz R. Selected single-nucleotide variants in GRIN1, GRIN2A, and GRIN2B encoding subunits of the NMDA receptor are not biomarkers of schizophrenia resistant to clozapine: exploratory study. Pharmacol Rep PR (2021) 73:309–15. doi: 10.1007/s43440-020-00165-4

CrossRef Full Text | Google Scholar

35. Allen AS, Berkovic SF, Cossette P, Delanty N, Dlugos D, Eichler EE, et al. De novo mutations in epileptic encephalopathies. Nature (2013) 501:217–21. doi: 10.1038/nature12439

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Yang A, Wang X, Hu Y, Shang C, Hong Y. Identification of hub gene GRIN1 correlated with histological grade and prognosis of glioma by weighted gene coexpression network analysis. BioMed Res Int (2021) 2021:4542995. doi: 10.1155/2021/4542995

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Altevogt P, Doberstein K, Fogel M. L1CAM in human cancer. Int J Cancer (2016) 138:1565–76. doi: 10.1002/ijc.29658

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Bednarikova M, Vinklerova P, Gottwaldova J, Ovesna P, Hausnerova J, Minar L, et al. The clinical significance of DJ1 and L1CAM serum level monitoring in patients with endometrial cancer. J Clin Med (2021) 10. doi: 10.3390/jcm10122640

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Gabrovska PN, Smith RA, Tiang T, Weinstein SR, Haupt LM, Griffiths LR. Semaphorin-plexin signalling genes associated with human breast tumourigenesis. Gene (2011) 489:63–9. doi: 10.1016/j.gene.2011.08.024

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Ayala GE, Dai H, Powell M, Li R, Ding Y, Wheeler TM, et al. Cancer-related axonogenesis and neurogenesis in prostate cancer. Clin Cancer Res an Off J Am Assoc For Cancer Res (2008) 14:7593–603. doi: 10.1158/1078-0432.CCR-08-1164

CrossRef Full Text | Google Scholar

41. Shergalis A, Bankhead A, Luesakul U, Muangsin N, Neamati N. Current challenges and opportunities in treating glioblastoma. Pharmacol Rev (2018) 70:412–45. doi: 10.1124/pr.117.014944

PubMed Abstract | CrossRef Full Text | Google Scholar