Single-cell insights into epithelial morphogenesis in the neonatal mouse uterus.

All experiments using CD1 mice ( Mus musculus ) were approved by the Institutional Animal Care and Use Committee of the University of Missouri. For the isolation of uterine epithelial tubes, neonatal uteri were isolated and digested in 1% trypsin (Sigma, T4799) in calcium and magnesium-free (CMF) HBSS (Gibco, 14175-095) for 5 to 20 min at 37 °C on orbital shaker. The digestion was stopped with 1% soyabean trypsin inhibitor (Gibco, 17075-029) in CMF-HBSS with 5 mM MgCl2 and 0.1 mg/mL Dnase-1 (Roche, 10104159001). Pasteur pipettes were carefully forged to create a diameter that allows the uterine fragment to enter inside it, forcing the epithelial portion to squeeze out while pulling with a mouth pipette. The epithelial tube fragments were rinsed gently with 0.04% BSA CMF-HBSS solution and placed into Eppendorf tubes. Single-cell RNA sequencing was conducted at the University of Missouri Genomics Technology Core facility as described previously ( 80 ). Briefly, single-cell droplet formation with a target count of 10,000 cells per sample was performed using the 10× Chromium system (10× Genomics; Pleasanton, CA). The generated scRNA-seq libraries were sequenced on an Illumina NovaSeq 6000. The base call files (BCL) were demultiplexed, processed to generate FASTQ files with Cell Ranger v6.0.2 command mkfastq , and then aligned to the mouse reference genome GRCm39 followed by filtering and barcode and unique molecular identifier (UMI) counting to generate count matrices for each sample. The count matrices were merged, and the dataset was then subjected to the standard quality control and clustering pipeline in Seurat v4.0 ( 81 ). Briefly, cells with unique feature counts over 2,500 or under 200, as well as cells with greater than 5% mitochondrial counts, were excluded from further analyses. Data were normalized with log normalization using the default scale factor of 10,000. Genes with high variability in expression across cells were identified with the FindVariableFeatures function. Data were then scaled, and linear dimension reduction was performed using principal component analysis (PCA). Cells were clustered with FindNeighbors and FindClusters with a resolution of 0.2, then subjected to UMAP nonlinear dimensional reduction to identify cell populations. Markers for each cluster were identified with FindAllMarkers and used to distinguish cell types which were determined by literature search and PangaloDB ( 82 ). Cell populations were integrated across time points with Harmony v0.1 ( 27 ). Gene ontology analysis was performed with the DEGs across pseudotime using ShinyGO v0.76.3 ( http://bioinformatics.sdstate.edu/go/ ). Prospective cell trajectories and lineages across pseudotime were predicted with Monocle3 ( 83 ) and differentially expressed genes (DEGs) (FDR P < 0.01) across the pseudotime trajectory determined. Transcription factor binding sites (TFBS) significantly (FDR P < 0.01) enriched near DEGs across pseudotime were identified with ChEA3 (ChIP-X Enrichment Analysis Version 3) using the ENCODE database ( 35 ). Transcription factors with significantly enriched binding sites near DEGs were then examined for RNA expression across cell populations. Additional supporting materials and methods can be found online.

Cellular heterogeneity in the epithelium of the developing neonatal mouse uterus was determined by scRNA-seq on PNDs 1, 5, 12, and 15 (n = 3 biological replicates per PND) ( Fig. 1 A ). The approach was to isolate epithelial sheets from the uterus after enzymatic digestion using physical separation ( 24 ). Note the smooth LE and nascent gland buds in the epithelium sheets recovered from a PND 15 mouse uterus ( Fig. 1 A ). Subsequently, the sheets were dissociated into single cells for sequencing on the 10X Genomic Chromium System platform. High-quality transcriptome profiles were obtained for 18,308 total cells with a mean sequencing depth of 43,689 reads per cell and 13,171 unique molecular identifiers (UMI) per cell. Importantly, systematic comparisons of technical batches and biological replicates confirmed high batch-to-batch or individual-to-individual concordance. As a reference, the developing neonatal mouse uterus was evaluated by immunofluorescent localization of cytokeratin 8 (KRT8) and FOXA2 ( Fig. 1 B ). Of note, FOXA2-expressing differentiated GE cells are present on PNDs 12 and 15 but not PNDs 1 and 5. Cell heterogeneity in the epithelium of the neonatal mouse uterus revealed by single-cell RNA-seq. ( A ) Schematic of experimental design. Each uterine horn was trimmed with scissors, enzymatically digested, and sheets of epithelium removed using manual pressure and a pulled pipette. Darkfield photomicrograph of an epithelium sheet recovered from a postnatal day (PND) 15 mouse uterus (center). Epithelial cells were isolated from the sheets, and single-cell transcriptomes were generated using a 10× Genomics and Illumina sequencing approach. ( B ) Histology of the developing mouse uterus with immunolocalization of cytokeratin 8 (KRT8) and forkhead box A2 (FOXA2). Nuclei are pseudocolored white. (Scale bars, 100 μm.) ( C ) UMAP plots of scRNA-seq data colored by cell type. ( D ) Annotation was performed based on expression of known marker genes presented in the dotplots. Epi, epithelium; GE, glandular epithelium; LE, luminal epithelium; Mes, mesenchyme; Myo, myometrium; PND, postnatal day; S, stroma. To assign each profiled cell to a specific type, a combination of automated unsupervised clustering ( Dataset S1 ) and manual annotation using established marker genes was employed ( SI Appendix , Table S1 ). Very few contaminating mesenchyme (Mes) or stroma cells were identified in the samples based on vimentin ( Vim ) expression. The undetermined Epi of the PND 1 and 5 uterus contained 7 and 8 unique epithelial cell clusters, respectively ( Fig. 1 C ). The PND 12 and 15 uterus contained 5 cell clusters consisting of determined LE, as well as one (PND 15) or two (PND 12) types of GE. Proliferating cells were present in all samples in discrete cell clusters ( Fig. 1 D ). Calb1 expression delineated the specification and determination of Epi to LE as it matured from PNDs 1 to 15. Likewise, Cxcl15 and Foxa2 expression marked the differentiating and developing GE on PNDs 12 and 15. As expected, datasets from PND 12 and 15 epithelium contained considerably more LE (98.3% and 95.7%) than GE (1.7% and 4.3%) cells, respectively. Interestingly, Sox9 expression progressively transitioned from 5 Epi cell types on PNDs 1 and 5 to only one type of LE and both GE types on PND 12 but was present in only the GE on PND 15. Next, pairwise analysis of differentially expressed genes (DEGs) was performed using the scRNA-seq data ( Dataset S2 and Fig. 2 A ). That analysis found (log2FC ≥ 0.5, P ≤ 0.05) considerable changes in gene expression between PNDs 1 and 5 (1,638 increased and 2,369 decreased), PNDs 5 and 12 (2,867 increased and 2,637 decreased), and PNDs 12 and 15 (2,638 increased and 2,623 decreased) in the uterine epithelium. Gene ontology (GO) and biological pathway analysis for genes increased between PNDs 1 and 15 revealed enrichment for pathways involve in cell adhesion, metabolism, and cancer and signaling pathways involving calcium, FOXO, tumor necrosis factor (TNF), mitogen-activated protein kinase (MAPK), and PI3K-Akt. The same analysis for genes decreased revealed enrichment for pathways including ECM-receptor interactions, focal adhesion, and cell cycle and DNA replication ( Fig. 2 B ). Combinatorial analysis revealed that 278 genes increased between PNDs 1 and 15, including Calb1 , Lpar3 , and Sult1d1 that are specifically expressed in the LE of the adult mouse uterus ( 25 ) ( Fig. 2 C ). Similarly, 308 genes decreased between PNDs 1 and 15, including Lef1 that is expressed in the GE of the neonatal mouse uterus ( 26 ) ( Fig. 2 D ). Analysis of differentially expressed genes in the epithelium of the neonatal mouse uterus. ( A ) Volcano plots displaying pairwise comparison of genes increased (red) and decreased (blue) across postnatal day (PND) using scRNA-seq data. ( B ) Biological pathways enriched in differentially expressed genes increased during development. ( C ) Venn diagram of genes increased in the epithelium during development. ( D ) Venn diagram of genes decreased in the epithelium during development. Next, undifferentiated (PND 1 and 5) and specified uterine epithelial samples (PND 12 and 15) were separately integrated, harmonized, and analyzed ( 27 ) ( Fig. 3 A and Dataset S3 ). Proliferation is essential for epithelial development and maintenance. In the undetermined Epi of PND 1 and 5 uteri, many of the proliferating cell marker genes ( Ccnd1, Mcm2, Mcm5, Mki67, and Pcna ) were enriched in cluster 4 ( SI Appendix , Fig. S1 ). In the determined LE of PND 12 and 15 uteri, most of the proliferating cell marker genes were enriched in clusters 9 and 12. Prior studies of the neonatal and adult mouse uterus have defined several candidate genes and biological pathways that potentially mark stem progenitor cells in the epithelium that generate or regenerate the LE and/or GE during development, the estrous cycle, or after parturition or injury ( 28 , 29 ) ( SI Appendix , Table S1 ). In this regard, the putative markers Aldh1a1 and Ald1ha3 transcripts were most abundant in clusters 10 and 11 of PND 12 and 15 uteri ( Fig. 3 B and SI Appendix , Fig. S2 ). In situ hybridization revealed that Aldh1a1 mRNA was most abundant in the epithelia of PND 1 and 5 uteri and then transitioned to predominantly the nascent GE in PND 12 and 15 uteri ( Fig. 3 B and SI Appendix , Fig. S3 ). It is noteworthy that a few LE cells express Aldh1a1 on PND 12 and 15. Gstm7 was enriched in clusters 1 and 4 of PND 1 and 5 uteri and clusters 11 and 12 of PND 12 and 15 uteri. Top2a was enriched clusters 4 and 5 of PND 1 and 5 uteri and cluster 12 of PND 12 and 15 uteri. Integration of single-cell transcriptome data for epithelial cells from PND 1 and 5 and PND 12 and 15 uteri. ( A ) Visualization of epithelial cell types in UMAP (Uniform Manifold Approximation and Projection). Dotplots of gene expression used to annotate cell types. Epi, epithelium; GE, glandular epithelium; LE, luminal epithelium. ( B ) UMAP visualization of select genes in the integrated samples and in situ localization of mRNA in the neonatal mouse uterus. Nuclei are pseudocolored white. (Scale bars, 200 μm.) Epi, epithelium; LE, luminal epithelium; GE, glandular epithelium; Mes, mesenchyme. The LE and GE cells were then subclustered from the integrated and harmonized PND 12 and 15 samples based on Foxa2 and Cxcl15 expression that mark differentiating GE ( 30 ) ( Fig. 4 A , SI Appendix , Table S1 , and Dataset S4 ). In the LE, 5 different cell populations were identified with differences in expression of Calb1 , Pcna , Tacstd2, and many other genes. Although the total number of GE cells was low (n = 222), two different cell populations were identified with differences in expression of Foxa2 , Cxcl15 , Pcna , Aldh1a1 , and several other genes ( Fig. 4 B and Dataset S4 ). Genes increased in GE compared to LE cells included Gstm7, Ccnd1, Aldh1a1, Gas6, Aldh1a3, Gpx2, Cited4, Sox9, Fzd10, Pcna, Lpl, and Mcm5 ( Dataset S5 ). Gas6 is a gene that encodes a gamma-carboxyglutamic acid (Gla)-containing protein implicated in cell proliferation ( 31 ). Although expressed in very few cells of PND 1 and 5 uteri, Gas6 mRNA was noticeably enriched in the Foxa2 -expressing cells of PND 12 and 15 uteri ( Fig. 3 B and SI Appendix , Fig. S1 ). In situ hybridization localized Gas6 mRNA to the developing GE of PND 12 and 15 uteri. Cited4 is encoded by an intronless gene that belongs to the CITED family of transcriptional coactivators that bind to several proteins, including CREB-binding protein (CBP) and p300, as well as transcription factor AP2 (Tfap2) ( 32 ). Cited4 mRNA was observed in very few cells of the epithelium on PND 1 and 5 but increased in epithelial cells on PND 12 and 15, particularly in the developing GE ( Fig. 3 B and SI Appendix , Fig. S3 ). Of note, those genes may also be expressed in the LE as the uterus matures ( 25 , 33 ). Analysis of epithelial cell types in PND 12 and 15 uteri. Integrated data from PND 12 and 15 samples were subclustered based on Foxa2 and Cxcl15 expression. ( A ) UMAP plots of LE and GE. ( B ) Ridgeplots of gene expression in LE or GE subclusters. ( C ) Biological pathways enriched in genes differentially expressed in GE cells. ( D ) Biological pathways enriched in genes differentially expressed in LE cells. Gene ontology (GO) and biological pathway analysis ( 34 ) for the 81 genes increased (average log2FC ≤ 0.5, adjusted P -value ≤ 0.01) in the GE as compared to LE cells on PND 12 and 15 revealed enrichment for DNA replication, cell cycle, and tight junction among other pathways ( Fig. 4 C ). In contrast, GO and pathway analysis for the 131 genes increased in the LE as compared to GE cells on PND 12 and 15 indicated enrichment for apoptosis, proteoglycans in cancer, and estrogen signaling pathway ( Fig. 4 C ). Next, TF target overrepresentation analysis for the genes differentially expressed between the LE and GE was conducted using ChIP-X Enrichment Analysis Version 3 (ChEA3) ( 35 ). The TF associated with the 81 genes increased in GE cells (i.e., ESR1, FOXA2, GATA, TCF7L2) were different from those in the 131 genes increased in LE cells (i.e., STAT3, JUN, FOS, CEBPD, REST, NR3C1). Of note, FOXA2 either directly interacts with or facilitates the DNA binding of several TFs including ESR1 and TCF7L2 as well as FOXA2 itself ( 36 , 37 ). Further, several genes, including Cited4, Cxcl15, Fzd10, Foxa2, Gas6, Gpx2, Gstm7, and Ihh, that are increased in the GE as compared to the LE were also identified as FOXA2 targets in a ChIP-seq study of the PND 12 and adult pseudopregnant mouse uterus ( 19 ). Monocle3 ( 38 ) was then implemented for trajectory analysis to reveal branching trajectories in the PND 12 and 15 integrated scRNA-seq data ( Dataset S4 ). First, trajectory analysis was anchored in the proliferative cells marked by Mcm2, Mcm5, Mik67, Pcna, and Top2a expression ( Fig. 5 A ). This analysis identified a preferred chosen path of the development of LE and GE cells from a proliferative cell population ( Fig. 5 B ). The chosen path of LE to GE cells was subset and reclustered followed by trajectory analysis anchored in the proliferative cells ( Fig. 5 C ). This analysis found that the proliferative cells bifurcated into both LE and GE cells. Note the increased expression of Foxa2, Cited4, Gas6, and Aldh1a1 in GE cells and Tacstd2 in LE cells as the cells become terminally specified and differentiated ( Fig. 5 D ). Collectively, this study provides a framework for the GRNs and biological pathways that govern epithelial development and differentiation in the neonatal uterus ( Fig. 6 ). Trajectory analysis of epithelial cells from PND 12 and 15 uteri. ( A ) Pseudotime analysis was conducted using Monocle 3 with proliferative cells as an anchor start site (asterisk). ( B ) Visualization of the chosen versus unchosen lineage route from the start site using proliferative progenitor cells. ( C ) Pseudotime analysis using only cells in the chosen lineage route with Monocle 3 and proliferative cells as an anchor start site (asterisk). ( D ) UMAP plots of select gene expression. Schematic summary of cell types and molecular features of uterine epithelial development, determination, and differentiation in the neonatal mouse uterus.

An important developmental event in postnatal uterine morphogenesis is maturation and specification of the undetermined epithelium into determined LE. This study highlighted the molecular features of this event that includes the upregulation of 278 genes in the determined LE cells by PND 15 including Calb1, Gsto1, Lpar3, and Sult1d1 , which are known to be expressed specifically in the determined LE cells of the adult mouse uterus. The genes up-regulated in the uterine LE between PNDs 5 and 15 included 31 TFs ( Atf7 , Nfil3 , Phf1 , Irx3 , Hoxb5 , Zbtb4 , Hoxb1 , Nr4a1 , Nfe2l3 , Dpf1 , Jdp2 , Sox13 , Nanog , Stat4 , Esrrg , Foxj2 , Sox30 , Ovol2 , Zkscan4 , Trerf1 , Thap3 , Arntl1 , Tef , Rora , Sox5 , Zscan2 , Casz1 , Hivep3 , Ar , Klf9 , and Klf15 ). Apart from Ar , Klf9, and Klf15 ( 39 ), none of those TFs have been investigated in the mouse uterus but may serve as components of specific GRNs governing epithelial maturation and determination in the neonatal uterus. Discovery of the critical GRNs involved in uterine epithelial morphogenesis will require additional studies utilizing techniques to study chromatin remodeling and TF binding in the developing epithelium such as single-cell ATAC-seq and ChIP-seq or Cut&Run-seq. A number of different TFs ( Foxa2, Pax8, Sox9, Sox17, and Trop2 ) and other genes ( Aldh1a1, Aldh1a3, Axin2, En2, Gstm7, Hist1hao, Hmga1, Lef1, and Lgr5 ) have been implicated in mouse and/or human endometrial epithelial stem progenitor cell differentiation and/or expansion ( 18 , 40 – 45 ). In the present study, many of those genes were found to be minimally expressed in the differentiating epithelium of the uterus. In contrast, Aldh1a1 was found to be most abundant in the developing GE on PNDs 12 and 15 by in situ hybridization and scRNA-seq analyses ( Fig. 3 ). The Aldh1a1 and Aldh1a3 genes encode an enzyme involved in retinoic acid metabolism that affects the proliferation and/or early differentiation of stem progenitor cells ( 46 ) and has been implicated in adenogenesis in the neonatal mouse uterus ( 47 ). A previous study found that Aldh1a1 expression was highly correlated with proliferation markers, stem cell markers such as Lgr5 , key TFs such as Pax2 and Emx2 , and Wnt signaling genes such as Wnt7a and Ctnnb1 ( 47 ). Indeed, self-renewal capacity of neonatal mouse uterine epithelial cells was correlated to Aldh1a1 activity based on an in vitro colony formation assay. Thus, Aldh1a1 expressing cells have classic hallmarks of stem progenitor cells in the epithelium of the neonatal mouse uterus ( 48 ); however, Aldh1a1 null mice are fertile and healthy ( 49 ). The Aldh1a family of enzymes includes Aldh1a1 , Aldh1a2 , Aldh1a3, and Aldh3a1 , and Aldh1a3 was coexpressed along with Aldh1a1 in different epithelial cell types in the present study. Thus, future studies should focus on inactivating other key genes such as Aldh1a3 to understand their function in uterine adenogenesis and epithelial regeneration. One interesting TF identified in the present study is Cited4 , which belongs to the CBP/p300-interacting transactivator with glutamic acid and aspartic acid-rich tail (CITED) family which is induced by various cytokines and participates in cytokine-induced proliferation and differentiation ( 32 ). The CITED4 TF regulates gene transcription, acts as transcriptional coactivator for TFAP2/AP-2, and enhances estrogen-dependent transactivation mediated by ESR1. Of note, ESR1 is an important regulator of postnatal mouse uterine development ( 50 ) and essential to maintain differentiated glands in the uterus ( 51 ). An important paralog of the Cited4 gene is Cited2 that has established roles in trophoblast cell differentiation in the placenta ( 52 ). The precise biological role of Cited4 in the neonatal and adult mouse is not known as Cited4 null mutants are embryonic lethal. A critical event in postnatal development of the uterus is the emergence of the GE lineage from the LE that begins around PND 7 coincident with epithelium determination and specification. Adenogenesis is critically dependent on Foxa2 expression in the differentiating GE ( 18 , 19 ); thus, the molecular mechanisms controlling Foxa2 gene activation in GE cells is of paramount importance. Moreover, the biological mechanisms underpinning its biological actions in adenogenesis is not known in the developing uterus. In other systems, Foxa2 facilitates signal-dependent cell lineage initiation and bifurcation via chromatin remodeling and enhancer priming, and then, it elicits cell type-specific gene expression via recruitment of lineage-specific TFs and activation of GRNs and biological pathways governing morphogenesis ( 36 , 53 , 54 ). A particular limitation of the present study is the relatively low number of total GE cells (n = 222) profiled by single-cell transcriptomics in the PND 12 and 15 mouse uterus. Indeed, other available single-cell studies of the neonatal mouse uterus also display this limitation and others ( 47 , 55 ). Single-cell transcriptome or multiome (ATAC-Seq + RNA-Seq) analysis of cells from the entire uterus or isolated epithelium is technically feasible but will be limited by the relatively low number of GE cells in the developing mouse uterus during initial adenogenesis. A potential solution is to use unique mouse genetic models to enrich for Foxa2-expressing cells such as Foxa2 venus fusion (FVF) reporter mice ( 56 ). FVF reporter mice possess a knock-in reporter (venus fusion that is a variant of yellow fluorescent protein) expressed under endogenous Foxa2 control but does not interfere with its function as a pioneer transcription factor. The FVF reporter accurately reflects FOXA2 protein levels and distribution, allowing for the determination of cells actively expressing FOXA2 in vivo. The FVF reporter would permit FOXA2 quantification and isolation of GE cells using a fluorescence-activated cell sorter, which has recently been successfully employed to elucidate how FOXA2 regulates developmental plasticity and cell lineage during embryogenesis and differentiation of organs such as the pancreas ( 36 , 53 ). Once isolated, the GE cells could be subjected to a variety of analyses including single-cell ATAC-Seq and RNA-Seq as well as bulk RNA-seq, ChIP-seq, and Cut&Tag-seq. In addition, three-dimensional organoids developed from the isolated GE cells would be an attractive model for in vitro mechanistic studies ( 42 , 50 , 57 ). In summary, this study resolved specific cell-type populations and their individual transcriptome profiles with scRNA-seq analyses to advance our understanding of uterine epithelial maturation and differentiation during the critical periods of epithelial determination and lineage bifurcation in the developing postnatal mouse uterus. This atlas provides a unique molecular foundation for investigating how uterine epithelial cell determination, identities, and lineages are achieved during development in the uterus ( Fig. 6 ). Future mechanistic studies will need to focus on the intrinsic GRNs governing epithelial cell morphogenesis as well as the instructive and essential role of the mesenchyme and stroma in epithelial cell morphogenesis ( 50 ). This knowledge is essential, as the success of gland development is a major determinant of the ability of the adult uterus to establish and maintain pregnancy. Recurrent pregnancy loss is a major cause of female infertility and observed in sheep and mice ( 18 , 58 – 63 ) as well as humans ( 64 ) that lack glands in the endometrium of their uterus. In addition, uterine glands are involved in the most prevalent type of uterine diseases (adenomyosis, endometriosis, and adenocarcinoma) that negatively affect fertility and health of women. Thus, an increased understanding of uterine epithelium biology is important for diagnosis, prevention, and treatment of fertility and pregnancy problems as well as endometrial-based diseases in women ( 65 – 68 ). This knowledge is also key for the development of novel medical therapies aimed at the endometrium to treat infertility and endometrial pathologies in women ( 69 ), such as Asherman’s syndrome (thin endometrium) ( 64 , 70 – 72 ). Excessive or ectopic growth of the glands is a facet of endometrial cancer, polyps, adenomyosis, and endometriosis ( 73 – 78 ). Consequently, elucidating the cellular and molecular mechanisms governing endometrial epithelial development is of paramount importance to understand and treat infertility and disease with the goal of improving women’s reproductive health ( 28 , 29 , 79 ).

Appendix 01 (PDF) Click here for additional data file. Dataset S01 (XLSX) Click here for additional data file. Dataset S02 (XLSX) Click here for additional data file. Dataset S03 (XLSX) Click here for additional data file. Dataset S04 (XLSX) Click here for additional data file. Dataset S05 (XLSX) Click here for additional data file. Dataset S06 (XLSX) Click here for additional data file.