The single-cell atlas of cultured human endometrial stromal cells.

Two primary HESCs, P1 and P2, that were isolated from the human endometrial biopsy of two donors at mid-secretory, were kindly provided by Dr. Steve Young at the University of North Carolina at Chapel Hill. One primary HESC, P3, was collected from one mid-secretory donor at Baylor College of Medicine. The Telomerase-immortalized endometrial stromal cell line (THESC) was derived from mid-secretory HESCs and generously provided by Dr. Charles J. Lockwood ( 17 ). The primary HESCs and the THESCs were cultured at 37 °C under 5% CO2 in Dulbecco’s modified Eagle’s medium (DMEM-F12, 11320033, Fisher), supplemented with 1X Penicillin-Streptomycin (15140122, Fisher) and 5% Fetal Bovine Serum (10–437-028, Fisher). Primary HESCs at passage 4–5 were collected for scRNA-Seq analysis. Both the primary cells and THESCs were rinsed with 1X PBS and incubated by 0.25% trypsin (2520056, Fisher) at 37 °C to dissociate the cells from the cell plates. After 2–4 min incubation when the majority cells were dissociated, the cell solutions were gently pipetted and collected in the 15ml tubes. Then, the trypsin was neutralized by culture medium at the ratio 1:3. The cell solutions were centrifuged at 300g for 5 min and the supernatant was removed. The cell pellets were resuspended in cell culture medium and filtered using a 40μm cell strainer. The filtered cell suspension was kept on ice until loaded onto the machine. The cells were counted and examined for viability with trypan blue staining using a TC-20 cell counter (Bio-Rad). Approximately 10,000 live cells at 1×10 6 cells/ml concentration were loaded into the Single Cell Chip to generate single cell emulsion in the Chromium Controller with Chromium Single Cell 3’ Library & Gel Bead Kit v3.1 (Cat. 1000268, 10x Genomics). Reverse transcription of mRNA and cDNA amplification were carried out following the manufacturer’s instruction (Cat. 1000268, 10x Genomics). The amplified cDNA was further fragmented to construct NGS libraries. Then, Real-time PCR Seq and Nova-Seq in NIEHS Epigenomics and DNA Sequencing Core Laboratory with the parameters recommended in the manufacture’s instruction manual. 2% PhiX control library was added in the pooled library for quality control. Raw FASTQ files were processed by Cell Ranger software (version 4.0.0, 10× Genomics) using “cellranger count” pipeline, which used STAR for alignment against human reference genome GRCh38/hg38. This pipeline used the barcodes to generate feature-cell matrix of unique molecular identifier (UMI) counting. All the primary HESCs and THESC were analyzed together by Seurat package (version 3.6.3) ( 18 ) with or without integration by Harmony ( 19 ). After filtering and normalization, the combined datasets were clustered following the standard unsupervised workflow. The marker genes of each cluster were identified by pairwise comparisons of each cluster against all other clusters using the FindAllMarkers function in the Seurat package with adjusted p-value less than 0.05. Expression of selected marker genes from each cluster was visualized by Dimplot in the Seurat package or Vision software ( 20 ). The heatmap of the top five marker genes from each cluster was generated according to the Seurat pipeline. The unsupervised hierarchy clustering based on the DEGs from all the clusters was generated using Partek® Genomics Suite® software, v7.0. The DEGs between primary HESCs and THESCs in the mature, proliferative and active fibroblasts were identified by pairwise comparisons using FindMarkers in the Seurat package. DEGs with adjusted p value 1.3 were used for Ingenuity Pathway Analysis (IPA, Qiagen) to identify the enriched pathways in the 13 cell clusters. After harmony integration, DEGs with adjusted p value 1.2 were used for Ingenuity Pathway Analysis (IPA, Qiagen) and DAVID Functional Annotation Bioinformatics Microarray Analysis to identify the enriched pathways ( 21 , 22 ). The cell cycle was predicted using CellCycleScoring in the Seurat package. Cell trajectory of each sample was analyzed using monocle3 ( 23 ). The overall communication map and dominant ligand-receptor pairs within the three primary HESCs and one THESC were analyzed using CellChat ( 24 ). Only the clusters with more than 50 cells were included for this analysis. The overall communication score was calculated by computeCommunProbPathway. The top ligand-receptor pairs were visualized by netVisual_bubble. The scRNA-Seq of our combined dataset was mapped to the published scRNA-Seq dataset from human endometrial tissues at mid-secretory ( 25 ) using the ClusterMap ( 26 ). The detailed script was in supplement material “HESC-scRNAseq-script”. The three primary HESCs and one THESC cell were seeded at 6-well plate, Cells from the same sample were pooled together for RNA isolation. The total RNA was isolated by Rneasy mini kit (74104, Qiagen). 1ug total RNA was reverse transcribed into cDNA using M-MLV reverse transcriptase (28025013, Fisher) following the manual. Quantitative real time PCR was performed using SsoAdvanced Universal SYBR Green Supermix (1725274, Bio-rad). SYBR green primers were designed using NIH Primer-blast or downloaded from the Primerbank ( 27 ) and synthesized by Sigma-Aldrich ( Excel S1 ). ΔΔCt values were calculated using the Ct number of 18S as amplification control and THESC as the base line to acquire relative mRNA levels per sample. The raw and processed data for scRNA-Seq of the three primary HESCs and one THESC were uploaded to GEO as GSE195918 .

In order to determine the cell atlas of the cultured human endometrial stromal cells, we performed single-cell RNA-Seq in three patient HESCs and one THESC. In total, 1,688, 2,189, and 2,219 HESCs from three different donors and 1,672 THESCs were collected for scRNA-seq analysis, respectively. On average, 6,734 genes and 54,610 UMI per cell, 6,480 genes and 49,368 UMI per cell, 6,400 genes and 46,638 UMI per cell, 6,588 genes and 51,024 UMI per cell were detected in the three primary HESCs and one THESC samples. Based on the number of cells, the average gene per cell and the average UMI per cell, the quality of our data in four samples was consistent and was comparable to the recently published single cell studies using human endometrial tissues and cells ( 11 , 13 , 15 ). First of all, we combined all the four samples for the scRNA-Seq analysis. To check the inter-sample differences, we reassigned one color to the cells from the same sample in the combined UMAP ( Fig. 1A ). Generally, the four samples were separated from each other suggesting the variations between the samples. Among them, the three primary HESCs were closer to each other with a small number of cells overlapped, while the THESC was further away from the primary HESCs with very few cells overlapped ( Fig. 1A ) suggesting that the primary cells still presented the individual differences with the bigger changes observed between primary and transformed HESCs. For the combined datasets, there were 13 distinct cell clusters labeled with different colors ( Fig. 1B ). Based on the expression patterns of the marker genes, we performed hierarchy cluster analysis and grouped the 13 cell clusters into three categories that we named mature, proliferative, and active fibroblasts ( Fig. 1C ). Each category contained cell clusters from different samples ( Fig. 1D ). The mature and proliferative fibroblasts had the cell clusters from all four samples, primary HESC P1, P2, P3 and THESC. Cluster 11 that belongs to two primary HESCs P1 and P3, and cluster 12 from THESC were assigned to the active fibroblasts. Each sample also contained the cells from different clusters. For example, primary HESC P1 was mainly composed of cluster 3, 4, 11, while P2 was consisted of cluster 1, 5, 9, 10. These results suggested that the four HESCs shared not only similar cell composition, but also maintained certain distinct cell clusters. To further understand the characters of the 13 cell clusters, we compared the marker genes and top altered pathways of each cluster ( Fig. 1E , table 1 , Fig, S1 and Excel S2 ). Clusters 0, 1, 2, 3 and 5 are grouped as mature fibroblasts. The mature fibroblasts were characterized as having an abundance of keratins, collagens and cellular matrix proteins, such as KRT17, ACTG2, VCAN. Besides, cluster 3 was enriched with DIO2, that has been associated with recurrent pregnancy loss ( 11 ) and pro-fibrotic factor IGFBP5 ( 28 ); Cluster 0 showed highly expressed retinoid synthesis related genes, such as Aldehyde Dehydrogenase ALDH1A1 and Retinol binding protein RBP1 ( 29 , 30 ); Cluster 2 had higher levels of decidual regulators PRRX2 ( 31 ) and PTX2 ( 32 ); Cluster 1 had abundant Keratin KRT17, that was also detected in the cultured human skin fibroblast ( 33 ), and contractile related smooth muscle actin ACTG2 ( 34 ); Cluster 5 exhibited marker genes, such as extracellular matrix Versican (VCAN) and cell adhesion protein ITGA1, that were critical for the differentiation of the myofibroblast ( 35 , 36 ); Clusters 4, 6, 7, 8, 9 and 10 were grouped as proliferative fibroblasts. In the proliferative fibroblasts, the genes related with cell proliferation were enriched in all the clusters while only the ranking of these genes still varied between different clusters. Specifically, cluster 4 showed relatively higher expression of mitosis markers CCNB1 and CDC20 ( 37 ); The top marker genes in cluster 8 were DNA replication regulators CLSPN and FEN1 ( 38 , 39 ); Cluster 6 was enriched with TM4SF1, that can facilitate cell migration ( 40 ) and the proliferative marker TOP2A; Cluster 7 had more ASPM and AURKA, that plays critical roles in spindle organization during mitosis; Cluster 9 was enriched with proliferative markers TOP2A and MKI67; Cluster 10 had relatively higher expressions of ribosome proteins RPL27A and RPL38A. Finally, clusters 11 and 12 were clusters of active fibroblasts, which expressed several genes critical for cell differentiation. For example, LRRC75A, that labels mesenchymal stromal cells with ossification ability ( 41 ), and EGR1, that was essential for decidualization ( 42 ), were more enriched in cluster 11; Chemokine CXCL12 that mediates epithelial proliferation ( 43 ), and fibroblast growth factor FGF5 were expressed at much higher levels in cluster 12. The expressions of 12 selected marker genes from different cell clusters were validated by Real-time PCR using the pooled total RNA from each sample ( Fig S1B ). Utilizing the marker gene expression in each cluster, differentially enriched pathways were predicted by Ingenuity Pathway Analysis ( Fig. 1E ). Indeed, the cardiac hypertrophy and fibrosis pathways that are related with mature fibroblast were more activated in the mature fibroblast group. The enhanced chromosome replication and inhibited G2/M checkpoint were mainly observed in the proliferative fibroblasts. Multiple signaling including EIF2 (translation) ( 44 ), mTOR ( 45 ), oxidative phosphorylation (mitochondria) ( 46 ), IGF1 ( 47 ), and Wnt/β-catenin ( 48 ) pathways, which are not only critical for maintaining the cellular homeostasis, but also have been associated with the decidualization ability of the HESCs, were more active in the two active fibroblast clusters. Endometrial mesenchymal stem cells (eMSC) were a small subset of colony forming cells, characterized as co-expressions of MCAM and PDGFRB , and single expression of SUSD2 ( 49 ). Previous studies indicated eMSC can be cultured similarly as endometrial fibroblast in vitro but gradually differentiate into stromal fibroblasts ( 50 , 51 ). We also checked eMSC marker expressions in the cultured HESCs and THESCs. To our surprise, 20–80% cells of the primary HESCs, while less than 20% cells of the THESCs in all the clusters except proliferative cluster 10 co-expressed MCAM and PDGFRB ( Fig. S2A ). And the highest co-expression levels were found in mature cluster 1 and 5 belonging to primary HESC P1 ( Fig. S2A ). Similar pattens were observed for SUSD2 with much lower percentages of positive cells ( Fig. S2B ). But the highest level of SUSD2 were detected in the mature cluster 3 which mainly contained the cells from primary HESC P2 ( Fig. S2B ). Eventually, less than 5% MCAM+/PDGFRB+/SUSD2+ cells were detected in most samples, except primary HESC P1 still maintained more than 10% triple positive cells in two mature fibroblast cluster 3 and 0 and one proliferative cluster 9 ( Fig. S2C ). In summary, the eMSC markers were still highly maintained in most cell subpopulations in the cultured primary HESCs and THESCs. Since several cell cycle genes, such as CCNB1, CCNB2, CENPF, CDC20, TOP2A, MKI67, PCNA, were enriched in the proliferative cell clusters including cluster 4, 6, 7. 8, 9, 10, respectively we hypothesized that the 13 subpopulations of fibroblasts may belong to the different stages of the cell cycle. To test it, we performed cell cycle prediction using the cell cycle score function in Seurat and labeled the cells at stage G1, G2/M and S with different colors ( Fig. 2A ). We found more than 80% of the mature and active fibroblasts were at the G1 stage; the proliferative fibroblasts were mainly at G2/M and S stages, expect cluster 10 showed a mixture of cells at all the three stages ( Fig. 2B ). These results unveiled the different cell stages in each cell cluster which may also contribute to the distinct characters of the cell clusters. According to the resemblance of the single cell transcriptome between each cluster, we predicted the cell trajectory in each sample ( Fig. 2C ). In general, if we assumed that the cell trajectory started from the proliferative fibroblast, then it branched into two routes in which one was toward the mature fibroblast, and the other pointed to the active fibroblast suggesting that the mature and active fibroblast may be two developing directions of these fibroblasts. Having observed transcriptomic differences between the primary HESCs and THESC, we next focused on the conserved cell composition. Our scRNA-Seq analysis indicated that the cell clusters in all the cultured HESC samples can be divided into mature, proliferative and active fibroblasts suggesting that a consistent organization of cell subpopulations may exist in the cultured HESCs. Integration across multiple scRNA-Seq datasets have been widely used to identify the conserved clustering of different samples ( 52 ). Therefore, we integrated our scRNA-Seq data from all the cultured HESCs using Harmony ( 19 ) that removed the individual differences of all the samples and only focused on the conserved cell organization across the samples. As expected, the three HESCs and THESC were largely overlapped with each other and shared similar cell clustering patterns after integration ( Fig. 3A ). Transcriptomic annotation identified 10 clusters ( Fig. 3B ). Among them, clusters 0 to 6 had a mixture of cells from the three HESCs and THESC. But cluster 7 was found in primary HESC P1 and P2 cells, cluster 8 was present in all three primary HESCs but not in THESC, and cluster 9 was found in primary HESC P2 and P3 cells, respectively ( Fig. 3C , D ). However, the clusters 7, 8, and 9 were comprised of a small number of cells between 1 to 52, that accounted for less than 2% cells in all the cells. Additionally, the numbers of cells from clusters 0 to 6 were comparable among HESCs and THESC ( Fig. 3D ). Similar to the previous analysis without Harmony integration, three categories of fibroblasts can be defined based on marker gene expression patterns ( Fig. 3E ), including mature fibroblasts (cluster 1, 8, 0, 9), active fibroblasts (cluster 3, 7, 6), and proliferative fibroblasts (cluster 2, 4, 5). The top marker genes of each cluster implied different biological characteristics ( Table 2 , Fig. S3 and Excel S3 ). In the mature fibroblasts, cluster 1 was enriched with muscle contraction gene KRT19 ( 53 ), stress response gene NUPR1 ( 54 ), and stem cell marker ALDH1A1 ( 55 ); Cluster 0 was enriched with IGF1 binding protein IGFBP5, basal cell marker KRT17 ( 56 ), and canonical smooth muscle marker ACTG2 ( 34 ); Cluster 8 only had six differentiated expressed marker genes, that mainly belong to the AP1 family including FOS, ATF3, and JUN ( 57 ); Cluster 9 shared the maker IGFBP5 with cluster 0, stem cell marker ALDH1A1 ( 55 ) with cluster 1, while cluster 9 also presented some unique markers related with aging including SRFBP1 and IGFBP2 ( 58 , 59 ); For the active fibroblasts, cluster 3 was enriched with myofibroblast transition markers, such as ACTA2, VCAN and CRYAB ( 35 , 60 , 61 ); Cluster 7 shared some cluster 3 marker genes including VCAN and ACTA2, but cluster 7 also had some unique marker, such as PTX3 that is a key immune mediator ( 62 ); Cluster 6 also had a unique marker gene LRR75A, that labels mesenchymal stromal cells with ossification ability ( 41 ). Cluster 2, 4 and 5 belongs to the proliferative cells with higher expression of cell proliferative genes. Since MKI67 expression starts to increase at S phase and drastically drops at the later stages of G2/M phase ( 63 ), cluster 2 that showed higher expression of multiple cell mitosis regulators including PTTG1, CDKN3 and CCNB1, but not MKI67 implied that they were more likely at the late stages of G2/M ( 64 – 66 ). Both cluster 4 and 5 had higher levels of MKI67 suggesting they were around S and G2/M stages. Additionally, cluster 4 had abundant expression of chromatin regulators during DNA synthesis, such as HIST1H1B, CLSPN, and ATAD2 ( 67 – 69 ) suggesting that they may be at S phase. In contrast, cluster 5 was highly enriched with spindle regulators, such as CKS2, CENPF, and ASPM ( 70 – 72 ) suggesting that it was at the M phases. Furthermore, we used the CellCycleScoring in the Seurat package to predict the cell cycle stages in the integrated data. It also suggested cluster 2 and 5 were more likely at G2/M phases, while cluster 4 was probably at S phase ( Fig. 3F ). Vision, a tool to describe the coordinated variation for scRNA-Seq data, has been used to interpret the functions of the cell clusters ( 20 ). The E2F and MYC signaling that is critical for DNA synthesis ( 73 , 74 ) were more enhanced at the proliferative fibroblasts, including cluster 2, 4, and 5. The hypoxia signaling and the hypoxia stimulated myogenic signaling ( 75 , 76 ) were more activated in the mature and active fibroblasts ( Fig. 4A ). Similarly, pathway analysis using the DEGs from each cluster also indicated that the cell cycle signaling was upregulated in the proliferative fibroblast. The fibrosis and its related ephrin signaling ( 77 ), the decidualization related integrin ( 78 ), RHOA ( 79 ), the IL8 signals that can be secreted by fibroblast to stimulate the immune cell motion ( 80 ), VEGF that can induce angiogenesis and fibroblast remodeling ( 81 , 82 ) were more activated in the active fibroblasts ( Fig. 4B ). In order to support the homeostasis of each cell culture system, cell communication plays a critical role. Having demonstrated that the four samples contained distinct cell clusters, we next examined the communications of these cell clusters within each sample using the CellChat program ( 24 ). In order to compare the communication between different samples, we still used the cell clusters that was generated without integration. All the samples contained at least one cell cluster from mature, active and proliferative fibroblasts, except primary HESC P2 that only had mature and proliferative fibroblasts ( Fig. 5A ). In total, strong communications between the mature, proliferative and active fibroblasts were observed in the all the samples. More interestingly, the cells located in cluster 10, one proliferative cluster in both P2 and P3, seems to have the least communication with other clusters. The communication probability between cell clusters was predicted based on the expression levels of the ligand-receptor pairs at different clusters and presented as a circle with color changed from red to blue which was proportional to the communication signal strength from strong to weak ( 24 ). We listed the top ligand-receptor pairs that were enriched between clusters and within each sample ( Fig. 5B ). First, we focused on the common communication signals of all four samples. All the ligand-receptor pairs can be divided into two major groups. One group is the extracellular matrix, adhesion and membrane proteins. The other is growth factors and cytokines. Since extracellular proteins were abundantly expressed in HESCs and mediated multiple cellular activities ( 83 ), it was not surprising that collagen ( 84 ), fibronectin (FN1) ( 85 ), laminin ( 86 ), semaphorin ( 87 ), and thrombospondin (THBSs) ( 88 ) and their receptors including integrin, CD44, CD47, Syndecan, Neuropilin, and Plexin were the signals that were mostly transmitted between different clusters within all the samples. Among the growth factors, fibroblast growth factors (FGFs) and their receptors FGFRs that modulate a broad range of cell functions ( 89 ) were the most active signals. HBEGF-EGFRs, an embryo implantation regulating ligand-receptor pair ( 90 , 91 ), Pleiotropin (PTN), a decidual mediator ( 92 ) that may interact with its receptor Nucleolin (NCL) to regulate inflammation ( 93 ), WNT-FDZs signaling, the key regulator in the uterus at both physiological and pathological conditions ( 94 , 95 ), TNFSF12 and its receptor TNFSF12A that promotes the proliferation and collagen synthesis of fibroblast ( 96 ), also showed higher communication scores in all the HESCs. In addition of these common ligand-receptor pairs, some specific communication signals within each sample were also detected. The VEGFs and its receptor VEGFR1 that may play complex roles in angiogenesis and inflammation, and the Delta-like 3 DLL3-NOTCH, an inhibitory NOTCH signaling ( 97 ) were only identified in the three primary HESCs, not the THESC. On the contrary, the ANGPTL2-Integtin ITGA/B signaling that promotes angiogenesis showed higher communication scores in the THESC than the primary HESCs, while the JAG1-NOTCH, an activating NOTCH signaling ( 98 ) was detected in the two primary HESCs P1 and P2, and the THESC. TGFβ-TGFBRs and PDGFs-PDGFRs that can coordinate to regulate the fibroblast proliferation and differentiation ( 99 , 100 ) showed more intensive communication in P2 and P3 than others. Cardiotrophin-1 (CTF1) that shares the same receptors LIFR and IL6ST as LIF and OSM and regulates the cardiac hypertrophy ( 101 , 102 ) was only observed in two primary HESCs P1 and P3. Regardless of the sample identity, the cell clusters from mature, proliferative and active fibroblast also showed variable number, identity and communication probability of the same ligand-receptor pairs. Similar to the pathway analysis ( Fig. 1E ), the mature and active fibroblast showed stronger communication signals of collagens than the proliferative fibroblast. Additionally, IL11 that labels a subset of proliferative fibroblasts in the colorectal tumor ( 103 ) were also found synthesized mainly by the proliferative fibroblasts in our study and may interact with its receptors IL11RA and IL6ST in the mature and proliferative fibroblast. These results indicated the distinct communication signals between different cell clusters and samples. Since the purpose of the in vitro culture of primary HESCs and THESCs is to mimic their functions in vivo , we compared our cultured HESCs from the mid-secretory donors with the published human endometrial scRNA also collected at mid-secretory ( 25 ) using the ClusterMap program ( 26 ). In the published human endometrium at mid-secretory, there were epithelium (Epi), endothelium (Endo), perictyes (Peri), natural killer / T cells (NK/T), and several clusters of stroma fibroblast (Fibr) ( Fig. 6A ). The Circos plot connected the matched groups between our cultured HESCs and the published mid-secretory human endometrium ( Fig. 6B ). Based on similarity of the single cell transcriptome, two major connections were established labeled by the ribbons with green and yellow colors which are the least transparent color in the figure. One was the proliferative cluster 10 from all four cultured HESCs showed higher similarity to the Fibr cluster 1 of mid-secretory endometrium; the other was the active cluster 11 and 12 from the cultured HESCs that were more correlated with Fibr cluster 5 of mid-secretory endometrium. Furthermore, we combined our cultured HESCs and the published human endometrium together, re-run the clustering analysis and identified five major groups in the combined datasets ( Fig. 6C ). Similar to Circos plot ( Fig. 6B ), HESC active cluster 11 and 12 and mid-secretory endometrial cluster 5, while HESC proliferative cluster 10, 2, 7, 4, 8, 6, 9 and Mid-secretory endometrial Fibr cluster 1 were grouped as two groups, respectively. Besides, HESC mature fibr 1 and 5 from P2 + Mid-secretory Fibr 4, HESC mature fibroblast cluster 0 and 3 from P1 and P3 + Mid-Secretory fibr 0, 2, 3, 5 were grouped with each other, respectively. Additionally, all the non-fibroblast clusters from Mid-secretory endometrium, including the epithelium, endothelium, pericytes and immune cells were grouped together, which were quite different from all the stromal clusters from human endometrial and our HESCs. These results suggested the cultured HESCs share similar characters with some sub-populations of human endometrial stromal fibroblasts.

The scRNA-Seq analysis of the three primary HESCs and one transformed THESC that were collected from the donors at mid-secretory phase identified individual differences between the donors and deviations between the primary and transformed HESCs. Each sample contained distinct cell clusters that can be characterized as mature, proliferative and active fibroblast with variable stages of cell cycles and intensive cell communications. More importantly, the cell clusters in the cultured primary and transformed HESCs showed higher correlation with the fibroblast clusters in the human endometrium in vivo suggesting the resemblance of the in vitro culture model and in vivo scenario.

Fibroblasts refer to the interstitial cells that are derived from the mesenchyme lacking the endothelial, immune and epithelial makers. Fibroblasts heterogeneity within and between organs has been extensively studied ( 34 , 104 , 105 ). The fibroblast subpopulations exhibited distinct characteristics that were largely affected by the microenvironment ( 106 ). In the human endometrium, subtypes of stromal fibroblasts have been proposed in multiple studies ( 11 , 12 , 25 , 107 , 108 ). Most studies have focused on the epithelial rather than the stromal subclusters and simply defined the stroma into decidual and non-decidual, and proliferative and non-proliferative clusters. Only one analysis that focused on the stroma compartment in the human endometrium reported more delicate structures within the endometrial stromal fibroblasts including three stromal clusters with no unique gene expression profiles, PAGE4+, ACAT2+, ECM, BMP+, CTNNB1+, THY1+, and ISG15+ clusters ( 107 ). In our study, we successfully identified multiple cell clusters in the primary and transformed HESCs that were cultured in vitro which can be correlated with several subtypes of stromal fibroblasts in the human endometrium in vivo . The application of primary HESCs is limited by the scarcity of the material ( 109 ). Therefore, investigators have been tried multiple methods to establish a HESC cell line to extend their life span. Transfection with SV40 T antigen can immortalize the HESCs but transform them into the neoplastic cell lines with marked karyotypic abnormalities ( 110 ). Ectopic expression of telomerase catalytic subunit (TERT), whose shortening can lead to the cell division inhibition ( 111 ), has successfully generated the stable HESC cell lines that display similar phenotypes of primary HESCs ( 17 , 112 , 113 ). One novelty of our study is to directly compare the scRNA-Seq of HESCs and THESCs. In our scRNA-seq, we still observed distinct differences between each of the three primary HESCs and the THESC, with the THESC being the most divergent of the four analyzed. However, both the primary HESCs and THESC contained the conserved mature, proliferative and active fibroblasts supporting that THESC mimic some features of the primary HESCs. These differences between the primary HESCs and the THESC may stem from the different culture passages or the original donors, or the aftermath of the immortalization methods. The immortalized THESC are a good model for all kinds of genetic manipulations since it can provide consistent results across multiple passages. But the results from the THESC must be further validated with the primary HESCs. Recently, the generation of human endometrial stromal cells from induced pluripotent stem cells is attracting a lot of attention for its therapeutic roles ( 114 ). But there are still extensive knowledge gaps between the fully functional endometrial stromal cells and stem cells. As expected, we also observed individual differences among the three primary HESCs, which were all collected at mid-secretory stages. These individual differences of the primary HESCs were most likely be inherited from the genetic and epigenetic differences of the donors who have distinct genetic polymorphism, experience variable environmental conditions and may undergo certain disease status that has not been checked during biopsy ( 115 ). Additionally, due to the temporal and spatial changes of endometrial tissues ( 25 ), the position of biopsy and the collection timing during mid-secretory may also influence the cell characters and compositions. Our observation of distinguished single-cell transcriptome among the primary HESCs further supported the importance of personal medicine in endometrial disorders, including infertility, endometriosis, fibrosis, endometrial cancers, and emphasized potential of the cultured HESCs as one valuable tool. scRNA-Seq has been widely used in identifying cell identity and comparing transcriptome from cell subpopulations. But due to the sample variance, experimental changes and sequencing platform differences, the comparison of multiple scRNA-Seq is still challenging ( 116 ). In our study, we performed the scRNA-Seq of all four samples in the same experimental condition, but we still observed the undeniable sample differences that stem from individual patients. Actually, similar observations have been reported before, in which the cell type differences between immune, mesenchymal and epithelial were well conserved across the samples, but when the study focused on the subpopulations of one cell type, the mesenchyme, they showed variable clusters according to different development days ( 117 ). Therefore, we think the different fibroblast clusterings that were defined without any integration represent both sample difference and intrinsic fibroblast heterogeneity. In contrast, the fibroblast clusterings identified after integration by Harmony method ( 19 ) can override the sample differences and focus on the fibroblast heterogeneity. As expected, both cell clusterings with or without integration can identify the clusters belonging to mature, proliferative, and active fibroblasts suggesting a similar cell subpopulation existed in the primary and transformed human endometrial stromal cells. Previous scRNA-Seq studies mainly detected the subpopulations in the cultured HESCs undergoing hormone induced decidualization or interacted with embryos ( 11 , 118 , 119 ). During decidualization, the senescent, activated and decidual fibroblasts have been detected simultaneously suggesting HESCs have distinct responses upon hormone treatment ( 11 , 118 ). But it is intriguing to know whether these subtypes of fibroblasts were generated by different cellular responses to the hormones or the variable intrinsic cell identity before hormone treatment. The novelty of our study is to determine the distinct cell populations have already of the primary HESCs and THESC without any hormone treatment. Based on the enriched genes, all these cell clusters can be attributed to three major categories: mature, proliferative and active fibroblast. Among them, the mature fibroblast had abundant collagens and extracellular matrix proteins, the proliferative showed higher proliferation ability, while the active fibroblast was enriched with multiple growth factor and cytokine pathways. Therefore, we hypothesize that the different characteristics of the cells before any hormone treatment may be critical to the distinct decidual responses that were observed later. Further analysis using lineage tracing of these specific cell types can provide more direct evidence. Primary human endometrial fibroblasts were routinely used to check its functions at multiple physiological and pathological conditions. The cultured HESCs normally inherit the donor’s status, such as decidualization defects ( 11 ), endometriosis ( 50 ), menopause stages ( 120 ). But the rapid adaption of fibroblasts with the culture environment often leads to loss or gain of certain characteristics. Along with the passages, it is inherent that the fibroblasts gradually lose proliferation ability ( 121 ). Additionally, ACTA2 expression was normally restricted to the myofibroblast and pericytes, not the fibroblasts, but its levels were elevated in the cultured fibroblasts suggesting a switch to the myofibroblast phenotype ( 122 ). Similarly, keratin, an epithelial marker, has also been reported in the cultured non-transformed skin fibroblast ( 33 ) and transformed skin and lung fibroblasts ( 123 ). Our scRNA-Seq also showed ACTA2, KRT17 and 19 were readily detected in the cultured HESCs and were expressed at a relatively higher level in some subpopulations of cells suggesting some inevitable changes in the cultured HESCs compared to the in vivo stromal cells. Therefore, 2D culture of stromal cells can reflect some characters of the cells, but cautions in interpretation should be taken, as it could not fully recapitulate the features in vivo . Recent advancement of 3D co-culture of endometrial stromal cells and endometrial epithelia organoid models may shed light in establishing the appropriate microenvironment for the fully functional human endometrial stromal cells. Human endometrium undergoes extensive regeneration during menstrual cycle, which are supported by a small subset of endometrial mesenchymal and epithelial stem cells. Among them, eMSC has been successfully isolated using the mesenchymal stem cell marker MCAM/PDGFRB or SUSD2 and cultured in vitro ( 50 , 51 ). In both studies, ordinary endometrial fibroblast medium has been used, thus, the differentiation of eMSC into endometrial fibroblast were continuously observed. Similarly, we also cultured primary HESCs and THESCs in ordinary endometrial fibroblast medium. In contrast to the 1–4% cells in vivo ( 49 ), we found the eMSC markers were highly expressed in almost all the subpopulations in the both the primary HESCs and immortalized THESCs at the percentage from 0 to 80%. A relative higher number of eMSC marker positive cells in the specific primary HESCs while a much lower number in the THESC were detected. Since a gradual loss of mesenchymal stem cell gene signature has been reported along the passages ( 50 ), the reduced eMSC marker expressions in the THESCs may be caused by the higher passage numbers. All these results suggested the presence of eMSC in the cultured HESCs and individual differences still exist among the samples. However, our scRNA-Seq analysis failed to distinguish these eMSC marker positive cells from the rest cells in the same clusters. Further studies using flow sorting to enrich these potential stem cells or scRNA-Seq with more in-depth sequencing may reveal the identity of these cells.

One limitation of our study is the low sensitivity of the scRNA-Seq technology. Therefore, it is possible that we may fail to detect the low expressed genes which may also play crucial roles in the cell homeostasis. The technology development is required to provide a more delicate picture of the single cell transcriptome. The other limitation of our study is we only capture the scRNA-Seq of the cultured cells at one timepoint. Although we predicted the cell trajectory using the software, it is possible that the cultured HESCs may have more dynamic changes that may be another critical point for the cell identity and functions. We propose to perform scRNA-Seq at different timepoints may provide a more comprehensive understanding of the cell behaviors and fates.

Issues involving fertility affect women worldwide ( 1 ). In addition to infertility, as defined as the inability to achieve pregnancy after trying for 1 year, and pregnancy complications, such as preeclampsia, gestational diabetes, and premature birth, represent a significant health concern. As the uterus functions to support embryo implantation, and fetal growth, understanding the molecular mechanisms regulating uterine biology Is critical to understanding issues that arise in pregnancy ( 2 , 3 ). However, the knowledge about the human uterus during pregnancy is very limited. Animal models, such as rodents, ruminants, primates, and especially the genetically modified mouse models are well-developed tools to study uterine function in vivo ( 4 – 6 ). Meanwhile, the in vitro culture of human endometrial cells, including stromal fibroblasts and epithelium, provides a valuable chance to directly investigate the functions of human uterus ( 7 ). Recently, the advancement of human endometrial organoid models assembled by both epithelium and stroma presents a promising avenue in establishing in vitro human endometrium to directly investigate their functions during pregnancy ( 8 , 9 ). The cultured human endometrial stromal cells (HESCs) is one fundamental of the in vitro human endometrial system. HESCs are fibroblasts isolated from the uterus of women and can be maintained in vitro . These cells can be stimulated with ovarian steroids and or cAMP to differentiate into decidual cells. This decidualization process is a requisite for successful embryo implantation, as the decidualization of stromal cells protect and nurture the implanting embryo. HESCs have been studied intensively to identify the complex signals during decidualization that is a critical transformation process that prepares the stromal cells to support the pregnancy ( 10 ). Additionally, distinct fibroblast subpopulations were identified in the human endometrium during the menstrual cycle ( 11 – 13 ), Although cell subpopulations of untreated human endometrial stromal cells have been proposed ( 14 , 15 ), our understanding about the cell composition of HESCs at the non-stimulated status is still very limited. In general, the stromal fibroblasts are easy to culture in vitro ( 16 ). The primary HESCs and the Telomerase-immortalized endometrial stromal cell line (THESC) ( 17 ) are both widely used. Primary HESCs retain the characters of human endometrial stromal cells but are limited by the passages of the cell line while the THESCs have been reported as immortalized cell lines that can decidualize normally upon hormone treatment ( 17 ). THESCs may be a good replacement of primary HESCs. However, a systematic comparison between HESCs and THESCs in the cell composition and transcriptome is prerequisite to unveil the conserved and differentiated characteristics between these two in vitro cell models. To test the hypothesis, we performed scRNA-Seq in three primary HESCs from different donors and one immortalized THESC.

Fig. S1 . Gene expressions in the primary HESCs and THESC. A. The top marker genes in each cluster. From grey to blue, the expression levels were from low to high. B. The real-time PCR of the marker genes in the primary HESCs and THESC. Fig. S2 . Stem cell marker expression in the primary HESCs and THESC. A. The expression patterns of MCAM (red), PDGFRB (green), and co-expression (yellow). The percentages (%) of MCAM +/ PDGFRB + cells in each cluster from each sample. B. The expression patterns of SUSD2 (blue). The percentages (%) of SUSD2 + cells in each cluster from each sample. C. The percentages (%) of MCAM +/ PDGFRB +/ SUSD2 + cells in each cluster from each sample. Fig. S3 . Top marker genes in the primary HESCs and THESC after integration. From grey to red, the expression levels were from low to high.