Transcriptional Profiling of SSEA-1+ Endometrial Epithelial Progenitor Cells Highlights Their Role in Endometrial Regeneration, Remodeling, and Homeostasis.

D.K.H. conceived the project. H.A.‐L., A.V., O.V., J.D., C.J.H., and D.K.H. designed the experiments. H.A.‐L., J.D., C.J.H., and A.V. performed the experimental work. H.A.‐L., J.S., D.G., O.V., A.M., and N.T. were involved in the data analysis. H.A.‐L. and D.K.H. wrote the first draft of the manuscript. All authors, H.A.‐L., J.S., D.G., O.V., A.V., C.J.H., J.D., A.M., A.D., N.T., and D.K.H., were involved in the editing and approval of the final manuscript.

The differential functional ability of EEC subtypes is demonstrated by the well‐established region‐specific functional differences in the endometrium (e.g., functionalis, the site for embryo‐implantation, and basalis, the site retaining regenerative capacity), and SSEA‐1 + EECs have demonstrated some functional properties in keeping with SPCs [ 2 , 10 ]. Therefore, we sought to examine the anticipated distinct transcriptome between SSEA‐1 + and the more differentiated SSEA‐1 − EEC population using freshly isolated and sorted cells from eight human endometrial biopsies. Transcriptomic analysis of the SSEA‐1 + versus SSEA‐1 − EEC populations was performed using a two‐color microarray‐based gene expression technology. Principal component analysis (PCA) demonstrated distinct clustering of SSEA‐1 + and SSEA‐1 − EECs (Figure  1A ). A total of 1054 (71 upregulated and 983 downregulated) differentially expressed genes (DEGs) were identified in the SSEA‐1 + EECs compared to the SSEA‐1 − EEC population (Figure  1B and Table  1 ). Hierarchical clustering was observed for SSEA‐1 + and SSEA‐1 − EEC populations based on their DEGs (Figure  1C,D ). (A) Principal component analysis (PCA) plot of SSEA‐1 + EECs (blue) and SSEA‐1 − EECs (red) across all eight endometrial samples showing distinct clustering. (B) Volcano plot of all 22 860 expressed genes. Blue points represent differentially expressed genes (DEGs) (FDR  1 and 983 downregulated DEGs with logFC < −1). Top upregulated and downregulated DEGs have been annotated within the volcano plot. (C) Heat map and hierarchical clustering of the mRNA expression of top 35 upregulated DEGs, showing distinct transcriptional clustering between SSEA‐1 + EECs (left) and SSEA‐1 − EECs (right). (D) Heat map and hierarchical clustering heatmap of the mRNA expression of top 35 downregulated DEGs, showing distinct transcriptional clustering between SSEA‐1 + EECs (left) and SSEA‐1 − EECs (right). List of differentially expressed genes (DEGs) in the SSEA‐1 + versus SSEA‐1 − endometrial epithelial cell populations (FDR  1 and top 100 downregulated genes with logFC < −1. GO pathways were analyzed according to three functional groups: biological processes, molecular functions and cellular components. The most significantly enriched pathways were related to structural organization and extracellular matrix (ECM). The top three biological processes included “external encapsuling structure organization” (76 genes, FDR = 1.08e−27), “ossification” (68 genes, FDR = 6.39e−15) and “regulation of vasculature development” (50 genes, FDR = 5.26e−12) (Figure  2A and Table  2 ). The top three most significantly enriched molecular functions were “extracellular matrix structural constituent” (47 genes, FDR = 4.48e−19), “collagen binding” (27 genes, FDR = 5.27e−16), and “extracellular matrix structural constituent conferring tensile strength” (18 genes, FDR = 1.24e−10) (Figure  2B and Table  3 ). The top three most significantly enriched cellular components were “collagen containing extracellular matrix” (89 genes, FDR = 3.12e−26), “endoplasmic reticulum lumen” (48 genes, FDR = 8.01e−09), and “collagen trimer” (23 genes, FDR = 3.48e−08) (Figure  2C and Table  4 ). Gene enrichment analysis. (A) Bar chart of top 10 enriched Gene Ontology (GO)‐Biological Processes (BP). (B) Enrichment plot of top 10 enriched GO‐Molecular Functions (MF). (C) Bar chart of top 10 enriched GO‐Cellular Components (CC). Top 10 most significantly enriched Gene Ontology (GO)‐Biological Processes (BP) pathways in the SSEA‐1 + versus SSEA‐1 − endometrial epithelial cell populations. FGF18, TNF, EMP2, THBS1, PLXND1, AGT, PDE3B, SERPINF1, HOXA5, COL4A2, ADAM12, DLL1, GATA2, WNT4, HGF, HMGA2, FLT1, CLDN5, GPR4, TIE1, NINJ1, PGF, ECM1, TWIST1, VASH1, CCBE1, NPPB, ANGPTL7, CHI3L1, RECK, TGFB2, ISM1, VEGF, SEMA5A, IL1B, ADM, AKT3, ANGPT2, PRKCA, TEK, KDR, THBS2, DCN, BMPER, ACVRL1, GREM1, GPNMB, SULF1, SPARC, WNT5A FGF18, TGFBI, FOXD1, HOXA5, SNAI1, EXT1, ARID5B, SERPINH1, SOX6, HMGA2, BGN, PPARGC1A, GPR4, RUNX2, HOXA3, EFEMP1, ECM1, HOXD3, WT1, ACTA2, GLI3, SOX5, RARB, STC1, MAF, NAMPT, CHI3L1, MGP, LOX, WNT5B, CHST11, HAND2, COL5A1, ANXA6, SNAI2, BMP2, COL3A1, HOXA11, ZEB1, COL1A1, CTSK, PRRX1, ACVRL1, GREM1, SULF1, LOXL2, PDGFRB, WNT5A MMP7, MMP11, MMP19, SERPINH1, WNT4, ADAMTS3, LARP6, MMP14, PRSS2, MMP16, RCN3, MMP9, MMP10, COL5A1, MMP12, F2R, COL1A2, MFAP4, COL1A1, CTSK, VIM, MMP1, MRC2, MMP2, FAP, PDGFRB, MMP3 AREG, TNF, SPINT, FOXD1, PLXND1, AGT, HOXA5, DCHS1, VDR, EXT1, WNT4, COL4A1, HGF, FGF7, TIE1, EDNRA, SALL1, WT1, GLI3, SMO, RSPO3, MMP14, WNT2, PKD2, NFATC4, PRDM1, LAMA1, PBX1, SNAI2, BMP2, ADM, LEF1, FOXF1, HOXA11, AR, KDR, TBX3, COL13A1, TBX2, GREM1, SULF1, WNT5A FOXF2, FOXD1, KITLG, TGFB1I1, HOXA5, HEYL, SNAI1, DCHS1, CITED2, ZFPM2, EXT1, WNT4, HGF, HMGA2, EDNRA, TWIST1, RBM24, WT1, ACTA2, ROBO1, SMO, FGFR1, HAS2, WNT2, ACTG2, PKD2, VASN, RADIL, TGFB2, HAND2, ANXA6, SNAI2, BMP2, SEMA5A, IL1B, CDH2, SOX11, LEF1, FOXF1, HTR2B, PDPN, TBX3, DAB2, DACT3, COL1A1, TBX2, ACVRL1, GREM1, LOXL2, PDGFRB, WNT5A MMP11, EXT1, SERPINH1, ADAMTS3, COL5A3, EFEMP2, LOX, FMOD, TGFB2, COL5A1, FKBP10, COL1A2, DDR2, COL3A1, AEBP1, COL12A1, COL13A1, COL1A1, GREM1, COL5A2, LUM, LOXL2 KRT16, FGF18, TNF, RAB25, RIN2, EMP2, THBS1, PTPRG, PLXND1, KITLG, AGT, SERPINF1, GJA1, ARID5B, GATA2, FGF7, FSTL1, DAAM2, KANK2, EDNRA, TWIST1, NR2F2, VASH1, ACTA2, PRSS3, ROBO1, SMO, CCBE1, PTN, FGFR1, STC1, HAS2, NANOS1, SYDE1, CDH13, MMP9, ZEB2, RADIL, TGFB2, HAND2, ANXA6, VEGFC, SEMA5A, DDR2, CDH2, APOE, SRPX2, AKT3, ANGPT2, HTR2B, PRKCA, TEK, KDR, ITGA4, DCN, BMPER, ACVRL1, GREM1, FAP, LOXL2, SPARC, WNT5A FGF18, AREG, TNF, THBS1, MTSS1, LAMB1, SERPINF1, HOXA5, VDR, GATA2, HGF, FGF7, LIMS2, FLT1, CEBPB, TIE1, PGF, ECM1, TWIST1, KLF9, NR2F2, GLI1, VASH1, ROBO1, NLRC3, SMO, PTN, CCL26, FGFR1, MMP14, HAS2, WNT2, CDH13, TGFB2, SNAI2, MMP12, BMP2, VEGFC, SEMA5A, COL8A1, APOE, SOX11, IGF2, AKT3, HTR2B, PRKCA, TEK, AR, KDR, ITGA4, DAB2, HTRA1, BMPER, TBX2, ACVRL1, SULF1, FAP, LOXL2, SPARC, WNT5A Top 10 most significantly enriched Gene Ontology (GO)‐Molecular Functions (MF) pathways in the SSEA‐1 + versus SSEA‐1 − endometrial epithelial cell populations. POSTN, TGFBI, VCAN, THBS1, LAMB1, COL4A2, CTHRC1, COL24A1, COL4A1, BGN, EFEMP1, ECM1, COL23A1, COL5A3, COL6A2, EDIL3, COL15A1, VWF, CHI3L1, EFEMP2, MGP, FBLN2, COL6A1, FMOD, LAMA1, LTBP2, COL6A3, COL5A1, PCOLCE, COL8A1, COL1A2, COL3A1, SRPX2, MFAP4, AEBP1, COL12A1, THBS2, DCN, FBLN5, COL13A1, COL7A1, COL1A1, SRPX, COL5A2, LUM, SPARC, NID1 TGFBI, THBS1, SERPINH1, C1QTNF1, ITGA10, COL5A3, COL6A2, CCBE1, LRRC15, VWF, ITGA11, MMP9, LOX, COL6A1, PCOLCE, MMP12, DDR2, SPARCL1, AEBP1, ITGA1, ANTXR1, DCN, CTSK, MRC2, LUM, SPARC, NID1 COL4A2, COL24A1, COL4A1, COL23A1, COL5A3, COL6A2, COL15A1, COL6A1, COL6A3, COL5A1, COL8A1, COL1A2, COL3A1, COL12A1, COL13A1, COL7A1, COL1A1, COL5A2 TGFBI, EMP2, THBS1, ITGA7, LAMB1, FCER2, ITGA10, ITGA8, ESM1, EGFL6, EDIL3, PTN, MMP14, VWF, ITGA11, JAM2, COL5A1, IL1B, COL3A1, IGF2, PRKCA, KDR, ITGA1, ADAMTS5, ITGA4, FBLN5, THY1, GPNMB, S1PR3, FAP POSTN, VCAN, THBS1, SLIT3, PTPRS, TNFAIP6, FGF7, FSTL1, BGN, PGF, ADAMTS3, COL23A1, CRISPLD2, COL5A3, RSPO3, PTN, LAYN, LPL, FGFR1, SERPINE2, EFEMP2, LXN, LTBP2, GREM2, COL5A1, PCOLCE, ANXA6, SEMA5A, APOE, ADAMTS5, THBS2, DCN, COL13A1, GPNMB, SULF1 COL4A1, COL6A1, COL5A1, COL1A2, COL3A1, COL1A1, PDGFRB, PDGFRA Top 10 most significantly enriched Gene Ontology (GO)‐Cellular Components (CC) pathways in the SSEA‐1 + versus SSEA‐1 − endometrial epithelial cell populations. Hallmark and KEGG pathways demonstrated that SSEA‐1 + EECs are important in pathways related to epithelial‐to‐mesenchymal transition and angiogenesis. The top three most significantly enriched Hallmark pathways relevant to regeneration and remodeling included “epithelial mesenchymal transition” (81 genes, FDR = 2.8e−41), “coagulation” (26 genes, FDR = 0.0001), and “angiogenesis” (12 genes, FDR = 0.0001) (Figure  3A and Table  5 ). The top three most significantly enriched KEGG pathways included “focal adhesion” (40 genes, FDR = 3.16e‐11), “ECM receptor interaction” (25 genes, FDR = 8.94e−11), and “arrhythmogenic right ventricular cardiomyopathy ARVC” (15 genes, FDR = 0.00047) (Figure  3B and Table  6 ). Gene enrichment analysis (A) Bar chart of top 10 enriched Hallmark pathways (B) Bar chart of top 10 enriched KEGG pathways. Top 10 most significantly enriched Hallmark pathways in the SSEA‐1 + versus SSEA‐1 − endometrial epithelial cell populations. Top 10 most significantly enriched KEGG pathways in the SSEA‐1 + versus SSEA‐1 − endometrial epithelial cell populations. Using IPA we were further able to biologically interpret the DEGs using the extensive knowledge base contained within the platform. Firstly, IPA generated associations between the DEGs and specific “diseases and disorders”, “molecular and cellular functions”, and “physiological system development and function”. The top five within each of these categories are highlighted in Table  7 . These suggest that SSEA‐1 + EECs are important for tissue homeostasis including functions such as tissue proliferation, development, maintenance, survival, and angiogenesis. Top diseases and biofunctions in the SSEA‐1 + versus SSEA‐1 − endometrial epithelial cell populations. Further exploration of “diseases and functions” associated with the DEGs, predicted “organismal death” ( z ‐score = 14.60, FDR = 4.85e−32) and “dysgenesis” ( z ‐score = 7.42, FDR = 8.41e−13) to be within the top activated biofunctions with a z ‐score > 2 and FDR < 0.05. The bottom z ‐scores (with z ‐scores < −2 and FDR < 0.05) predicted inhibition of “size of body” ( z ‐score = −9.28, FDR = 9.09e−15), “cell movement” ( z ‐score = −8.17, FDR = 1.13e−52), “cell movement of tumor cell lines” ( z ‐score = −7.97, FDR = 1.15e−32), “migration of cells” ( z ‐score = −7.81, FDR = 2.54e−49) and “migration of tumor cell lines” ( z ‐score = −7.67, FDR = 1.44e−28). The IPA downstream effects analysis table (Table  S6 ) shows the tabular output that predicts “organismal death” as the predicted top activated biofunctions, based on the observed changes in gene expression within the SSEA‐1 + EECs. Using the QIAGEN IPA knowledge base, well‐characterized metabolic and cell signaling pathways were explored based on overlap from the DEGs. IPA's canonical pathway analysis suggests SSEA‐1 + EECs to have a tumor suppressor and homeostatic function and a reduced hormone responsive phenotype. Significant overlap was found between the genes involved within five canonical pathways which were predicted to be activated within the SSEA‐1 + EECs. These included “RhoGDI signaling” ( p ‐value = 6.52e−04, z ‐score = 3.21), “endocannabinoid cancer inhibition pathway” ( p ‐value = 6.79e−03, z ‐score = 3.21), “PTEN signaling” ( p ‐value = 1.81e−02, z ‐score = 1.74), “chaperone mediated autophagy signalling pathway” ( p ‐value = 0.04, z ‐score = 2.83), and “CDX gastrointestinal cancer signalling pathway” ( p ‐value = 8.69e−05, z ‐score 2.29). An example of the “PTEN signalling pathway” and “endocannabinoid cancer inhibition pathway” are shown in Figure  S1A,B , highlighting the tumor suppressor pathways that are predicted to be activated within the SSEA‐1 + EECs. Overall, an overwhelming majority of 83 significantly overlapping canonical pathways were predicted to be inhibited in the SSEA‐1 + EECs. The top five predicted to be most inhibited within the SSEA‐1 + EECs were “molecular mechanisms of cancer” ( p ‐value = 5.55e−06, z ‐score = −7.54), “FAK signaling” ( p ‐value = 1.44e−08, z‐score = −7.39), “pulmonary fibrosis signalling pathway” ( p ‐value = 8.15e−22, z ‐score = −7.18), “phagosome formation” ( p ‐value = 6.77e−04, z ‐score = −6.93) and “CREB signalling in neurons” ( p ‐value = 2.5e−04, z ‐score = −6.87). The canonical “estrogen receptor signalling” pathway is also seen to significantly overlap with the DEGs ( p ‐value = 4.79e−03); however, using the IPA Knowledge Base, this pathway is predicted to be inhibited within the SSEA‐1 + EEC population (z‐score = −4.16). This pathway is seen to control functions such as “tumor cell proliferation”, “tumor EMT”, “metastasis”, “cell proliferation” and “migration of tumour cells”, which are all predicted to be inhibited due to the downregulation of DEGs, mainly matrix metallopeptidases (MMPs) involved in this pathway. The pathway “estrogen receptor signalling” and genes involved are detailed in Figure  S2 .

“Upstream analysis” within IPA was used to predict significantly activated and inhibited upstream regulators within the SSEA‐1 + EEC population and their associated mechanistic and causal networks, which may explain the observed expression changes. Upstream regulator analysis allows determination of likely upstream regulators that are connected to the DEGs either through direct or indirect relationships. To understand the relationship between upstream regulators and DEGs within the dataset, mechanistic and causal networks were used. As upstream regulators are not necessarily independent of one another, mechanistic networks allowed us to hypothesize other connecting regulators that are likely to be part of the same signaling or causal mechanism, whereas causal networks allowed the detection of novel master upstream regulators [ 21 ]. Within our dataset, the transcription factor SPDEF ( p ‐value of overlap = 1.05e−10, z ‐score = 4.30) was predicted to be one of the top‐most activated upstream regulators, with significant overlap of their target molecules seen within this microarray DEGs dataset. On the contrary, TGFB1 is predicted to be one of the top‐most inhibited upstream regulators with the most significant number of target molecules of overlap within our dataset ( p ‐value of overlap = 8.64e−40, z ‐score = −8.20). The mechanistic networks of SPDEF and TGFB1 and their downstream target molecules within our DEGs are displayed within Figure  4A,B . When exploring causal networks for a master upstream regulator, TGFB1 was also found to be predicted as the top‐most inhibited master regulator for the DEGs ( p ‐value of overlap = 8.64e−40, z ‐score = −8.20) (Figure  4C ). In addition, upstream regulator analysis revealed overlap between several upstream regulators predicted to control the biological mechanisms involved in SSEA‐1 + EECs and regulators known to be involved in endometrial carcinogenesis [ 26 , 27 , 28 ]. These regulators include ARID2, KRAS, CTNNB1, PIK3CA , and TP63 . IPA's upstream regulator analysis. Mechanistic networks for (A) SPDEF and (B) TGFB1 displaying other predicted connecting upstream regulators likely to be part of the same signaling mechanism. (C) Causal network for TGFB1 predicted to be a master regulator and its downstream targets. Orange nodes represents predicted activation. Green nodes represent significantly downregulated DEGs within the dataset. Solid lines represent a direct relationship and dashed lines represent indirect relationships between nodes. Orange lines predict a relationship resulting in activation and blue lines predict a relationship resulting in inhibition. Microarray data was validated using RT‐qPCR to assess the mRNA expression levels of a selection of significantly up‐and downregulated genes. MMP7 ( p  = 0.016), MMP26 ( p  = 0.016), SPRR2A ( p  = 0.016) and FUT3 ( p  = 0.03) showed significantly higher mRNA expression levels within the SSEA‐1 + EECs compared with the SSEA‐1 − EECs (Figure  5A ). In contrast, significantly lower gene expression levels for TMEM158 , AXIN2 , ST3GAL2 , WNT5A , TCF4 , and ZEB1 were confirmed within the SSEA‐1 + EECs compared to the SSEA‐1 − EECs ( p  = 0.016; Figure  5B ). Dual‐IF staining of human endometrial tissue sections confirmed coexpression of SSEA‐1 with proteins encoded by significantly upregulated genes identified by microarray analysis, including MMP7, MMP26, C11orf52, and CD47 (Figure  6 ). Quantitative real‐time PCR analysis of transcripts validating microarray gene expression data for (A) upregulated genes: MMP7 , MMP26 , SPRR2A , FUT3 and (B) downregulated genes: TMEM158 , AXIN2 , ST3GAL2 , WNT5A , TCF4 , ZEB1 . Relative normalized mRNA expression was confirmed to be significantly greater in SSEA‐1 + EECs for MMP7 ( p  = 0.016), MMP26 ( p  = 0.016), SPRR2A ( p  = 0.016) and FUT3 ( p  = 0.031). Relative normalized mRNA expression was confirmed to be significantly lower in SSEA‐1 + EECs for TMEM158 ( p  = 0.016), AXIN2 ( p  = 0.016), ST3GAL2 ( p  = 0.016), WNT5A ( p  = 0.016), TCF4 ( p  = 0.016) and ZEB1 ( p  = 0.016). significant differences p  < 0.05 are depicted with “*”. Representative micrographs show co‐localization of (A) MMP7 , (B) MMP26 , (C) C11orf52 and (D) CD47 with SSEA‐1 within the basal glands of the endometrium using dual immunofluorescence. Scale bars = 50 μm. DEGs MMP7, MMP26, SPRR2A , and FUT3 were found to be significantly upregulated within SSEA‐1 + EECs by both microarray gene expression analysis and RT‐qPCR. They were mapped to the publicly available single‐cell transcriptomic dataset published by Garcia et al. [ 22 ] as a further method of external validation to assess the epithelial cell types that express these genes (Figure  7A,B ). At an expression level > 2, MMP7 is predominantly expressed by SOX9 + ( SOX9 + EECs: 418 cells; SOX9 + proliferative EECs: 382 cells) and SOX9 + / LGR5 + EECs mainly within the proliferative phase ( SOX9 + / LRG5 + EECs: 398 cells) and within a much smaller number of ciliated LGR5 + EECs (ciliated EECs: 10 cells; ciliated LGR5 + EECs: 16 cells) (Figure  7C ). At an expression level > 2, MMP26 was found to be predominantly expressed by glandular SOX9 + EECs within the proliferative phase (glandular EECs: 270 cells; SOX9 + EECs: 550 cells; SOX9 + proliferative EECs: 125 cells) (Figure  7D ). When combined, SOX9 + cells within the proliferative phase were found to predominantly express both MMP7 and MMP26 at an expression level > 2 ( SOX9 + EECs: 74 cells; SOX9 + proliferative EECs: 47 cells; SOX9 + /LRG5 + EECs: 8 cells), suggesting co‐expression of SSEA‐1 + and SOX9 + EECs, in keeping with previous literature [ 10 ] (Figure  7E ). UMAP of single‐cell EEC transcriptomic data published by Garcia et al. [ 22 ]. (A) UMAP of all subclustered and subsampled epithelial populations. (B) UMAP mapping MMP7 , MMP26 , SPRR2A and FUT3 to single‐cell EEC clusters. (C) Violin plots of MMP7 and (D) MMP26 expressing single‐cell EECs and number of epithelial cell types expressing each marker at an expression level > 2. (E) UMAP of cells co‐expressing MMP7 and MMP26 at an expression level > 2, mapped to single‐cell EEC clusters. The ovarian hormone milieu in the secretory phase of the menstrual cycle induces endometrial cell differentiation [ 29 ]. To simulate this in vitro, human endometrial epithelial organoids were exposed to hormonal stimuli. The expression of proposed EEC SPC markers SSEA‐1, SOX9 and N‐cadherin, PR as a marker of estrogen response, and the proliferation marker Ki‐67 were assessed (Figure  8A,B ). Compared with controls (standard media), hormone treatment resulted in a higher mean number of organoids expressing SSEA‐1 (34.3 out of 73.3 (46.8%) versus 55 out of 82 (67.1%)), PR (4 out of 73.7 (5.4%) versus 37 out of 62.3 (59.4%)), and Ki‐67 (58.7 out of 74 (79.3%) versus 73.3 out of 80.7 (90.8%)), while a lower proportion expressed N‐cadherin (4.7 out of 70.7 (6.6%) versus 2.7 out of 79.3 (3.4%)). SOX9 was ubiquitously present across all organoids exposed to both culture conditions (69 out of 69 (100%) versus 98.33 out of 98.33 (100%)). Immunohistochemical staining for PR, Ki‐67, SSEA‐1, SOX9 and N‐Cadherin in endometrial epithelial organoids exposed to either standard media or standard media plus hormones. (A) Representative micrographs of immunohistochemical (IHC) staining for each protein. Scale bar is 200 μm. (B) IHC scoring displayed as the mean number of positively or negatively stainedg endometrial epithelial organoids, Error bars represent the standard error of the mean (SEM).

This study reports the transcriptional profile of the first postulated human endometrial epithelial SPC population, SSEA‐1 + EECs, which agrees with their expected SPC function. Although it has been hypothesized that the basalis EECs have a role in endometrial regeneration, our data show for the first time that the transcriptional signature of SSEA‐1 + EECs is in keeping with the expected SPC phenotypic and functional pathways. Relative to other EECs, SSEA‐1 + EECs possess a distinct transcriptional signature that highlights their importance in regeneration, remodeling, and vascularization of the human endometrium. We have employed multiple approaches in pathway and gene set analysis (GO, KEGG, Hallmark and IPA canonical pathways) which cover a wide range of overlapping comprehensive data repositories and produce robust, corroborative results suggesting that the SSEA‐1 + EEC transcriptional profile is consistent with SPC function. These findings are pivotal to guide future research to understand and target the specific pathways for diagnosis and therapy for gynecological conditions. For example, the SPC function, responsible for endometrial regeneration, is postulated to be perturbed in persistent endometrial pathologies; thus, examinations of the SSEA‐1 + EEC transcriptome in diseased endometrium will provide clinically impactful results. Pathway analysis revealed important processes related to “vascular and mesenchymal development”, “extracellular matrix structural constituent”, “angiogenesis”, “TNFα and Notch signaling” related to SSEA‐1 + EECs. These were associated with biological functions, such as apoptosis, the inhibition of cell movement, cell proliferation and cell growth, supportive of the quiescent SPC phenotype postulated of the SSEA‐1 + EECs. Corresponding with the findings of Valentijn et al. [ 2 ], significantly reduced PR (a well‐established downstream target of estrogen action) expression levels and a downregulation of estrogen receptor pathway were also confirmed. Abnormal endometrial SPC populations have been implicated in gynecological pathologies such as endometriosis, endometrial cancer, and recurrent reproductive failure [ 8 , 10 , 30 ]. Understanding the biological processes and functional pathways specific to SSEA‐1 + EECs in the normal endometrium is essential for improving our understanding of endometrial physiology as well as in the diseased state. In this respect, the identification of important biofunctions in SSEA‐1 + EECs, such as the inhibition of tumor cell invasion and migration, in addition to maintaining endometrial homeostasis, is notable. The activation of functional pathways such as “PTEN signaling” and “endocannabinoid cancer inhibition” in IPA's “diseases and functions” and “canonical pathways” analysis demonstrates the fine equilibrium that is expected by a progenitor cell population. These would support rapid and robust proliferation at the time of regeneration, while providing protection from hyperplastic or carcinogenic transformation [ 31 , 32 , 33 , 34 , 35 ]. Inhibition of the “estrogen receptor signalling pathway” within SSEA‐1 + EECs also confirms their restrained hormone responsiveness, a feature that is important to prevent abnormal endometrial growth and carcinogenesis [ 36 ]. Our microarray data was validated using multiple methods. Ten significantly up‐or downregulated genes from the microarray analysis were internally validated using RT‐qPCR on the same sample set, demonstrating concordance. Additionally, co‐localization of the gene products of four selected genes with SSEA‐1 + EECs was demonstrated in an external biological cohort of endometrial samples. A further external in silico validation was conducted using a publicly available single‐cell transcriptomics dataset, which maps the epithelial subtypes within the human endometrium [ 22 ]. The significantly upregulated DEGs ( MMP7, MMP26, and SPRR2A ) of SSEA‐1 + EECs that we had already validated with RT‐qPCR were mapped in silico with FUT3 . Since there is no direct gene coding for the glycoprotein SSEA‐1, FUT3 may be an important surrogate marker for SSEA‐1 + EECs [ 37 , 38 ]. While SPRR2A and FUT3 expression were sparse and at low levels within the Garcia‐Alonso et al. dataset [ 22 ], MMP7 and MMP26 showed much stronger expression mainly localized to SOX9 + and SOX9 + /LGR5 + EECs. This also corresponds to the previously described progenitor cell phenotype of SSEA‐1 + /SOX9 + co‐localization within basal glandular epithelial cells by Valentijn et al. [ 2 ]. MMP7, MMP26, and SPRR2A are seen to be important for wound and mucosal repair, and therefore may contribute to epithelial barrier restoration following menstruation [ 39 , 40 ]. In 2012, Nguyen et al. used a gene microarray approach to assess the transcriptome of highly purified EpCAM + cells from post‐menopausal endometrium (proposed to represent the pre‐menopausal basalis only) versus epithelial cells from pre‐menopausal full‐thickness endometrium (which contained functionalis and basalis cells). Our studies are not comparable methodologically, and thus, expectedly, there was no overall correlation found between the two datasets (spearman correlation coefficient = −0.04). However, we did find an overlap between seven upregulated and 33 downregulated DEGs across the two microarray datasets (Figure  S3 ). Of the genes validated by RT‐qPCR, the significantly upregulated genes MMP7 and FUT3 and significantly downregulated genes TMEM158 and ZEB1 correlated with the post‐menopausal samples included by Nguyen et al. [ 12 ]. Interestingly, while Nguyen et al. found nuclear AXIN2 to be upregulated in post‐menopausal endometrium compared to pre‐menopausal endometrium, we found AXIN2 , another proposed endometrial progenitor marker, to be downregulated in the pre‐menopausal SSEA‐1 + EECs. Further exploration of EECs of pre‐and postmenopausal endometrium is required in the future to complete our understanding of the discerning features relevant to SPC populations exposed to these two different hormonal milieus. The effect of hormones in influencing SSEA‐1 + EECs were examined using an endometrial epithelial organoid model. Exposure to a combination of oestradiol and progestin demonstrated stronger immunostaining of PR and Ki‐67, thus, estrogen‐responsive proliferation. Although the inhibited “estrogen receptor signalling pathway” within SSEA‐1 + EECs in our in silico study suggest their restrained hormone responsiveness, in a monocellular organoid system that we employed, the SSEA‐1 + cells appear to respond to oestrogenic signals. However, the hormonal regulation of individual cell types in the intact endometrium is more complex, and the epithelial response to progesterone for example, is exerted via its action on the stromal cells [ 41 ]. Our organoid system did not contain stromal component, therefore was not suitable to fully assess the hormone regulation of SSEA1 + EECs that occur in vivo. Differential expression analysis showed significantly higher expression of VIM , PDGFRA , and PDGFRB within the SSEA‐1 − EECs. This is likely to have been secondary to stromal contamination obtained from the flow‐through within the MACS cell separation process and is a recognized limitation of this technique. We have previously published the typical cell purity from MACs to be > 75% [ 2 ]. Our work has provided critical information regarding the transcriptional profile of premenopausal SSEA‐1 + EECs, supporting their pivotal role as progenitor cells driving endometrial regeneration and remodeling, which may play a role in protecting and regulating healthy endometrial homeostasis. Although SSEA‐1 + EECs generally reside in the basalis endometrial epithelium of healthy women, SSEA‐1 + EECs are also seen in the luminal epithelium [ 30 ] and in the functionalis of some women with gynecological pathologies such as endometriosis [ 10 ]. Further studies are now warranted using novel techniques such as spatial transcriptomics to assess in situ transcriptomics of different SSEA‐1 + EEC subpopulations. This would allow us to understand the functional differences of SSEA‐1 + EEC subtypes based on their region‐specific location within the endometrium. Further functional studies are required to explore SSEA‐1 + EEC mechanistic pathways in model systems that contain all main cell types recapitulating the in vivo endometrium and by manipulation of their upstream regulators.

The human endometrium is an indispensable organ in the female reproductive system, with an unparalleled regenerative ability. It is organized into two functionally and phenotypically different layers: (1) a superficial functionalis, which appears only during the reproductive years of a woman's life, to undergo an iterative monthly cycle of shedding, regeneration, remodeling, and transformation according to the ovarian steroid cue, and (2) a permanent deeper basalis, which persists throughout a woman's life. Endometrial stem/progenitor cells (SPCs) are postulated to reside in the deeper basalis and hypothesized to be responsible for the extraordinary capacity for full phenotypic and functional restoration of the endometrial functionalis layer [ 1 , 2 ]. This efficient and scarless regeneration of the endometrial functionalis happens after menstrual shedding, parturition, and in the revived post‐menopausal endometrium, when stimulated with exogenous hormones [ 2 , 3 , 4 , 5 , 6 , 7 ]. Several types of SPCs have been identified within the endometrium, including endometrial mesenchymal stem cells (eMSCs), endometrial epithelial progenitor cells (eEPCs), and side population cells [ 8 ]. In 2013, we described an endometrial epithelial cell (EEC) sub‐population that express surface marker stage‐specific embryonic antigen‐1 (SSEA‐1). SSEA‐1 + EECs assume the postulated SPC location, the basalis, while demonstrating some expected SPC functional properties. SSEA‐1 + EECs exhibit higher telomerase activity, longer telomere lengths, quiescence, and resistance to hormonal stimuli (lower estrogen receptor α ( ESR1 ) and progesterone receptor ( PR )), suggesting a less differentiated phenotype compared to SSEA‐1 − EECs [ 2 , 8 , 9 ]. In an in vitro three‐dimensional (3D) culture system, SSEA‐1 + EECs demonstrated a higher propensity to generate endometrial gland‐like spheroids and growth similar to that observed after endometrial denudation [ 2 ]. SSEA‐1 + EECs also express nuclear SRY‐box transcription factor 9 (nSOX9) and nuclear β‐catenin, suggesting activated Wnt signaling [ 2 ]. The involvement of SSEA‐1 + /nSOX9 + EECs in endometrial pathologies has been implied by their aberrant location in women with endometriosis [ 10 ]. Subsequent studies have proposed other EEC progenitor markers including leucine‐rich repeat containing G protein‐coupled receptor 5 ( LGR5 ) [ 11 ], axis inhibition protein 2 ( AXIN2 ) [ 12 ] and cadherin 2 ( CDH2 ) [ 13 ], with a suggested epithelial SPC hierarchy based on the exact location within the basalis glands [ 13 ]. Since the EEC SPC population makes a vital contribution to routine endometrial regeneration and pregnancy establishment, understanding the transcriptional profile and the associated functional pathways of these cells is essential to understanding their role within normal endometrial physiology [ 8 , 9 ]. Furthermore, such data will be invaluable for the identification of pathology‐related aberrations in SPC populations. In this study, we characterize the transcriptional profile of isolated SSEA‐1 enriched (SSEA‐1 + ) EECs compared to SSEA‐1 depleted (SSEA‐1 − ) EECs and identify differences in key functional pathways. Functional work employing an in vitro endometrial organoid model demonstrated hormonal regulation in endometrial regeneration and differentiation and the role SSEA‐1 + EECs play in this process.

University of Liverpool has received honoraria for consultancy from Theramex and has received payment for presentations from Theramex and Gideon Richter for D.K.H.'s work. The remaining authors have no Conflicts of interests to report.

Raw microarray data is available from the Gene Expression Omnibus (GEO at the National Center for Biotechnology Information, http://www.ncbi.nlm.nih.gov/geo/ ; accession number GSE280323). Endometrial samples were collected with informed, written consent from women attending the Liverpool Women's Hospital for benign gynecological conditions. For microarray and RT‐qPCR experiments, endometrial tissue samples were obtained from eight patients. Three further endometrial samples were obtained for epithelial organoid culture. All patients had regular menstrual cycles and did not take any hormonal therapy in the 3 months prior to enrollment. Age of participants ranged between 30 and 47 years. Further clinical characteristics of the participants are detailed in Table  S1 . Collection of human endometrial samples was approved by the Adult Ethics committees (REC references; 09/HI005/55, 19/WA/0271 and 19/SC/0449). From those attending for a hysterectomy, a full thickness wedge endometrial biopsy was taken, while a pipelle endometrial sample was obtained from those attending for a laparoscopy. Freshly harvested endometrial tissue was mechanically and enzymatically digested and filtered into stromal and epithelial cell fractions [ 2 , 10 ]. The epithelial single‐cell suspension was labeled with anti‐SSEA1 (CD15) MicroBeads (#130‐094‐530, Miltenyi Biotec, UK) and magnetic bead (MACS) sorted into SSEA‐1 + and SSEA‐1 − populations using MACS separation columns (MS columns, Miltenyi Biotec, UK) according to the manufacturer's instructions. Cell purity using this technique has previously been reported [ 2 ]. Total RNA from SSEA‐1 + and SSEA‐1 − fractions was extracted using the TRIzol Plus RNA Purification System (Invitrogen Life Technologies, Paisley, UK) and quantified using a NanoDrop 1000 (ThermoFisher Scientific, UK) and was DNAse treated (Promega, UK) to remove genomic DNA according to the manufacturer's instructions. A portion of DNAse treated RNA was retained for microarray analysis from each sample before the remaining RNA was used for cDNA synthesis. 1 μg of DNase treated RNA was reverse transcribed using iScript cDNA Synthesis Kit (Bio‐Rad Laboratories Ltd., UK) as per the manufacturer's instructions. Microarray analysis was performed at the Centre for Genomic Research (University of Liverpool). DNase treated RNA was quantified and the quality assessed using the Agilent 2100 Bioanalyser RNA 6000 Nano and Pico chips (Agilent Technologies, Santa Clara, US). Samples were prepared for hybridization onto Agilent SurePrint G3 format 8 × 60k Human Genome Array (Agilent Technologies). 40 ng RNA was used in the initial input and the Agilent Low Input Quick Amp Labelling Kit (Agilent Technologies) was used for target preparation, using the Two‐Color Microarray‐Based Gene Expression Analysis protocol, according to the manufacturer's instructions. Labeled amplified cRNA was purified using the Qiagen RNeasy Minikit (Qiagen, Hilden, Germany), according to the manufacturer's instructions. The cRNA was quantified by Nanodrop and the yield and specific activity calculated. A total of 600 ng (300 ng of both Cy3 and Cy5) of labeled cRNA was fragmented and, following the addition of 2 × HiRPM hybridization buffer, samples were loaded onto arrays and hybridized for 17 h at 65°C in an Agilent hybridization oven (Agilent Technologies). Following hybridization, the arrays were washed using the Agilent gene expression wash buffer kit and scanned using the Agilent G2505C scanner with Agilent G3_GX_2colour settings (Agilent Technologies). Data was extracted using the Agilent feature extraction software version 11.0.1.1. Data processing and statistical analysis of microarray data was performed using R (Version 4.2.3) (R Core Team, 2024. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing. Available at: https://www.r‐project.org/ ). Data was background corrected using the “normexp” function [ 14 ], normalized within arrays using “vsn” [ 15 ] and “loess” [ 16 ] methods, and normalized between arrays using the “Aquantile” method [ 17 ] in limma. Log2 normalized expression data was centred and visualized with principal component analysis. Microarray data was analyzed in limma by applying a linear model using the empirical bayes moderated t test, and p ‐values were adjusted to correct for multiple tests using false discovery rate (FDR) using the Benjamini and Hochberg (BH) method [ 18 ]. Differentially expressed genes (DEGs) were filtered using FDR  1 or < −1. Enriched Hallmark, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed using the Over Representation Analysis method [ 19 ] via the ClusterProfiler “enricher” function [ 20 ] on the DEGs using all genes inputted into the differential expression analysis as the background. Multiple testing correction was performed using BH, and a FDR cutoff of 0.05 was applied. Canonical pathway analysis, upstream regulators, mechanistic and causal network analysis, and associations with diseases and functions were generated using QIAGEN Ingenuity Pathway Analysis (IPA) (QIAGEN Inc. https://digitalinsights.qiagen.com/IPA ) [ 21 ] by inputting the DEGs into the “core analysis” function and a cut‐off set to logFC  1. The significance of biofunctions and canonical pathways was tested by the Fisher's Exact Test p ‐values and analysis cutoffs set to p ‐value  2 (to assess activation) or < −2 (to assess inhibition). External validation of microarray data against publicly available endometrial single cell spatial transcriptomic data published by Garcia‐Alonso et al. [ 22 ] was performed using R. Garcia‐Alonso et al. pre‐processed single cell data was visualized using the previously reported UMAP coordinates and cell type annotations [ 23 ]. To assess epithelial cell type distribution for expression of matrix metallopeptidase 7 (MMP7) and matrix metallopeptidase 26 (MMP26) a normalized log2 expression level > 2 was set. 1 μL cDNA was amplified in triplicate for 40 cycles in a final reaction volume of 10 μL using iTaq Universal SYBR Green Supermix (Bio‐Rad Laboratories Ltd., Hemel Hempstead, Hertfordshire, UK) and a Biorad CFX Connect Real‐Time System (Bio‐Rad Laboratories Ltd.). No template and no reverse transcriptase controls were included for each target in each experiment. Relative gene expression was calculated and normalized to the reference genes glyceraldehyde 3‐phosphate dehydrogenase ( GAPDH ), tyrosine 3‐monooxygenase/tryptophan 5‐monooxygenase activation protein zeta ( YWHAZ ) and actin beta ( ACTB ) and normalized to Ishikawa endometrial cancer cell line (ISK) as an internal control using Biorad CFX manager (version 3.1, Bio‐Rad Laboratories Ltd). The amplification products were verified using gel electrophoresis. The PCR primers are detailed in Table  S2 . Dual immunofluorescence (IF) was performed on 3 μm formalin‐fixed paraffin‐embedded (FFPE) tissue sections. Primary anti‐rabbit MMP7, MMP26, cluster differentiation 47 (CD47) or chromosome 11 open reading frame 52 (C11orf52) antibodies were applied to sections alongside Alexa Fluor 488 conjugated SSEA1 and incubated at 4°C overnight. A combination of anti‐mouse IgG and anti‐rabbit IgG was used as a negative control. The secondary antibody anti‐rabbit Alexa Fluor 555 (Cell Signaling Technology, Hitchin, UK) was used at the recommended concentration. Details of antigen‐retrieval conditions and antibody concentrations are listed in Table  S3 . To diminish unwanted background autofluorescence, slides were incubated at room temperature with TrueVIEW autofluorescence quenching kit (Vector Laboratories, Peterborough, UK) as per the manufacturers instructions, before mounting in Vectashield with DAPI (Vector Laboratories). IF‐stained slides were visualized on a Zeiss LSM 800 confocal microscope fitted with GaAsP detectors. Laser lines used were of 405, 488, and 561 nm excitation wavelengths. Zen Blue software was used for image capture and Image J [ 24 ] for processing. Endometrial tissue was minced into small pieces (< 1 mm) using a scalpel blade and digested with 1 mg/mL Dispase II (Gibco), 2 mg/mL collagenase type I (Gibco) and 80  μg/ml deoxyribonuclease (DNase) 1 (Merck) for ~1 h at 37°C in a shaking water bath. Digests were periodically triturated to enhance tissue breakdown and observed under a microscope to check for the presence of free epithelial glands. The digests were passed through a 40 μm cell sieve (Falcon) to separate glandular (retentate) and stromal (flow‐through) elements. Glandular elements were collected by backwashing the sieve into a collection dish. Endometrial epithelial organoids were generated from the glandular fraction following an established protocol [ 25 ]. Briefly, glands were centrifuged at 500x  g for 5 min and resuspended in phenol‐free DMEM/F12 (Gibco). Glands were partially fragmented by trituration and centrifuged at 500x  g for 5 min. The supernatant was removed and glands resuspended in ice‐cold phenol‐free Matrigel (Corning) at a 1:20 ratio (glands: Matrigel). The gland‐Matrigel mixture was deposited as 20 μL droplets in 48 well tissue culture plates and allowed to set at 37°C for 15 min. Droplets were overlaid with expansion medium (Table  S4 ) and maintained at 37°C under 5% CO 2 . Organoids were cultured to passage 2/3 before hormone treatment; organoids were passaged and allowed to establish for 48 h. Media was supplemented with 10 nM β‐oestradiol (E2, Sigma Aldrich) for 24 h, followed by 10 nM E2, 1 μM medroxyprogesterone acetate (MPA, Sigma Aldrich) and 500 μM cyclic adenosine monophosphate (cAMP, Sigma Aldrich) for 5 days. Organoids were collected using Cell Recovery Solution (Corning) and fixed with 10% neutral buffered formalin for 1 h. Fixed organoids were embedded in HistoGel (Fisher Scientific) and embedded in paraffin. Immunohistochemistry (IHC) was performed as previously described [ 2 ]. FFPE organoids were sectioned at 3 μm, dewaxed, and rehydrated prior to heat‐induced epitope retrieval in a pressure cooker. Endogenous peroxidase activity was blocked with 0.3% hydrogen peroxide (Thermo Fisher Scientific, Runcorn, UK) for 10 min, prior to incubation with diluted primary antibody overnight at 4°C (Table  S5 ). The appropriate ImmPRESS polymer‐based system was applied for 30 min at room temperature, and visualization was with ImmPACT DAB, following the manufacturer's instructions (Vector Laboratories). Sections were counterstained using Gill II Hematoxylin (Thermo Fisher Scientific), dehydrated, cleared, and mounted using Consul‐Mount (Thermo Fisher Scientific). Matching isotype IgG antibodies replaced the primary antibody as a negative control, with an internal positive control in each staining run. Slides were digitalized using an Aperio CS2 slide scanner (Leica Biosystems, Milton Keynes, UK). Scoring was assessed by counting the proportion of positively stained organoids for the primary antibody versus the proportion of negatively stained organoids. GraphPad Prism software (version 9.0, GraphPad Software, San Diego, CA, USA) was used to statistically analyze paired qPCR data, using the Wilcoxon matched‐paired signed rank test and the level of statistical significance set at p ‐value < 0.05.

Figure S1. Figure S2. Figure S3. Table S1. Table S2. Table S3. Table S4. Table S5. Table S6.