Reconstructing Lineage Hierarchies of Mouse Uterus Epithelial Development Using Single-Cell Analysis.

The endometrial layer comprises luminal and glandular epithelia that both develop from the same simple layer of fetal uterine epithelium. Mechanisms of uterine epithelial progenitor self-renewal and differentiation are unclear. This study aims to systematically analyze the molecular and cellular mechanisms of uterine epithelial development by single-cell analysis. An integrated set of single-cell transcriptomic data of uterine epithelial progenitors and their differentiated progenies is provided. Additionally the unique molecular signatures of these cells, characterized by sequential upregulation of specific epigenetic and metabolic activities, and activation of unique signaling pathways and transcription factors, were also investigated. Finally a unique subpopulation of early progenitor, as well as differentiated luminal and glandular lineages, were identified. A complex cellular hierarchy of uterine epithelial development was thus delineated. Our study therefore systematically decoded molecular markers and a cellular program of uterine epithelial development that sheds light on uterine developmental biology.

Introduction

Uterine epithelia play key roles in sperm migration, embryo implantation, and fetal nutrition through their physiological functions in nutrient synthesis and transportation (Jeong et al., 2010). Traumatic injury or infection of the endometrium (Al-Inany, 2001), pre-pubertal exposure to hormones such as progesterone (Filant et al., 2012), and dysregulation of the uterus (Dahm-Kahler et al., 2016) could cause absence of uterine glandular epithelia in adulthood and even infertility (Kobayashi and Behringer, 2003). Uterine epithelial adenogenesis begins postnatally, involving budding, tubulogenesis, coiling, and branching of luminal epithelia (Gray et al., 2001). Mechanisms orchestrating these processes are complex. Only a few genes were reported to be involved in the development of uterine epithelia (Kobayashi and Behringer, 2003). It is important to systematically map the molecular and cellular dynamics during early gland development.

Stem/progenitor cells have been reported to be responsible for tissue development and homeostasis, while dysregulation of tissue-resident stem cells may result in dysplasia or diseases (Berry et al., 2016). Evidence for the existence of endometrial epithelial stem cells has been reported recently (Chan et al., 2004, Maruyama, 2014), but little is known about epithelial stem cells during development of the uterus. In-depth understanding of uterine epithelial stem cells could facilitate the development of alternative strategies to manipulate and treat uterus-associated diseases such as dysregulation of the window of implantation and Asherman's syndrome (Dahm-Kahler et al., 2016).

Single-cell analysis technologies have been increasingly utilized to identify rare cell subpopulations (Lanctot, 2015). In this study, single-cell RNA sequencing (RNA-seq) was applied to uterine epithelial cells at five developmental stages, ranging from neonatal to mature stages, to characterize the transcriptional profiles of single cells and to analyze the molecular and cellular dynamics during development of mouse uterus.

Results

Ultra-Shallow Single-Cell RNA-Seq Reveals Cell Heterogeneity during Uterine Epithelial Development

A previous study has shown that development of mouse uterine epithelia is initiated immediately after birth (Gray et al., 2001). In this study, we used CD326 as the epithelial marker (confirmed by immunocytochemistry with anti-CD326 antibody, Figures 1A and 2E) to isolate mouse uterine epithelia from five developmental stages (postnatal day 0 [P0], P7, P14, P28, and P56), and plated each cell within the wells of a 96-well plate using flow-cytometry sorting (Figure 1A). In total, 960 cells were collected. The CEL-Seq protocol was then used to perform single-cell RNA-seq (Jaitin et al., 2014) at an average depth of approximately 20,000 reads per cell (Figure S1). The sequence data were recorded and mapped using TopHat, with the average mapping rate onto annotated genes being about 66%–68% (Figure S1). The final gene expression data were normalized using Cuffnorm according to standard procedures (Trapnell et al., 2012) and represented by tags per million. Genes expressed by less than four cells were filtered (GEO: GSE98451).

Figure 1

Ultra-Shallow Single-Cell RNA-Seq Reveals Cell Heterogeneity during Uterine Epithelial Development

(A) Schematic of the single-cell analysis of five developmental stages, Firstly, the tissues were digested into single-cell suspension; secondly, cells were labeled with CD326-APC antibody and sorted by FACS; thirdly, single-cell RNA-seq was conducted by CEL-Seq; Finally, the data were analyzed and validated.

(B) Heatmap of 714 single cells from the five developmental stages by transcriptome analysis revealed four distinct populations. D, postnatal day.

(C) t-SNE analysis identified cell heterogeneity during the development of uterine epithelia. PND, postnatal day.

(D) Unique marker genes of the four main cell populations. D, postnatal day.

(E) Representative markers for each cell population.

Figure 2

Cell Heterogeneity during Development Exhibits Characteristics of Progenitor and Mature States

(A) Gene ontology (GO) analysis showed elevated cellular response to estrogen stimulus and estrogen receptor signaling during late development of uterine epithelia.

(B) Expression of esr1 gene during development.

(C) Gene ontology analysis showed that positive regulation of cyclin-dependent protein kinase activity was high during the early development of uterine epithelia.

(D) Representative cyclin-dependent protein kinase gene expression during development.

(E) Expression of Ki67 protein in the epithelia (CD326+) during development. Scale bar, 20 μm.

(F) Gene ontology analysis showed that EMT was high during early development of uterine epithelia.

(G) Expression of EMT-related genes during development.

(H) Gene ontology analysis showed that telomerase activity was high during early development of uterine epithelia.

(I) Expression of telomerase activity-related genes during development.

(J) Specific telomerase inhibitors decreased cell proliferation of epithelial cell from early development. Data are presented as mean ± SEM, n = 6 independent experiments. ∗∗p < 0.01, ∗∗∗p < 0.001.

PND, postnatal day.

Both the unbiased clustering analysis and t-distributed stochastic neighbor embedding (t-SNE) analysis showed that single cells from P0 to P56 could be classified into four different clusters (Figures 1B and 1C). These clusters of cells were assigned to different developmental stages of the uterus: the first cluster including cells from P0 and P7 (cluster 1), the second from P14 (cluster 2), the third from P28 (cluster 3), and the fourth from P56 (cluster 4). The data showed that each cell cluster possessed a unique group of expressed genes (Figure 1D), with cluster 1 expressing Wnt7a and Pax2, and clusters 2, 3, and 4 expressing Gata2, Prss29, and Muc4, respectively (Figure 1E).

Cell Heterogeneity during Development Exhibit Characteristics of Both Progenitor and Mature States

To further obtain a functional insight into each cell cluster, we performed gene set enrichment analysis (GSEA) (Table S1). According to the GSEA results, different cell clusters possessed unique gene ontology according to their respective differentially expressed genes. GSEA revealed that cells in cluster 1 possessed characteristics of stem/progenitor cells with low hormone receptor expression, high proliferative capacity, endothelial-to-mesenchymal transition (EMT) characteristics, and higher telomerase activity.

Stem/progenitor cells in hormone-responsive tissues, such as the mammary glands and uterine epithelia, are always negative for expression of the hormone receptor, estrogen receptor α (Esr1) (Giraddi et al., 2015, Janzen et al., 2013, Kato et al., 2007, Masuda et al., 2010). GSEA showed that cells in the early developmental stages (P0 and P7, cluster 1) displayed low levels of cellular response to estrogen stimulus and activation of the intracellular estrogen receptor signaling pathway (Figure 2A), as well as low expression of the Esr1 gene (Figure 2B). However, the cellular response to estrogen stimulus and intracellular estrogen receptor signaling pathway began to increase from P14 (cluster 2), along with the increased expression of Esr1 gene, which correlates with the beginning of uterine epithelial differentiation.

Stem cells are in two phases with respect to the cell cycle, with some undergoing rapid proliferation and others remaining quiescent (Rumman et al., 2015). Gene ontology associated with the cell cycle such as positive regulation of cyclin-dependent protein kinase activity was highly enriched at the early developmental stage (P0 and P7, cluster 1) (Figure 2C), concomitant with the high expression of Cdk genes (Figure 2D). Immunofluorescence staining with anti-Ki67 antibody also confirmed the result that cells at the early developmental stages are highly proliferative (Figure 2E).

EMT is widely regarded to be a key characteristic of stem cells (Battula et al., 2010), particularly in some epithelial tissues, i.e., mammary stem cells (Guo et al., 2012). Gene ontology associated with EMT was highest in P0 cells (cluster 1), but decreased to the lowest level in P14 cells (cluster 2) (Figure 2F). Genes involved in the EMT process such as Smad4, Sox9, Fgfr2, and Tgfbr1 were highly expressed during the early developmental stages (Figure 2G).

High telomerase activity is another crucial characteristic of stem cells (Wong et al., 2010) that support their long-term self-renewal and proliferation (Kong et al., 2014). Gene ontology associated with telomerase activity was highest in cluster 1 (Figure 2H), and genes included in the ontology associated with telomerase activity such as Ptges3 and Dkc1 were found to be highly expressed during the early developmental stages (Figure 2I). Small molecules that inhibit telomerase activity significantly inhibited the proliferation of P7 uterine epithelia (Figure 2J). These results thus indicated the involvement of telomerase during the early development of uterine epithelia. Telomerase activity also decreased dramatically by P14 (cluster 2) implying that uterine epithelia began to differentiate from P14.

Molecular Cascades Regulating the Development and Maturation of Uterine Epithelia

To further investigate the molecular cascade regulating the proliferation, self-renewal, and differentiation of uterine epithelial cells, we used GSEA to identify transcription factors (TFs), signaling pathways, epigenetic activities, and metabolic states highly enriched at the early and late developmental stages.

TFs from the differentially expressed genes of the respective cell clusters were selected (Table S2) and visualized using a heatmap (Figure 3A). Results showed that each cell cluster possessed a unique cluster of TFs. This study thus confirmed the known TFs, as well as identified TFs that are not studied during the development and differentiation of uterine epithelia. For example, Pax2 is one of the known TFs to be involved in the early development of the uterus (Kobayashi and Behringer, 2003). Esr1 is also a known nuclear hormone receptor responsible for estrogen stimulation during maturation. Most of the other TFs identified herein are not studied in the development of uterine epithelia. We found that Gata2 (Rubel et al., 2012, Rubel et al., 2016) was uniquely expressed by P14 cells. Immunocytochemistry was performed to confirm the expression of GATA2 at the protein level, with the results showing the highest expression of GATA2 protein in the P14 cells, with particularly strong expression in the cell nuclei (Figure 3B).

Figure 3

Molecular Cascades Regulating the Development and Maturation of Uterine Epithelia

(A) Differentially expressed transcription factors (TFs) of the respective cell cluster were selected and visualized using a heatmap.

(B) Expression of GATA2 protein in uterine epithelia (CD326+) during development. Scale bar, 20 μm.

(C–K) Sequential upregulation of signaling pathways (C), epigenetic activity (D), and metabolic activity (E) during development. Involvement of WNT (F) and RA metabolic processes (I) during early development of uterine epithelia. Expression of genes related to WNT (G) and RA metabolic processes (J) during development. Specific WNT (H) and ALDH1 (K) inhibitors decreased proliferation of epithelial cells from early development.

(L) Expression of ALDH1A1 protein in the epithelia (CD326+) during development. Scale bars, 10 μm or 20 μm.

PND, postnatal day. Data are presented as mean ± SEM, n = 6 independent experiments. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001. See also Figure S2.

Distinct signaling pathways (Figure 3C) were involved during the early development and maturation of uterine epithelial development. Genes associated with functional entities of WNT (Figures 3F and 3G), retinoic acid (RA) receptor signaling, small guanosine triphosphatase (GTPase), bone morphogenetic protein (BMP), Nodal, fibroblast growth factor (FGF), epidermal growth factor receptor, Insulin, and Hippo signaling pathway were highly enriched in cluster 1 (P0–P7) (Figure 3C). On the other hand, genes associated with innate immune response-activating signal transduction, pattern recognition receptor signaling pathway, Toll-like receptor and apoptotic signaling pathway, NIK/nuclear factor κB signaling, JAK-STAT, and vascular endothelial growth factor receptor signaling pathway were distinctly upregulated in the mature uterine epithelial cells (clusters 2–4) (Figure 3C). Small-molecule antagonists against WNT (Figure 3H), small GTPase, BMP, FGF, JNK, and Insulin signaling (Figures S2A–S2E) inhibited the proliferation of P7 uterine epithelial cells. These results thus confirmed the involvement of these signaling pathways during the early development of uterine epithelia.

During the development and maturation of uterine epithelial cells, dynamic epigenetic changes occur at various levels, with RNA, DNA, and histone modifications (Figure 3D). RNA methylation activity peaked in P7 uterine epithelial cells, and dramatically decreased to the lowest level in P14 uterine epithelial cells. Both DNA methylation and demethylation were observed to be highly active during the early developmental stages of uterine epithelial cells, with DNA methylation peaking in P7 cells and DNA demethylation peaking in P0 cells. Both histone methylation and demethylation were highly active in P7 uterine epithelial cells (Table S3), while histone H3 and H4 acetylation were highly active in P0 uterine epithelial cells (Table S3). Histone H3 deacetylation was observed to be highly active in P0 and P28 uterine epithelial cells, while histone H4 deacetylation peaked in P7 cells (Table S3). Small molecules that inhibit histone methyl transferase and histone acetyltransferase significantly inhibited the proliferation of P7 uterine epithelia (Figures S2F and S2G).

Metabolic activities associated with the glycolytic process, the tricarboxylic acid cycle, and fatty acid oxidation were predominant during the early development of uterine epithelia, which correlated with rapid cell proliferation (Figure 3E). Nucleotide acid metabolic processes peaked during the early development of uterine epithelia, which may be due to active DNA replication in proliferating P0 and P7 cells. Processes associated with carbohydrate metabolism, proteoglycan metabolism, and extracellular exosomes were highly active in mature uterine epithelial cells (P28 and P56) (Figure 3E), which may be due to the fact that mature uterine epithelial cells secrete nutrients and growth factors, such as carbohydrates and glycoproteins that are essential for embryo implantation and growth (Filant and Spencer, 2013). In addition, fatty acid metabolic processes, triglyceride metabolic processes, lactate metabolic processes (Figure 3E), RA metabolic processes (Figure 3I), and related gene expression (Figure 3J) peaked in P14 cells. Immunofluorescence staining with anti-ALDH1A1 antibody also confirmed these results (Figure 3L). Small molecules that inhibit ALDH1 of the retinoic acid metabolic processes significantly inhibited the proliferation of P7 uterine epithelia (Figure 3K).

Aldehyde Dehydrogenase 1A1 Enriches for an Epithelial Stem/Progenitor Cell-like Subpopulation at the Early Developmental Stage of Mouse Uterus

To further investigate the heterogeneity of uterine epithelia at early developmental stages, we performed single-cell qPCR with cells from P7 mouse uterus. Ninety-six pairs of primers of candidate stem cell markers, cell surface marker genes, uterus-specific TFs, housekeeping genes, and epithelial markers were selected (Table S4). Ninety-six CD326+ epithelia were randomly isolated from the P7 mouse uterus by fluorescence-activated cell sorting (FACS) and subjected to single-cell pre-amplification with primers specific for the aforementioned genes. Subsequently, single-cell qPCR was performed.

Cluster analysis was performed with the threshold-crossing (Ct) value of each gene from individual cells (GEO: GSE98438). Postnatal endometrial epithelia could be clustered into different subpopulations (Figure 4A). A unique small cluster of cells expressed high levels of the Aldh1a1 gene (Figure 4B). Furthermore, Aldh1a1 expression was highly correlated with stem cell markers such as Lgr5, Cd133, and Msi1, key TFs such as Pax2 and Emx2, and Wnt signaling genes such as Wnt7a and Ctnnb, as well as other genes such as Ezh2 and Ki67 (Figure 4C). Gene ontology analysis using genes correlated with Aldh1a1 showed association with stem cell proliferation, Wnt signaling, positive regulation of mesonephros development, and morphogenesis of branching structures (p < 0.05) (Table S5). These results thus indicated the key roles of these genes and Aldh1a1+ cells in gland formation during early uterine development. Immunohistochemistry showed that ALDH1A1 is highly expressed in the glands as well as in the invagination sites responsible for gland formation in P7 mice. These sites also expressed high levels of aldh1a1-correlated proteins such as Ki67, EZH2, and PAX2, thus implying that these cells were rapidly proliferating (Figure 4I). We concluded that Aldh1a1 was strongly associated with endometrial epithelial stem/progenitor cells in the early P7 mouse uterus.

Figure 4

Aldehyde Dehydrogenase 1A1 Enriches for an Epithelial Stem/Progenitor Subpopulation during Early Development of the Mouse Uterus

(A) P7 endometrial epithelial cells were clustered into different subpopulations by hierarchical clustering.

(B) Aldh1a1 is highly expressed in a unique cluster of cells.

(C) Genes correlated with Aldh1a1 expression were ordered according to the respective correlation coefficients.

(D) FACS isolation of CD326+ ALDH1A1-high cells from a single-cell suspension of P7 mouse uterus.

(E) CD326+ ALDH1A1-high cells display high expression of aldh1a1-correlated genes.

(F) CD326+ ALDH1A1-high cells possess high self-renewal capacity with more colonies formed.

(G) Comparison of the colony formation efficiency of CD326+ ALDH1A1+ and CD326+ ALDH1A1− cells. Data are presented as mean ± SEM, n = 3 independent experiments. ∗∗p < 0.01.

(H) ALDH1A1 is highly expressed by some epithelial cells within the colony in vitro. Scale bar, 100 μm.

(I) ALDH1A1 is highly expressed in the glands as well as at the invagination sites responsible for gland formation, and is co-expressed together with high levels of Ki67, EZH2, PAX2, and VIMENTIN. Scale bar, 50 μm.

ALDH1A1+ cells were isolated by FACS using the commercial ALDEFLUOR kit (STEMCELL Technologies). We chose the marker CD326 to label the uterine epithelia. The results showed that the ALDEFLUOR kit could be used to label, detect, and isolate ALDH1A1-high cells from endometrial epithelia in P7 mice. FACS analysis showed that 8% of cells are ALDH1A1+CD326+ cells in the P7 mouse uterus, with 42.78% of epithelial cells being ALDH1A1-high cells (Figure 4D). qPCR showed that CD326+ ALDH1A1-high cells expressed higher levels of Aldh1a1-correlated genes such as Vimentin, Lgr5, Cd133, and Ki67 (Figure 4E) compared with the CD326+ ALDH1A1-low cells. High Vimentin expression implied that CD326+ ALDH1A1-high cells of P7 mouse uterus may be involved in EMT.

Next, we compared the self-renewal capacity of the ALDH1A1-low CD326+ with those of the ALDH1A1-high CD326+ cells. The results showed that CD326+ ALDH1A1-high cells have higher self-renewal capacity, exhibiting more colony formation within the culture dish compared with the ALDH1A1-low CD326+ cells (Figures 4F and 4G), with a small number of ALDH1A1+ cells still remaining in the colony after culture (Figure 4H).

These results thus showed that ALDH1A1+ epithelial cells are highly proliferative with high self-renewal capacity, which are the classic hallmarks of stem/progenitor cells (Chan et al., 2004), and these cells are comparable with previously reported endometrial epithelial stem/progenitor cells (Table S6). Hence, we concluded that ALDH1A1 enriches for a population with stem/progenitor properties that are mainly located in the developing uterine gland.

Differentiation of Luminal and Glandular Epithelia during the Maturation of Uterine Epithelia

The initiation of uterine epithelial cell differentiation in P14 cells is associated with decreased cell proliferation rate and a significant increase in the expression of ESR1 gene (Figures 2A–2C and 5D). To further investigate cellular heterogeneity during maturation of uterine epithelial cells, we employed t-SNE analysis to further cluster cells from P14 to P56 into subgroups (Figures 5A–5C). The heatmap of distinct cell clusters (Figure 5A) showed the respective marker genes (Figure 5B), among which cluster-1 cells displayed high expression of Hist1h2ao, a gene essential for cell proliferation, cluster-2 cells expressed high levels of Gata2 and Lpar3, cluster-4 and cluster-5 cells expressed high levels of Prss28 and Prss29, and cluster-6 cells expressed high levels of Lif (Figure 5B). According to a previous study (Filant and Spencer, 2013), Lpar3 gene is a known marker of luminal epithelia (Figures 5E and 5F), while Prss28, Prss29, and Lif are markers of glandular epithelia (Figures 5G–5I). Hence, cellular heterogeneity during the differentiation of uterus cells arose from the differentiation of luminal and glandular cells since P14 (Figure 5J).

Figure 5

Differentiation of Luminal and Glandular Epithelia during the Maturation of Uterine Epithelia

(A) Heatmap of single cells from the late developmental stages revealed six populations.

(B) Unique marker genes of the six main cell populations, cluster 1 (C1) to cluster 6 (C6).

(C) t-SNE analysis identified cell heterogeneity along the late developmental stages of uterine epithelia. PND, postnatal day.

(D–I) Expression of representative markers selected for each cell population.

(J) Differentiation of luminal (LE) and glandular (GE) epithelia occurs during the maturation of uterine epithelia.

Mapping the Cellular Hierarchy of Uterine Epithelial Cells during Development by Lineage Progression Analysis

Cell lineage hierarchy was reconstructed using the SPADE software suite. TF genes, proliferation genes, luminal and glandular marker genes, Wnt genes, stem cell markers, and differentiation genes (Table S7) that were differentially expressed were selected from each cell cluster to reconstruct the lineage SPADE tree, so that different cell clusters from distinct differentiation states were spontaneously separated into their respective branches. The distance between different branches shows the lineage distance between different unique clusters of cells (Qiu et al., 2011).

Branches that have high expression of both Hmga1 and Tead2 genes could be defined as the early progenitor clusters (Figures 6B and 6C), among which the branch that highly expressed Hist1h2ao gene could be defined as the fast-proliferating cell subcluster of the early progenitors (Figure 6D). Moreover, within the fast-proliferating Hist1h2ao progenitor branch, the one that highly expressed Aldh1a1 was a subpopulation of stem/progenitor cells (Figure 6E). As the P0 mouse uterus epithelial cells have a very limited number of ALDH1A1+ cells, after the gland epithelial cells began to invaginate and bud into the mesenchymal layer, the number of ALDH1A1-positive cells started to increase. Thus, we speculated that the ALDH1A1-negative subpopulations of the early progenitors (marked with Hmga1+Tead2+Hist1h2ao+ALDH1A1−Esr1− on the SPADE map) might give rise to ALDH1A1+ progenitors (marked with Hmga1+Tead2+Hist1h2ao+ALDH1A+Esr1− on the SPADE map). In addition to the ALDH1A1− progenitor and the ALDH1A1+ fast proliferative progenitor cells, there is a subpopulation in the early development of mouse uterus marked with Hmga1+Tead2+Hist1h2ao−ALDH1A1−Esr1− in the SPADE map, which we term “non-cycling progenitor,” that might be an intermediate between ALDH1A1+ progenitor and Esr1+ fast proliferative differentiated cells (Figure 6J).

Figure 6

Mapping Cellular Hierarchy of Uterine Epithelial Cells during Development and Maturation by Lineage Progression Analysis

(A) Esr1+ lineage represents mature uterine epithelial lineage.

(B and C) Hmga1 (B) and Tead2 (C) represent early progenitor lineage.

(D) Hist1h2ao represents fast-proliferating lineage.

(E) Aldh1a1 represents a subpopulation of stem/progenitor lineage of the early uterine epithelia.

(F) Lpar3 represents differentiated luminal epithelial lineage during the maturation of uterine epithelia.

(G) Gata2 is highly enriched in the mature Lpar3 lineage.

(H) Prss29 represents differentiated glandular epithelial lineage during the maturation of uterine epithelia.

(I) Reconstructed cellular lineage hierarchy from early stem/progenitor to late differentiated luminal and glandular epithelia.

(J) A work model of the uterine epithelial lineage map.

Branches that strongly expressed the hormone receptor gene Esr1 could be defined as the maturation/differentiation cluster (Figure 6A). This cluster could be further separated into sub-branches: the proliferative subcluster branches displayed high expression of Hist1h2ao (Figure 6D), while high expression of Lpar3 and Prss29 was observed in the luminal epithelial subcluster branches (Figures 6F and 6G) and glandular epithelial subcluster branches, respectively (Figure 6H). These results also indicated that in addition to the early progenitors and differentiated luminal and glandular epithelia, there is a lineage that expressed both Hmga1 and low levels of Esr1 (Figures 6A, 6B, and 6D), which may represent a transitional lineage between the early progenitors and the late differentiated luminal and glandular epithelia.

In summary, the P0 cells mainly comprised the early ALDH1A1− progenitor (Hmga1+Tead2+Hist1h2ao+ALDH1A1−Esr1−) and a few ALDH1A1+ progenitor cells (Hmga1+Tead2+Hist1h2ao+ALDH1A1+Esr1−). Most of the ALDH1A1+ progenitor cells were mainly located in the gland of the P7 mouse uterus epithelial cells, while some of the other cells of P7 mouse uterus were early ALDH1A1− progenitor and late quiescent progenitor (Hmga1+Tead2+). The P14 cells mainly comprised the differentiated luminal cells (Esr+Lpar3+) and the rest of the cells were proliferative differentiated cells (Esr1+Hmgal+Hist1h2ao+) and differentiated glandular cells (Esr1+Prss29+). The P28 and P56 cells mainly comprised the proliferative differentiated cells (Esr1+Hmgal+Hist1h2ao+), differentiated luminal cells (Esr1+Lpar3+), and differentiated glandular cells (Esr1+Prss29+).

Finally, we reconstructed a theoretical model of the uterine epithelial lineage map from SPADE (Figure 6I), which may arise from fast-proliferating ALDH1A1-negative progenitors that transit into fast-proliferating ALDH1A1-positive progenitors, and through “non-cycling progenitors” into the fast-proliferating Esr1+ transitional lineage, ending at the differentiated luminal and glandular cells (Figure 6J).

Discussion

In this study we generated four useful datasets: we provided a single-cell transcriptome profile of the uterine epithelia from neonatal to sexually mature mice in vivo, identified a subpopulation of stem/progenitor cells during early development, revealed the molecular cascades orchestrating their development and maturation, and reconstructed a complex temporal and spatial cellular hierarchical map.

Single-Cell Transcriptome Profile of Uterine Epithelial Development In Vivo

Single-cell analysis technology has been increasingly applied to study the molecular and cellular functions of rare cell populations, particularly stem/progenitor cells (Krieger and Simons, 2015). This study altered the sequencing depth limit for the single-cell RNA-seq study. For full-length cDNA library construction and sequencing, the limit of sequencing depth was previously shown to be 50,000 reads per single cell to accurately classify different cell types (Pollen et al., 2014). However, the 3′ end sequencing adopted in this study showed that shallower sequencing (20,000 reads per cell on average) would be sufficient to classify cells at distinct developmental stages.

Our study provides a single-cell RNA-seq dataset and comprehensive overview of transcriptome dynamics underlying uterine epithelial development and maturation in vivo. We validated the single-cell RNA-seq dataset at multiple levels using single-cell qPCR, qPCR, immunofluorescence, FACS, cell culture, and small molecules. We confirmed the known markers and identified more markers during the development of uterine epithelia.

Epithelial Stem/Progenitors during Early Development of Uterine Epithelia

Stem/progenitor cells have been reported to be responsible for tissue development and homeostasis (Berry et al., 2016). Since the report of clonogenic and highly proliferative endometrial epithelial stem/progenitor cells in 2004 (Chan et al., 2004), various subpopulations have been identified, such as side population cells (Cervello et al., 2011, Chan and Gargett, 2006), label retaining cells (Chan and Gargett, 2006, Patterson and Pru, 2013), SSEA1+ cells (Valentijn et al., 2013), and CD44+ epithelial cells (Janzen et al., 2013). Nevertheless, all of these studies focused on the characteristics of uterus epithelial stem/progenitor cells in adulthood. In this study, we have identified and characterized ALDH1A1 as a marker that enriches for a population with stem/progenitor properties during early development of the uterine epithelia (P7) with high proliferative potential, self-renewal capacity, ability to undergo EMT, and low hormone-response activity. The Aldh1a1 gene identified as the marker for epithelial stem/progenitor cells in the P7 uterus is also known to be a stem cell-associated marker in many other normal tissues and cancers (Ma and Allan, 2011).

Most of the tissue-specific stem cells have two states with regard to cell cycle to maintain tissue homeostasis (Rumman et al., 2015); for example, in the small intestine there are Bmi1+ quiescent stem cells and Lgr5+ proliferative stem cells (Li et al., 2014a). Epithelial stem/progenitor cells in the uterus tissue may also have two states; the SSEA1+ epithelial stem cells in human endometria were reported to show a significant reduction in Ki67 expression compared with SSEA1− cells (Valentijn et al., 2013), while CD44+ epithelial stem cells in the mouse uterus were shown to express high levels of cell division-related genes (Janzen et al., 2013). The ALDH1A1+ subpopulation in this study was mainly located in the newly formed glands of mice during the fast glandular genesis of the uterus, indicating that the ALDH1A1+ cells would be highly proliferative in supporting the formation of new glands.

Molecular Cascades Orchestrating Uterine Epithelial Development

Molecular mechanisms orchestrating uterine development and differentiation are highly complex (Kobayashi and Behringer, 2003). Our study provided a systematic and continuous analysis of molecular cascades during development and maturation of the uterine epithelia. We identified developmental stage-specific signaling pathways, epigenetic and metabolic activities, as well as the respective TFs involved during different stages of uterine epithelial development. In this study, we have identified molecular cascades that were previously uncharacterized. Among these we found that TF Gata2 was uniquely expressed by P14 cells. GATA2 regulated the expression of the gene regulatory network required during progesterone-induced differentiation of healthy endometrial cells by direct binding to their promoter during pregnancy (Rubel et al., 2012, Rubel et al., 2016). In addition, GATA2 was reported to directly bind to the CDKs Esr1 and Esr2 genomic regions by chromatin immunoprecipitation sequencing using endothelial cells (Kanki et al., 2011), with the decrease in Gata2 expression promoting proliferation, migration, and invasion of HepG2 cells (Li et al., 2014b). It was also reported that when Gata2 expression was repressed by hypermethylation, endometriosis occurred (Dyson et al., 2014). A possible explanation for this is that GATA2 may inhibit cell proliferation during uterine epithelial development and induce maturation of the uterus through direct binding to the promoter regions (Cdks and Esr1). This is in accordance with our observation that cells in P14 ceased to proliferate and began to express Esr1, which may imply the involvement of Gata2.

Overall, the molecular cascades described here are defined by unique TFs and markers, stemness gene ontologies, signaling pathways, and epigenetic and metabolic regulation (Figure 7), which interact and cross-talk to ensure proper development and maturation of uterine epithelia.

Figure 7

Schematic Summary of Molecular Signatures Regulating Mouse Uterine Epithelial Development and Maturation

The molecular cascades described here include developmental specific markers and TFs, stemness gene ontologies, signaling pathways, and epigenetic and metabolic regulation. PND, postnatal day; EMT, epithelial-to mesenchymal transition; GPCR, G-protein-coupled receptor; cGMP, cyclic guanosine monophosphate; BMP, bone morphogenetic protein; TLR, Toll-like receptor; NFkB, nuclear factor κB; TCA, tricarboxylic acid; RA, retinoic acid.

Cellular Hierarchical Map of Uterine Epithelial Development

Our study provides a complex cellular hierarchical map of uterine epithelial development and maturation (Figure 6J). We reconstructed the lineage hierarchical map of uterine epithelial development using SPADE. The map represents both temporal (progenitor and differentiation/maturation lineage) and spatial (luminal and glandular lineage) relationships among epithelial cells during the development and maturation of the mouse uterus.

Experimental Procedures

Preparation of Cell Suspension from Mouse Uterus

Female ICR mice (purchased from the Zhejiang Academy of Medical Science) of different ages (P0, P7, P14, P28, and P56) were utilized in this study. Uterus tissues were prepared as described previously (Masuda et al., 2007). Mice were euthanized by cervical dislocation, and whole uterus were explanted by microdissection and immersed immediately into ice-cold PBS (Sigma), followed by removal of the surrounding tissues such as blood vessels, ovaries, and oviducts. All procedures were performed with approved protocols in accordance with guidelines of the animal experimental center of Zhejiang University.

Single-cell suspensions of uterus tissues were then prepared by enzymatic digestion overnight at 4°C with type I collagenase (Life Technologies) diluted in low-glucose DMEM (Gibco), and the digested tissues were triturated into a single-cell suspension with a 1-mL pipette. The enzymes and cell debris were removed by centrifugation at 200 × g for 5 min and washed with PBS, followed by filtering the single-cell suspension through a 70-μm strainer (Corning). The resultant single-cell suspension could then be utilized for flow-cytometry analysis, single-cell analysis, and cell culture.

Fluorescence-Activated Cell Sorting of Uterine Cells

Uterine epithelial cells at different stages (P0, P7, P14, P28, and P56) were isolated by FACS of the uterus single-cell suspension labeled with an anti-CD326 antibody (CD326-APC, eBioscience #17-5791), and cells without labeling antibody were used as negative control. Single CD326+ cells were sorted into each well of a 96-well plate (Bio-Rad) on BD Influx (Becton Dickson) at Hangzhou Normal University. The sorted cells were further used for single-cell analysis.

CD326+ ALDH1A1-high cells from P7 mouse uterus were isolated by FACS of the mouse uterus single-cell suspension labeled with the anti-CD326 antibody and the ALDEFLUOR kit (Stem Cell Technologies); cells treated with an ALDH1 inhibitor provided in the kit were used as a negative control and sorted on BD Influx (Becton Dickson) at HangZhou Normal University. The sorted cells were further used for cell culture.

High-Throughput Microfluidic Real-Time Single-Cell qPCR

Single-cell qPCR was performed according to standard protocols provided by Fluidigm. In brief, 96 genes were chosen for analysis by single-cell qPCR (Table S4). Multiplexed primers for single-cell qPCR were designed with Primer5 and synthesized by Generay. Single-cell sequence-specific pre-amplification individual primer sets (a total of 96) were pooled to a final concentration of 0.1 mM for each primer. Individual CD326+ uterine epithelial cells sorted in the 96-well PCR plate were loaded with 5 μL of RT-PCR master mix (2.5 μL of Cells Direct reaction mix [Invitrogen]; 0.5 μL of primer pool; 0.1 μL of RT/Taq enzyme [Invitrogen]; 1.9 μL of nuclease-free water) in each well. The plates were immediately placed on the PCR machine after centrifugation. Cell lysis and RT were performed followed by reverse transcriptase inactivation, and Taq polymerase activation that was achieved by heating to 95°C for 3 min. Subsequently, the resultant cDNA went through 20 cycles of sequence-specific amplification by denaturing at 95°C for 15 s, annealing, and elongation at 60°C for 15 min. Pre-amplified products were diluted prior to analysis. Amplified single-cell samples were analyzed with the Universal PCR Master Mix (Applied Biosystems), EvaGreen Binding Dye (Biotium), and individual qPCR primers using 96.96 Dynamic Arrays on a BioMark System (Fluidigm). Ct values were calculated using the BioMark Real-Time PCR Analysis software (Fluidigm).

Single-Cell RT, Pre-amplification, Library Preparation, and Sequencing

Sorted cells in the plates were immediately lysed or frozen on dry ice. Single-cell mRNA RT, pre-amplification, and library construction for sequencing were performed according to a previous CEL-seq protocol (Jaitin et al., 2014). In brief, single uterine epithelial cells from the 96-well plate were lysed and reverse transcribed (Invitrogen) by addition of a lysis buffer and a barcoded RT mix, respectively. After RT, the single-cell reaction samples were pooled and digested by ExoI nuclease (NEB), followed subsequently by second-strand synthesis (NEB) and in vitro transcription (NEB). Finally, RNAs were fragmented (Ambion) and the sequencing library was constructed by ligation (NEB), RT (Invitrogen), and PCR (KapaBiosystems) by addition of Illumina primers (with sequencing primer, index, and linkers) to both ends of the fragments. The library was then sequenced using the Illumina HiSeq2500 system.

Computational Bioinformatics Analyses of Single-Cell Data

Single-cell RNA-seq reads data were mapped using TopHat and transcripts were assembled using Cufflinks, with the final gene expression data being normalized using Cuffnorm according to standard procedures (Trapnell et al., 2012). Gene expression was represented by tags per million. Genes expressed by less than four cells were filtered. Correlation coefficients of Aldh1a1expression with other genes were calculated using Spearman’s correlation analysis. Enrichment scores were calculated by GSEA (Subramanian et al., 2005) executed in MATLAB R2013a (MathWorks). Single-cell qPCR data with Ct values were analyzed by the SINGuLAR Analysis Toolset 3.0 software (Fluidigm). Principal components analysis, hierarchical clustering, correlation, and visualization were performed according to the toolset by using R software. Gene ontology enrichment analysis and GSEA were performed with DAVID informatics resources. Spanning-tree Progression Analysis of Density-normalized Events (SPADE) was conducted with SPADE V3.0 software (Qiu et al., 2011) (http://pengqiu.gatech.edu/software/SPADE/index.html).

H&E Staining

The mouse uteri at different developmental stages were fixed in 4% (w/v) paraformaldehyde and then dehydrated in an ethanol gradient, prior to embedding in paraffin and sectioning at 10-μm thickness. The paraffin sections were then stained with H&E.

Immunofluorescent Staining

The series of 10-mm-thick sections were rehydrated, fixed with 4% (w/v) paraformaldehyde for 30 min, rinsed three times with PBS, and treated with blocking solution (1% BSA) for 30 min, prior to incubation with primary antibodies at 4°C overnight. The primary antibodies rabbit anti-mouse monoclonal antibodies against ALDH1A1 (Abcam, ab52492), Ki67 (Abcam, ab16667), EZH2 (CST, 5246), PAX2 (Abcam, ab79389), VIMENTIN (DakoCytomation, 103465-002), and GATA2 (Abcam, ab173817), and rat anti-mouse monoclonal antibodies against CD326 (eBioscience, 17-5791) were used to detect the expression of selected proteins within the mouse uteri. Goat anti-rabbit secondary antibody (Invitrogen, A11008), goat anti-rat secondary antibody (Beyotime Institute of Biotechnology, A0507), and DAPI (Beyotime Institute of Biotechnology, C1002) were used to visualize the respective primary antibodies and the cell nuclei. All procedures were carried out according to the manufacturer's instructions.

Cell Culture

The digested single-cell suspension as well as the CD326+ epithelial cells sorted by FACS was further cultured in the epithelial defined CnT-PR medium (CELLnTEC) at 37°C and 5% CO2. To validate the results of bioinformatics GSEA and further investigate the molecular cascade controlling the self-renewal and proliferation of uterine epithelial cells, we investigated epigenetic and signaling pathways that were highly enriched at the early developmental stage. We further investigated small molecules targeting the gene ontology, which are highly enriched in the first cell cluster. The effects of these small molecules on cell proliferation of developing uterine epithelia from P7, which were cultured in epithelial cell defined medium (for a list of small molecules and concentrations see Table S8), were evaluated with absorbance readings at 450 nm, utilizing the cell counting kit-8 (CCK-8; Dojindo). The colony-forming unit (CFU) efficiency of the ALDH1A1+ epithelial cells and ALDH1A1− epithelial cells were compared by seeding at the density of 10,000, 5,000, and 1,000 cells per well of a 24-well plate after FACS, and cultured in the epithelial defined CnT-PR medium (CELLnTEC) at 37°C and 5% CO2.

Statistical Analysis

The CFU efficiency of the ALDH1A1+ cells and ALDH1A1− cells after FACS was compared by using an unpaired t test in Prism 5.0 software; normality was tested by Kolmogorov-Smirnov test, with all p values less than 0.05 being considered statistically significant. The effects of small molecules against the respective gene ontology and pathways on cell proliferation were compared using Dunnett's one-way ANOVA in Prism 5.0 software, with all p values less than 0.05 being considered statistically significant.

Author Contributions

B.W.: acquisition of data, data analysis and interpretation, manuscript writing; C.A.: data analysis and interpretation; Y. Li: acquisition of data, manuscript writing; Z.Y., L.G., Y. Liu: acquisition of data; Z.L., D.Z.: data analysis; B.H.: manuscript writing; H.O.: conception and design, manuscript writing; X.Z.: conception and design, manuscript writing.