Mapping Single-Cell Transcriptomes of Endometrium Reveals Potential Biomarkers in Endometrial Cancer

The endometrium is the layer forming the inner wall of the uterus of mammals. It reacts to both estrogen and progesterone, and therefore changes significantly with the estrus cycle and menstrual cycle. Dysfunctions of endometrium had been associated with various human diseases, including abnormal uterine bleeding, infertility, pre-eclampsia, endometriosis and endometrial carcinoma (EC). 1 , 2 Thus, understanding the gene expression regulation during menstrual cycle in humans is crucial for understanding normal functions of endometrium and the mechanism of EC. It also provides important insights for precision medicine in endometrium-related diseases. Significant efforts have been made to understand the pathogenesis and identify biomarkers for endometrium-related diseases, such as endometrial cancer. Anna et al found that ncRNA plays a significant role in endometrial cancer and is associated with different clinical prognostic outcomes. 3 These results suggested that, in addition to clinicopathological stratification, the inherent genetic and epigenetic characteristics of EC patients are also crucial factors in the precision medicine of oncology. 4 With the development of high throughput sequencing technologies, single-cell RNA sequencing (scRNA-seq) as a revolutionary technology can uncover novel cell types, explore genetic and functional heterogeneity in various cellular contexts. Compared to traditional bulk RNA sequencing data, which only assesses the average expression levels of all cells in a mixed tissue, scRNA-seq performs sequencing at the level of individual cells, offering higher resolution. This not only reveals cell types and heterogeneity but also offers a deeper understanding of specific cell functions and dynamic changes, delivering more detailed and precise insights than traditional approaches. Recently, scRNA-seq has been used to understandthe potential mechanisms of reproductive diseases. Wang et al utilized scRNA-seq to analyze the cellular and molecular signatures of decidual and peripheral leukocytes in normal and unexplained recurrent miscarriage. 5 Liu et al performed a scRNA-seq analysis of adenomyosis and supported the theory of adenomyosis derived from the invasion and migration of the endometrium. 6 Ma et al performed single-cell analysis and identified nine cell types, and determined a potential developmental trajectory associated with endometriosis. 7 One recent study mapped the temporal and spatial dynamics of the human endometrium in vivo and in vitro. 8 Guo et al presented a view of endometrial carcinoma at single-cell resolution and revealed the characteristics of endometrial epithelial cells in the endometrium. 9 However, the dynamic immune microenvironment in human endometrium, including cellular composition, cellular states, and cell–cell interactions, remains unclear. Moreover, EC is the most commonly diagnosed gynecologic malignancy. 10 Traditional classification of EC is primary based either on clinical and endocrine features. 10 , 11 Histological subtyping is commonly used in clinical to guide prognosis and treatment decisions for EC patients, while ongoing researches are evaluating the potential molecular subtyping. An integrated genomic analysis had resulted in the molecular classification of endometrioid and serous carcinomas into four distinct subgroups. 12 The subtypes identified by the different classification systems correlate to some extent; however, the potential molecular pathways of different subtypes of EC are still unknown. Therefore, we interrogate the immune microenvironment of human endometrial cells during the proliferative and secretory phases of women menstrual cycle. We identified the diverse cell types in human endometrial and characterization of the RNA expression patterns in different cell types. In particular, we explored the transcriptional regulation and revealed the transcription factors (TFs) that exhibited high activities in cell types. Cell–cell communications were examined and further uncovered five RNA-based molecular subtypes of EC.

To analyze the scRNA-seq cells from human endometrium through four menstrual cycle stages, we performed principal-component analysis (PCA) using the top 2000 most variably expressed genes across 59,397 cells. Cells were clustered into transcriptionally distinct clusters with top 20 principal components (PCs). The cells were visualized using UMAP plot and revealed eight clusters that could be annotated to known cell types ( Figure 1A ). Moreover, we found that cells were clustered together based on cell types but not based on patient identify ( Figure 1B ). We then used well-known marker genes to define the identity of each cell cluster. For example, epithelial cells expressed KRT8 and PAEP, endothelial cells expressed VWF, fibroblast cells expressed APOD, DCN and COL3A1, perivascular (PV) cells expressed RGS5, smooth muscle cells expressed ACTG2 and MYH11, multi-potent stromal cells (MSC) expressed TOP2A and UBE2C, lymphoid and myeloid cells expressed CD74, NKG7 and GNLY ( Figure 1A , and E and Figure S1 ). Figure 1 Diversity of cell types in human endometrium delineated by single-cell transcriptomic analysis. ( A and B ) UMAP projection of cells from 10 human endometrium samples. ( A ) colored by cell types and ( B ) colored by patients. ( C ) The cell fractions of different cell types originating from different phases of menstrual cycle. ( D ) The cell fractions of different cell types originating from different patients. ( E ) Bubble plots showing the expression levels of genes highly expressed in corresponding cell types. Diversity of cell types in human endometrium delineated by single-cell transcriptomic analysis. ( A and B ) UMAP projection of cells from 10 human endometrium samples. ( A ) colored by cell types and ( B ) colored by patients. ( C ) The cell fractions of different cell types originating from different phases of menstrual cycle. ( D ) The cell fractions of different cell types originating from different patients. ( E ) Bubble plots showing the expression levels of genes highly expressed in corresponding cell types. Next, we calculated the proportion of cells during four menstrual cycle stages and different patients. We found that the proportions of cell types among stages and patients were significantly different ( Figure 1C and D , p-values < 2.2E-16). PV cells were predominated in early-secretory stage and in A30 patient. Fibroblast cells decreased in early-secretory stage and immune cells were enriched in proliferative stage and late-secretory stage ( Figure 1C ). MSC decreased during the menstrual cycle stages. Moreover, we identified the highly expressed marker genes in each cell types. We found that C1QA, C1QB and C1QC were highly expressed in myeloid cells ( Figure 1E ), indicating their phagocytic ability. 28 Together, our comprehensive analysis provided a comprehensive catalog of the major cell types together with their cellular position in endometrium. We next explored the cell type-specific RNAs that could help explain distinct biological states of these cell types. Functional enrichment analysis revealed that genes highly expressed in fibroblast cells were significantly enriched in wound healing, regulation of vasculature development and regulation of angiogenesis ( Figure 2A ). Genes highly expressed in smooth muscle cells were enriched in wound healing, extracellular matrix organization and extracellular structure organization, whereas PV cell-specific genes were enriched in response to corticosteroid, steroid hormone and glucocorticoid ( Figure 2A ). Epithelial cell-specific genes were significantly enriched in epithelial cell proliferation and tissue migration, and endothelial cell-specific genes were enriched in regulation of vasculature development and angiogenesis ( Figure 2A ). MSCs are a population of self-renewing multipotent cells in the perivascular regions of the endometrium in both the basalis and functionalis. 29 , 30 We found that genes highly expressed in MSCs were enriched in sister chromatid segregation and nuclear division ( Figure 2A ). A large proportion of MSCs were in G2M stage of cell cycle ( Figure 2B ). In particular, we found that two proliferative marker genes, MKI67 and TOP2A, were highly expressed in MSC cells ( Figure 2C and D ). These results suggested that enriched functional analysis of each cluster supported their functions in human endometrium. Figure 2 Expression and functions of genes highly expressed in diverse cell types. ( A ) Heat map showing the expression levels of genes highly expressed in corresponding cell types. Enriched functions were indicated on the left. ( B ) UMAP plot, color-coded for cell cycle phases. ( C and D ) UMAP plot, color-coded for relative expression (lowest expression to highest expression, white to red) of marker genes. ( C ) for MKI67 and ( D ) for TOP2A. Expression and functions of genes highly expressed in diverse cell types. ( A ) Heat map showing the expression levels of genes highly expressed in corresponding cell types. Enriched functions were indicated on the left. ( B ) UMAP plot, color-coded for cell cycle phases. ( C and D ) UMAP plot, color-coded for relative expression (lowest expression to highest expression, white to red) of marker genes. ( C ) for MKI67 and ( D ) for TOP2A. We analyzed the PV population (n=8,120) based on the known markers and re-clustered it into two distinct populations as indicated in the UMAP ( Figure 3A and Figure S2A ). We analyzed the relative proportion of cells in each cluster and noted that two clusters exhibited similar proportion during menstrual cycle stages and patients ( Figure 3B , p=1.862E-8 and Figure S2B ). These observations were consistent with the results PV-MYH11+ are characteristic of myometrium, while PV-STEAP4+ are only present in the endometrium. 8 We next analyzed the gene expression profiles and identified the top differentially expressed genes in two PV populations ( Figure 3C ). In particular, STEAP4 and MYH11 were separately expressed in two PV populations ( Figure 3D and E ). We found that several collagen-related genes (eg COL3A1, COL1A2, COL4A1 and COL1A1) were highly expressed in PV-STEAP4 populations ( Figure 3C ), which might be correlated with their roles in repair of endometrial damage and induced angiogenesis. 31 Moreover, IGFBP5 was highly expressed in PV-STEAP4+ subtypes, which is consistent with its roles in promoting angiogenic and neurogenic differentiation. 32 , 33 In contrast, the PV-MYH11+ populations were characteristic of myometrium and highly expressed MUSTN1 and MYH11. Figure 3 PV cell clusters in human endometrium. ( A ) UMAP plot of PV cells color-coded according to the assigned cell subtypes. ( B ) The cell fractions of different PV cell types originating from different phases of menstrual cycle. ( C ) Bubble plots showing the expression levels of genes highly expressed in corresponding PV subtypes. ( D and E ) UMAP plot, color-coded for relative expression (lowest expression to highest expression, white to red) of marker genes. ( D ) for STEAP4 and ( E ) for MYH11. ( F ) Heat map showing the expression levels of genes highly expressed in corresponding PV cell subtypes. Enriched functions were indicated on the left. ( G ) Lollipop plot showing the differences of pathway activities between two PV subtypes. PV cell clusters in human endometrium. ( A ) UMAP plot of PV cells color-coded according to the assigned cell subtypes. ( B ) The cell fractions of different PV cell types originating from different phases of menstrual cycle. ( C ) Bubble plots showing the expression levels of genes highly expressed in corresponding PV subtypes. ( D and E ) UMAP plot, color-coded for relative expression (lowest expression to highest expression, white to red) of marker genes. ( D ) for STEAP4 and ( E ) for MYH11. ( F ) Heat map showing the expression levels of genes highly expressed in corresponding PV cell subtypes. Enriched functions were indicated on the left. ( G ) Lollipop plot showing the differences of pathway activities between two PV subtypes. To further determine the specific roles of two PV populations that might contribute to human endometrium, we performed functional enrichment analysis based on the differentially expressed genes. We found that genes highly expressed in two PV populations were both significantly enriched in wound healing, response to oxidative stress and extracellular matrix organization ( Figure S2C ). However, the proportions of genes in PV-STEAP4+ were much higher. In particular, genes highly expressed in PV-STEAP4+ were significantly enriched in extracellular matrix organization and structure organization, while genes in PV-MYH11+ were significantly enriched in muscle development-related functions ( Figure 3F ). Cancer hallmark-related pathways also exhibited distinct activities in two PV populations. PV-STEAP4+ cells exhibited higher activities in immune and metabolism-related functions and PV-MYH11+ cells exhibited higher activities in signaling and proliferative pathways, such as PI3K-AKT-mTOR ( Figure 3G ). It has been demonstrated that glycolysis is beneficial to angiogenesis and plays an important role in endometrial decidualization. 34–36 In addition, PI3K-AKT-mTOR pathway plays an important role in the decidualization of endometrium. The epithelial cell cluster (n=5,758) was then assigned to four epithelial subtypes based on known markers obtained from the published literatures. 37 , 38 Three secretory glandular cells (secretory LGALS1+/PAEP+/MT+) and ciliated cells were identified ( Figure 4A ). We calculated the proportion of four subpopulations across menstrual cycle and found that ciliated epithelial cells were enriched in proliferative phase, secretory LGALS1+ cells were enriched in early and mid-secretory phases ( Figure 4B , p<2.2E-16). By examining the gene expression patterns of four epithelial subpopulations, we revealed that ciliated cells expressed high levels of markers, such as the ciliated marker TPPP3, PIFO and FAM183A ( Figure 4C and D ). Figure 4 Epithelial cell clusters in human endometrium. ( A ) UMAP plot of epithelial cells color-coded according to the assigned cell subtypes. ( B ) The cell fractions of different epithelial cell types originating from different phases of menstrual cycle. ( C ) Bubble plots showing the expression levels of genes highly expressed in corresponding epithelial subtypes. ( D ) UMAP plot, color-coded for relative expression (lowest expression to highest expression, white to red) of marker genes. ( E ) Heat map showing the expression levels of genes highly expressed in corresponding epithelial cell subtypes. Enriched functions were indicated on the left. ( F ) Heat map showing the differences of pathway activities among epithelial subtypes. Epithelial cell clusters in human endometrium. ( A ) UMAP plot of epithelial cells color-coded according to the assigned cell subtypes. ( B ) The cell fractions of different epithelial cell types originating from different phases of menstrual cycle. ( C ) Bubble plots showing the expression levels of genes highly expressed in corresponding epithelial subtypes. ( D ) UMAP plot, color-coded for relative expression (lowest expression to highest expression, white to red) of marker genes. ( E ) Heat map showing the expression levels of genes highly expressed in corresponding epithelial cell subtypes. Enriched functions were indicated on the left. ( F ) Heat map showing the differences of pathway activities among epithelial subtypes. The gene set enrichment analyses showed that genes up-regulated in ciliated cells were mainly enriched for cilium organization, cilium assembly, cilium movement functions ( Figure 4E and Figure S3 ). This result is similar with one previous study. 9 The secretory MT+ cell populations highly expressed genes in metallothionein (MT) family, such as MT2A, MT1G, MT1E, MT1X and MT1H ( Figure 4E ), which might contribute to impaired endometrial receptivity. 39 Moreover, higher proportion of secretory MT+ cells was observed in late-secretory phase ( Figure 4B ). This is consistent with previous observations that MT genes were highly expressed in late-secretory phase. 40 The secretory PAEP+ cell populations were enriched during the menstrual cycle and highly expressed CLDN3, CLDN4 and CXCL family genes ( Figure 4E ). Functional analysis revealed that genes up-regulated in this cell population were enriched in regulation of cell–cell adhesion, cell migration and reproductive development ( Figure S3 ). We next performed functional analysis based on cancer hallmark-related pathways. By examining the pathway activities, we found that secretory populations exhibited higher pathway activities in numerous pathways ( Figure 4F ). The ciliate epithelial cells exhibited higher pathway activities in bile acid metabolism and mitotic spindle. The Wnt signaling pathway, DNA repair and G2M checkpoint pathways were enriched in secretory LGALS1+ cell populations ( Figure 4F ). In particular, epithelial mesenchymal transition was enriched in secretory LGALS1+, which was consistent with their roles in regeneration or differentiation of endometrium. 29 The inflammation-related pathways were significantly enriched in secretory PAEP+ cell populations, such as IL6-JAK-STAT3 and IL2-STAT5 signaling pathways ( Figure 4F ). These results extended the transcriptional signature and potential pathways underlying the human endometrial epithelial cells. Fibroblasts (n=21,865) were further separated into four clusters and annotated based on the highest expressing genes ( Figure 5A ). SPARACL1, ID4, MMP11 and EGR1 were highly expressed in corresponding clusters. We found that fibroblast ID4+ and SPARCL1+ cells were primarily observed in early- and late-secretory phases ( Figure 5B ). Fibroblast MMP11+ cells were primarily in the proliferative phase and EGR1+ cells were in mid-secretory phase ( Figure 5B , p<2.2E-16). The proportions were variable in different patients ( Figure S4A ), suggesting the high heterogeneity among patients. Gene expression analysis revealed fibroblast SPARCL1+ highly expressed fibroblast-related genes, such as DCN, APOD, COL1A2 and COL15A1 ( Figure 5C and Figure S4B ). In particular, four clusters highly expressed corresponding marker genes ( Figure 5D ). We next performed the functional enrichment analysis and found that genes up-regulated in fibroblast SPARCL1+ were significantly enriched in extracellular structure organization and negative regulation of locomotion ( Figure 5E ). Genes highly expressed in fibroblast ID4+ were enriched in epithelial cell proliferation, positive regulation of cytokine production and regulation of angiogenesis ( Figure 5E ). Numerous of transcription factors (TFs) were highly expressed in fibroblast EGR1+ sub-populations, such as JUNB, EGR1, FOS and JUN. We also observed high expression of GADD45B in this population, which plays important roles in DNA repair. 41 We also investigated the cancer hallmark pathway activities in different fibroblast populations and found that Wnt, MYC and peroxisome exhibited higher activities in ID4+ cell populations ( Figure S4C ). Notch and PI3K-AKT signaling pathways were active in SPARCL1+ cell populations ( Figure S4C ). Together, these results uncovered the potential pathways underlying different fibroblast populations in human endometrium. Figure 5 Fibroblast cell clusters in human endometrium. ( A ) UMAP plot of fibroblast cells color-coded according to the assigned cell subtypes. ( B ) The cell fractions of different fibroblast cell types originating from different phases of menstrual cycle. ( C ) Bubble plots showing the expression levels of genes highly expressed in corresponding fibroblast subtypes. ( D ) UMAP plot, color-coded for relative expression (lowest expression to highest expression, white to red) of marker genes. ( E ) Heat map showing the differences of pathway activities among fibroblast subtypes. Fibroblast cell clusters in human endometrium. ( A ) UMAP plot of fibroblast cells color-coded according to the assigned cell subtypes. ( B ) The cell fractions of different fibroblast cell types originating from different phases of menstrual cycle. ( C ) Bubble plots showing the expression levels of genes highly expressed in corresponding fibroblast subtypes. ( D ) UMAP plot, color-coded for relative expression (lowest expression to highest expression, white to red) of marker genes. ( E ) Heat map showing the differences of pathway activities among fibroblast subtypes. It has been demonstrated that genes usually interact with each other to form a complex interaction network. 42 We thus analyzed the correlation between gene expressions and we found that genes highly expressed in cell types were co-expressed with each other ( Figure 6A ), indicating the modular programs associated with basic cellular functions of cells. TFs are important regulators and play important roles in regulating gene expression. 43 In addition, we calculated the module activities as the average expression of genes. We found that the expressions of several TFs (such as IRF8, SOX17, FOS and KLF2) were significantly correlated with the module activities ( Figure S5 ). We thus identified the cell type-specific TFs based on the pySCENIC pipeline. In total, we identified 336 TFs exhibited high activities in different cell types and 62 TFs also exhibited differential expression ( Figure 6B ). We found that TFs exhibited high cell type- and phase-specificity. For example, SOX4 exhibited higher activity in proliferative phase ( Figure 6B ). SOX4 is closely associated with the development and progression of many malignant tumors and plays an important role in the cell growth and proliferation. 44 We also observed that the expressions of SOX4 in proliferative phase were significantly higher than secretory phases. In addition, MAF, KLF4, JUN, FOS and EGR1 exhibited higher activities in secretory phases ( Figure 6C ). Figure 6 TF regulators of cell types. ( A ) Co-expression of genes differentially expressed in cell type. ( B ) Heat maps showing the TF activities in each cell types. Each row represents a TF and column represents a cell. Numbers in the brackets were the potential target genes of TFs. ( C ) Violin plots showing the expression of TFs during menstrual cycle. ( D ) Bar plots showing the enriched functions of TF target genes. TF regulators of cell types. ( A ) Co-expression of genes differentially expressed in cell type. ( B ) Heat maps showing the TF activities in each cell types. Each row represents a TF and column represents a cell. Numbers in the brackets were the potential target genes of TFs. ( C ) Violin plots showing the expression of TFs during menstrual cycle. ( D ) Bar plots showing the enriched functions of TF target genes. Moreover, we found that several TFs exhibited cell type-specific activities in human endometrium. CD59 exhibited higher activities in endothelial cells, which is highly expressed in endothelial cells of human endometrium ( Figure S6 ). It has been demonstrated that CD59 may be important in protection of endothelial cells against C-mediated damage at local sites of inflammation, thereby maintaining the vascular integrity. 45 E2F1 and EZH2 exhibited higher activities in MSC cell populations ( Figure 6B ). E2F1 plays a pivotal role in driving cells out of a quiescent state and into the S phase of the cell cycle 46 and the activity of EZH2 influences cell fate regulation. 47 Functional analysis of the targets of EZH2 revealed that they were significantly enriched in DNA replication and cell cycle, which is consistent with the proliferative state of MSC cells. EGR1 exhibited higher activities in fibroblast cells ( Figure 6B ), which had been demonstrated to contribute to inflammatory factors and fibrosis reduction. 48 The potential targets of EGR1 were significantly enriched in Wnt signaling pathway and cell growth ( Figure 6D ). In addition, FOS exhibited higher activities in fibroblast and smooth muscle cell populations ( Figure 6B ) and its targets were significantly enriched in muscle tissue development, Wnt signaling and cell differentiation ( Figure 6D ). It has been demonstrated that FOS plays an important role in control of fibroblast senescence and activation of programmed cell death. 49 Several TFs also exhibited higher activities in immune cells, such as STAT4, BATF, RUNX3 and EOMES ( Figure 6B ). These results suggest that the dynamic transcriptome during menstrual cycle was strictly regulated by TFs that were with dynamic activities. The development of tissue and progression of complex diseases relies on a complex network of cell–cell interactions. 50 We thus inferred intercellular communications for human endometrium based on CellChat. 26 As a result, we found that the myeloid cells frequently interact with MSC and epithelial cells, while the fibroblast cells frequently interact with others ( Figure 7A ). Immune cells mainly function as signal input cells. As the signal output cells, epithelial and endothelial cells mainly interact with immune cells as targets ( Figure 7A ). Next, the 15 signaling pathways associated with inferred networks were mapped onto a two-dimensional manifold and clustered into four groups ( Figure 7B and C ). MIF, IGF, PTN and MK pathways were grouped into cluster 1, while group 4 was formed by pathways of VEGF, CXCL and VISFATIN ( Figure 7C ). Figure 7 Cell–cell communications in human endometrium. ( A ) Cell–cell communications. ( B ) Jointly projecting and clustering signaling pathways into a shared two-dimension manifold according to their functional similarity. The size is proportional to the total communication probability. ( C ) Magnified view of each pathway group. ( D ) The inferred MIF signaling pathway network in human endometrium. Left and right portions show the autocrine and paracrine signaling, respectively. Circle sizes are proportional to the number of cells in each cell group and edge width represents the communication probability. ( E ) Relative contribution of each ligand-receptor pair to the overall MIF signaling pathway network. ( F ) Expression distribution of MIF signaling genes. ( G ) The inferred SPP1 signaling pathway network in human endometrium. Left and right portions show the autocrine and paracrine signaling, respectively. Circle sizes are proportional to the number of cells in each cell group and edge width represents the communication probability. ( H ) Relative contribution of each ligand-receptor pair to the overall SPP1 signaling pathway network. ( I ) Expression distribution of SPP1 signaling genes. ( J and K ) The dot plot showing the comparison of incoming and outgoing signaling patterns of secreting cells. The dot size is proportional to the contribution score computed from pattern recognition analysis. Higher contribution score implies the signaling pathway is more enriched in the corresponding cell group. Cell–cell communications in human endometrium. ( A ) Cell–cell communications. ( B ) Jointly projecting and clustering signaling pathways into a shared two-dimension manifold according to their functional similarity. The size is proportional to the total communication probability. ( C ) Magnified view of each pathway group. ( D ) The inferred MIF signaling pathway network in human endometrium. Left and right portions show the autocrine and paracrine signaling, respectively. Circle sizes are proportional to the number of cells in each cell group and edge width represents the communication probability. ( E ) Relative contribution of each ligand-receptor pair to the overall MIF signaling pathway network. ( F ) Expression distribution of MIF signaling genes. ( G ) The inferred SPP1 signaling pathway network in human endometrium. Left and right portions show the autocrine and paracrine signaling, respectively. Circle sizes are proportional to the number of cells in each cell group and edge width represents the communication probability. ( H ) Relative contribution of each ligand-receptor pair to the overall SPP1 signaling pathway network. ( I ) Expression distribution of SPP1 signaling genes. ( J and K ) The dot plot showing the comparison of incoming and outgoing signaling patterns of secreting cells. The dot size is proportional to the contribution score computed from pattern recognition analysis. Higher contribution score implies the signaling pathway is more enriched in the corresponding cell group. We specifically examined how macrophage migration inhibitory factor (MIF) communications among cell populations ( Figure 7D ). All cells function as signal output cells, although the intensity of interaction was different. When MSC and epithelial cells acted as signal input cells, the output intensity of cells was low. In this pathway, immune cells were mainly used as signal input cells, and the strength of intercellular interaction was high ( Figure 7D ). Macrophage migration inhibitory factor had been identified as a potential biomarker of endometriosis. 51 , 52 We also examined ligand receptor pairs that play a major role in this pathway and identified two ligand-receptor pairs MIF−(CD74+CD44), and MIF−(CD74+CXCR4) ( Figure 7E ). MIF exhibited high expression in all cell types, while CXCR4, CD74 and CD44 exhibited higher expression in immune cells ( Figure 7F ). We also found that the SPP1 signaling pathway only took myeloid cells as signal output cells, and exported signals to fibroblasts, PV cells, SMC cells and lymphocytes in paracrine mode, and played a role in autocrine mode to a large extent ( Figure 7G ). The SPP1-(ITGA4+ITGB1) and SPP1-CD44 plays important roles in cell–cell communications ( Figure 7H ). SPP1 and CD44 exhibited higher expression in immune cells, while ITGB1 was highly expressed in all cell types ( Figure 7I ). Next, we compared the information flow for each signaling pathway. In signal input, several pathways only functions in one cell type, such as EGF and EDN in fibroblasts, VEGF and CXCL in endothelial cells, and ncWNT in PV cells ( Figure 7J ). MIF signaling pathway only participates in the signal input of immune cells and most of these pathways are involved in signal input patterns of myeloid cells, fibroblasts, and endothelial cells ( Figure 7J ). In the signal output part, we found that IL6, EGF, SPP1 and CCL only played a role in the signal output of immune cells. EDN only plays a role in signal output of endothelial cells. As with signal input patterns, most of these pathways were involved in signal output by fibroblasts and myeloid cells ( Figure 7K ). Together, cell–cell communication analysis enables multifaceted assessment of intercellular communication patterns in human endometrium. Endometrial cancer (EC) is the most common gynecologic malignancy. We next explored the expression patterns of the ligand-receptor in EC. We found that the patients can be grouped into five clusters based on the ligand-receptor of the 15 pathways identified in cell–cell communication ( Figure 8A ). The survival rates for patients in five clusters were significantly different, and patients in cluster-5 were with poor survival ( Figure 8B , Log rank test p=0.00025). We found that patients in cluster-5 were with distinct expression of 17 genes, which were involved multiple cytokine-related pathways, including MIF, CXCL, SPP1 and VEGF. We next calculated the pathway activities for these pathways and found that the majority of pathway activities in cluster-5 were significantly higher than other clusters. For example, epidermal growth factor (EGF) was a conventional mitogenic factor that can stimulate the proliferation of various types of cells including epithelial cells and fibroblasts. 53 The activities of EGF in cluster-5 were significantly higher than that in cluster 1, 3, and 4 ( Figure 8C ). The EGFR family also plays an important role in maintaining epithelial homeostasis 54 and overexpressed in EC and ovarian cancers. 55 , 56 Compared with other clusters, we found that ERBB2 was highly expressed in patients of cluster-5 ( Figure 8D ). Figure 8 Molecular subtypes of EC. ( A ) Heat map showing the consensus matrix of EC patients. ( B ) The overall survival curves based on EC patients data, stratified by the expression patterns of ligand-receptor pairs. ( C ) Boxplots showing the EGF pathway activities in patients of different EC clusters. ( D ) Boxplots showing the expressions of ERBB2 in patients of different EC clusters. ( E ) Boxplots showing the IGF pathway activities in patients of different EC clusters. ( F ) Boxplots showing the expressions of IGF2 in patients of different EC clusters. Molecular subtypes of EC. ( A ) Heat map showing the consensus matrix of EC patients. ( B ) The overall survival curves based on EC patients data, stratified by the expression patterns of ligand-receptor pairs. ( C ) Boxplots showing the EGF pathway activities in patients of different EC clusters. ( D ) Boxplots showing the expressions of ERBB2 in patients of different EC clusters. ( E ) Boxplots showing the IGF pathway activities in patients of different EC clusters. ( F ) Boxplots showing the expressions of IGF2 in patients of different EC clusters. In addition, as a major ligand of IGF, IGF2 overexpression has been shown to play a role in many cancers. 57 However, there seem to be few reports about its role in endometrial cancer. We found that the activities of IGF pathway and expressions of IGF2 in patients of cluster-5 were significantly higher than cluster-2 and cluster-4 ( Figure 8E and F ). In contrast, several pathways (ie CXCL) had lower pathway activities in cluster-5 and the expressions of CXCR4 were also lower in patients of cluster-5. These results suggested that the genes identified here can be used to predict the prognostic survival of patients with endometrial cancer.

Single-cell RNA-sequencing was used to create a cell census of the human endometrium. We downloaded the raw sequencing data from ArrayExpress under the accession number E-MTAB-10287. 8 In total, 11 samples obtained from five patients were sequenced, including four menstrual cycle stages (proliferative stage, early-secretory stage, mid-secretory stage and late-secretory stage). We selected 10 samples obtained from endometrium for further analysis. Gene expression profiles and clinical information of human endometrial cancers were obtained from The Cancer Genome Atlas (TCGA) project ( https://portal.gdc.cancer.gov/ ). The expressions of genes were measured by Fragments Per Kilobase of exon model per Million mapped fragments (FPKM). In total, 535 human endometrial cancers were included in our analysis. The 10x Genomics scRNA-seq data were first analyzed using Cellranger 6.1.1 with the raw fastq files as input. The GRCh38 genome was used as the reference genome and default parameters were used. For each sample, the feature-barcode matrix was then converted into a Seurat object using the Seurat R package. 13 To enrich for high-quality cells in each sample, we performed quality control (QC) for each sample dataset individually. First, we filtered cells that expressed less than 200 genes. For A13 and A30 patients, we excluded the cells with >10% mitochondrial reads. For E1, E2 and E3, cells with >20% mitochondrial reads were excluded. In addition, we filtered cells with hemoglobin protein-related reads >5%. Next, we used the “isOutlier” function in scater package to detect the outliers. The outlier cells were defined in each of the following metrics: log(UMI counts) (>2 MADs, both), log(number of genes expressed) (>2 MADs, both) and log(percent mitochondrial read count) (>2 MADs, high end). We used the “CellCycleScoring” in Seurat to calculate the cell cycle score for each cell and regression with the “vars.to.regress”. Cell cycle genes (G2/M and S) of were obtained from Seurat package. To reduce the false positive rate in doublet calling, only cells marked as doublets by both scDblFinder and doubletFinder 14 were removed from our analysis. The quality control was performed for single sample and the “CCA” function in Seurat was used to integrate all samples. The read count matrices were normalized using “NormalizedData” with “LogNormalize” as the normalization method. Feature selection was performed by “FindVariableFeatures” using the “vst” method and the top 2000 variable genes were identified. The top 2000 most variable genes were summarized by principal component analysis (PCA). To identify groups of distinct cells, graph-based Louvain clustering was performed based on top 20 PCs. The “FindClusters” function with a resolution of 0.2 was used to identify the cell clusters and UMAP plots were generated in R for visualization. We performed cell types annotation based on two methods. One was gene signature enrichment and another one was reference-based annotation with the SingleR package. 15 First, the cell types were annotated based on signatures from ESTIMATE 16 and PangladoDB. 17 The “AddModuleScore” was used to calculate the signature scores. The median scores of each cell clusters were calculated and if the median >0.1, we considered the clusters as corresponding cell types. The primary annotations were performed by SingleR using the “HumanPrimaryCellAtlasData” as reference dataset. The marker genes of cell types were obtained from literature or CellMarker database. 18 Differential gene expression analysis was performed by “FindAllMarkers” in Seurat with the min.pct set to 0.25 and Wilcoxon’s rank sum test. Genes with p.adiust0.25 were considered as up-regulated. For functional enrichment analysis, we selected top 50 highly expressed genes in each main cell type. For the cell subtypes, we used all differentially expressed genes. The functional enrichment analysis was performed by clusterProfiler, 19 and heat maps were generated by ComplexHeatmap. 20 To calculate the gene set functional scores, we used AUCell to perform this analysis. 21 The functional gene sets were obtained from the Molecular Signatures Database (MSigDB) hallmark gene set collection. 22 Wilcoxon’s rank sum test was used to compare the functional scores in different cell subtypes. P-values were adjusted by false discovery rate (FDR). The gene sets with log2FoldChange>0 and p.adjust<0.05 were considered as enriched in corresponding cell subtypes. To investigate the expression correlation among differentially expressed genes in each cell type, we calculated the Spearman correlation coefficient (SCC) among top 50 differentially expressed genes based on the expression across cell types. The SCC matrix was visualized by heat map. We used pySCENIC to identify the TF regulators in each cell types identified in human endometrium. 21 , 23 We first downloaded the motifs that allow using RcisTarget (mc9nr) from cisTarget database. 24 The input UMI count matrix was normalized (CPM) and log-transformed. Only genes in RcisTarget were included in further analysis. GENIE3 was performed to identify TF-modules in each cell type. 25 TF-modules having less than 10 genes were filtered out. The top 1 percentile of the number of detected genes per cell was used to calculate the AUCell enrichment of each TF regulon in each cell. ComplexHeatmap was used to generate the heat map of the activity matrix. Cell–cell communication analysis was performed using CellChat (Version 1.1.0), based on the known ligand-receptor pairs in CellChatDB. 26 Briefly, the normalized genes expression matrix and cell type labels generated by Seurat were subjected as input for CellChat. For the main analyses, the core functions “computeCommunProb”, “computeCommunProbPathway” and “mergeCellChat” were applied using default parameters. The “computeNetSimilarityPairwise” function was used to calculate the similarity between pathways and pathways were clustered into different groups based on functional similarities. To identify the molecular subtypes of human endometrial cancer, we performed clustering based on the gene expression. First, the ligand-receptor pairs prioritized in single-cell data analysis were used in the clustering. Nonnegative matrix factorization (NMF) was performed based on R package. 27 The parameters “rank=2:6, method=brunet” were used in this analysis. Five molecular subtypes were identified based on the cophenetic curve. The survival analysis was performed by R packages (survival and survminer). Log rank test was used to evaluate the difference of survival rates among subtypes.

In this study, our single transcriptome analysis yields a comprehensive catalog of the major subsets of uterine cells in endometrium. We comprehensively identified the diverse cell types in human endometrium and found that there were a higher proportion of fibroblast, PV and epithelial cells. We further characterized the gene expression patterns in different cell types and identified the cell-type specific genes. In particular, we explored the transcriptional regulation and revealed the TFs that exhibited high activities in specific cell types. In addition, cell–cell communications were examined and revealed the critical ligand-receptor pairs. The expression patterns of identified ligand-receptor further uncovered five molecular subtypes of EC. All these results increased our understanding of the dynamic microenvironment of uterine cells in endometrium, as well as potential biomarkers for EC. In particular, we identified two perivascular cell subtypes, four epithelial subtypes and four fibroblast cell types in endometrium. We showed that the ciliated epithelium is a typical endometrial cell type, which is consistent with the results of one recent study. 58 Although the existence of ciliated cells in human endometrium was observed long times ago, 59 the gene expression patterns and marker genes were not identified previously. We identified numerous of genes that showed higher expression in ciliated cells, such as TPPP3, CAPS, AGR3 and PIFO. AGR3 is a specialized member of the protein disulfide isomerase (PDI) family, and it has been found to play an unexpected role in the regulation of ciliary beat frequency and mucociliary clearance in the airway. 60 In total, we identified 4,038 immune cells, including monocyte, NK1, NK2 and T cells ( Figure S7A ). Immune cells were mainly enriched in proliferation and late secretory stage, NK2 cells were mainly enriched in late secretory stage, and T cells were mainly distributed in proliferative stage ( Figure S7B ). Monocyte cells highly expressed CD74 and LYZ, T cells highly expressed CD3D and IL7R, and NK cells highly expressed NKG7 and GNLY ( Figure S7C ). In addition, the high expression of antigen presentation and phagocytic genes were observed in monocyte, such as HLA−DRA, HLA−DRB1, HLA−DPA1, C1QA, C1QB, and C1QC, indicating that they play an important role in the clearance of apoptotic cells after endometrial shedding ( Figure S7D ). Functional enrichment analysis suggested that genes highly expressed in monocyte were significantly enriched in antigen processing and presentation, cell–cell adhesion and leukocyte migration ( Figure S7E ). Moreover, genes highly expressed in T cells were enriched in T cell differentiation, T cell activation and leukocyte cell−cell adhesion ( Figure S7E ). These results suggest that immune cells play an important role in the immune response, clearance of apoptotic cells and endometrial tissue remodeling after endometrium shedding. In addition, we found that genes highly expressed in cell types showed expression perturbations in cancer. For example, APOD, DCN, LUM and CFD were significantly highly expressed in fibroblast cells ( Figure S8A ). DCN has been demonstrated to play tumor suppressive functions in cancer, 61 , 62 and we found that it was lowly expressed in EC. CAV1 and VWF were highly expressed in endothelial cells but showed lower expression in EC ( Figure S8B ). Epithelial-related genes, such as SCGB2A1 and CLDN4, were highly expressed in cancer ( Figure S8C ). Moreover, we found that CCL3, C1QA, GZMB and IL32 were highly expressed in cancer ( Figure S8D and E ). SCGB2A1 was found to be perturbed in various cancer types, such as breast and ovarian cancers. 63 , 64 We also found the expressions of several genes were associated with survival of cancer patients ( Figure S8F ). These results highlight the critical roles of these genes in endometrial cancer. In summary, deconvolution of immune microenvironment that drives transcriptional programs throughout the menstrual cycle is important for understanding RNA regulatory biology of endometrium. Our dynamic immune microenvironment analyses provide novel insights into future development of RNA-based treatments for endometriosis and endometrial carcinoma. Our findings provide a theoretical foundation for personalized medicine research on endometrial diseases. However, further experimental validation is needed, including studies using mouse models or cell lines, to minimize the impact on human patients, especially those who wish to preserve their fertility. This emphasizes the importance of tailored approaches in improving treatment efficacy and patient outcomes, while also effectively reducing the occurrence of potential medical incidents.