Single cell RNA-seq study of wild type and Hox9,10,11 mutant developing uterus

The uterus is a remarkable organ that must guard against infections while maintaining the ability to support growth of a fetus without rejection. The Hoxa10 and Hoxa11 genes have previously been shown to play essential roles in uterus development and function. In this report we show that the Hoxa9,10,11, Hoxc9,10,11, Hoxd9,10,11 genes play a redundant role in the formation of uterine glands. In addition, we use single cell RNA-seq to create a high resolution gene expression atlas of the developing wild type mouse uterus. Cell types and subtypes are defined, for example dividing endothelial cells into arterial, venous, capillary, and lymphatic, while epithelial cells separate into luminal and glandular subtypes. Further, a surprising heterogeneity of stromal and myocyte cell types are identified. Transcription factor codes and ligand/receptor interactions are characterized. We also used single cell RNA-seq to globally define the altered gene expression patterns in all developing uterus cell types for two Hox mutants, with 8 or 9 mutant Hox genes. The mutants show a striking disruption of Wnt signaling as well as the Cxcl12/Cxcr4 ligand/receptor axis.

glands showed defective development post pregnancy 22 . Mice with multiple combinations of mutations of either the Hox10 or Hox11 paralog groups have interesting kidney and limb abnormalities, but no described infertility phenotypes extending beyond those associated with the Hoxa10 and Hoxa11 genes 23,24 . These results suggest unique roles for Hoxa10 and Hoxa11 in uterus development and function.
In this report we extend this approach to search for possible female fertility functions for the Hoxc9,10,11 genes. We observed that while Hoxa9,10,11+/−, Hoxd9,10,11+/− mice showed reduced fertility, with about half normal litter size, mice that were in addition heterozygous for mutation of the Hoxc9,10,11 genes were almost completely infertile. In this report we show that Hoxc9, 10,11 genes have redundant function with Hoxa9, 10,11 and Hoxd9,10,11 in uterine gland formation.
Single cell RNA-seq (scRNA-seq) is a powerful tool for the dissection of normal and mutant development 26 . It can define the global gene expression states of the multiple cell types present in a developing organ. Analysis using scRNA-seq can help determine how early lineage decisions are made 27 . It can characterize the transcription factor codes that define different cell types, and it can provide a comprehensive analysis of potential ligand-receptor interactions 28,29 . In this report we used scRNA-seq to examine the wild type developing uterus at PND12, as early gland formation is taking place. In addition, we used scRNA-seq to examine the perturbed gene expression patterns of all cell types of the PND12 Hoxa9,10,11+/−, Hoxc9,10,11+/−, Hoxd9,10,11+/− (ACD+/−) mutant uteri, with heterozygous mutation of 9 Hox genes, as well as Hoxa9,10,+/−, Hoxc9,10,11+/−, Hoxd9,10,11+/− (ACD+/−WTA11), with two wild type Hoxa11 genes and heterozygous mutation of 8 Hox genes. The results create a single cell resolution atlas of the gene expression patterns of the wild type developing uterus and also define the changed gene expression profiles of all cell types in mutant uteri. We observed a compelling disruption of Wnt signaling in the Hox mutants, as well as strongly reduced stromal expression of Cxcl12 (SDF1), which encodes a ligand for Cxcr4, expressed by the LE. In multiple developing contexts Cxcr4/Cxcl12 interaction has been previously shown to drive cell migration and tubulogenesis [30][31][32][33][34] , key steps in uterine gland formation.
Hoxc10 and Hoxc11 are also expressed in the human female uterus, suggesting possible fertility function 35 . Frame shift mutations were simultaneously introduced into the first exons of Hoxc9, Hoxc10, and Hoxc11 using recombineering 36 . It is important to note that this approach leaves shared enhancer regions intact, thereby reducing impact on the expression of the remaining Hox genes not mutated. Female Hoxc9,10,11−/− mice mated to wild-type males had similar size litters (9.33 ± 0.88 pups/litter, n = 3) to those from wild-type matings (9.92 ± 0.35, n = 13). All Hox mutant female fertility assays were carried out through mating to wild type males. To further test for possible Hoxc9,10,11 fertility function we then used a "sensitized" mouse that was heterozygous for the six Hoxa9,10,11+/− Hoxd9,10,11+/− (AD+/−) genes, which results in reduced fertility with 5.3 pups per litter 25 . When we made mice heterozygous for all nine Hoxa9, 10,11,Hoxc9,10,11,Hoxd9,10,11 (ACD+/−) genes we observed a dramatic decrease in female fertility. From a total of 41 vaginal plugs only two litters were born, with a total of 6 pups, giving about 0.15 pups per vaginal plug. The further reduction of fertility in the ACD+/− mice compared to AD+/− gives evidence for fertility function for the Hoxc9,10,11 genes.
We also generated ACD+/− mice but with two wild type Hoxa11 alleles (Hoxa9,10+/−, c9,10,11+/−, d9,10,11+/−) (ACD+/−WTA11). As noted previously, Hoxa11 is known to play a critical role in both male and female fertility 19,20,37 . Four Hoxa9,10+/−, c9,10,11+/−, d9,1011+/− females were mated with wild type males, giving twelve litters totaling 125 pups, or 10.41 ± 0.77 pups per litter, which is comparable to wild type. It is remarkable that the addition of one wild type Hoxa11 allele to the ACD+/− mice, with nine mutant Hox alleles, can rescue a wild type phenotype. This shows the relative importance of Hoxa11 in female fertility. It should be noted that the most severe Hox mutant genotype resulting from any of the Hox mutant female fertility testing matings would be ACD+/−, since a wild type allele would be inherited for each Hox gene from the father, and ACD+/− mice are fully viable. Hence, progeny genotypes should not impact the fertility testing results.
In this study we focus on the ACD+/− genotype, which shows a highly penetrant female infertility, yet the mice are quite healthy, without the severe kidney malformations associated with removal of additional Hox10,11 genes. ACD+/− mice show a dramatic reduction in the number of uterine glands. Histologic analysis of the ACD+/− uteri showed a reduction in the number of uterine glands ( Fig. 1D-F) compared to WT uteri ( Fig. 1A-C). At PND10, PND20 and adult the luminal epithelium (LE), mesenchymal stroma, and the muscle layers of ACD+/− mice appeared normal ( Fig. 1). At all stages examined, however, there was a dramatic reduction in the number of uterine glands. The most striking difference was seen in the adult uterus ( Fig. 1). Glands were counted from transverse sections of uteri. Wild-type adult uteri contained 32.63 ± 2.46 glands/section (n = 13), while ACD+/− uteri contained 3.44 ± 0.66 glands/section (n = 7). These results give evidence for the sub-fertility being the result of a drastic reduction in the number of uterine glands, which supply vital proteins, including Lif 38 and calcitonin 3 , which are required for embryo implantation. Interestingly, the adult ACD+/−WTA11 mice, with two wild type Hoxa11 alleles, had wild type numbers of uterine glands (37.0 ± 6.66 glands/section, n = 3), again showing the importance of normal Hoxa11 expression in female fertility.
The ACD+/− mice produce normal numbers of embryos. WT and ACD+/− mutant females were mated with WT males. The detection of a vaginal plug equaled gestational day 1 (GD1). Mice were euthanized at GD4 and one horn of each uterus was gently flushed. Similar numbers of blastocysts were isolated from wild-type (4.45 ± 0.72 blastocysts/uterine horn) and ACD+/− mice (3.38 ± 0.54 blastocysts/uterine horn). Recovered blastocysts from ACD+/− mice were morphologically normal (data not shown).
Single cell RNA-Seq analysis of the wild type and ACD+/− and ACD+/−WTA11 mutant developing uteri. scRNA-seq is a powerful technology that can give a high resolution definition of the gene expression programs that drive development 26,27,[39][40][41][42] . To better understand the normal and perturbed gene expression patterns of the developing uterus we carried out scRNA-seq on the WT, ACD+/−, and ACD+/−WTA11 uteri at PND12. We were particularly interested in uterine gland formation and how this process might be disrupted in the ACD+/− mutants. Gland formation begins with the budding of the luminal epithelium at approximately PND6 and continues into adulthood during endometrial regeneration following menstruation in humans and following birth in both humans and mice 4,5,8,43,44 . Wild type. scRNA-Seq was performed using Drop-Seq 45 . A total of 6,343 WT cells passed quality control filters. Unsupervised clustering of cells was carried out using the Seurat R package as previously described 46 . A t-SNE plot showing the resulting clusters is shown in Fig. 2A. Major cell types were defined using the expression of known cell type-specific genes, with epithelial cells (n = 1487) expressing Epcam (epithelial cell adhesion molecule) 47 and Cdh1, endothelial cells (n = 373) expressing Pecam1 (platelet and endothelial cell adhesion molecule) 48 and Emcn, stromal cells (n = 2241) expressing Pdgfra 49 , Dcn, and Col15a1, myocytes (n = 1749) expressing Mef2a (myocyte enhancer factor 2A) 50 and Pdlim3, lymphocytes (n = 106) expressing Cd52 51 , and Ptprc, myeloid cells (n = 39) expressing Lyz2 52 , Pf4 53 , and C1qc, pericytes (n = 269) expressing Rgs5 54 and Cspg4, and mesothelial cells (n = 79) expressing Upk3b 55 and Lrrn4. A heatmap of the top 15 marker genes from each cluster is shown in Fig. 2B. Lists of genes with compartment enriched expression are presented in Table S1.   Table S2. These results define a surprising number of distinct cell types present in the developing uterus, including four distinct categories of stromal cells and four different myocytes.
Genes with elevated expression in GE versus LE included Foxa2, Wif1, Gas6, Ccnd1, Cxcl15, BC048679, F5, and Egfl6, which were previously shown to have elevated GE expression in a microarray study 56 . Genes more strongly expressed in LE cells compared to GE included Tmsb4x, Wfdc2, 1600014C10Rik, Calb1, S100a6 and Ctsd. Mesothelial cells (n = 79) are specialized squamous cells that form a monolayer that lines the body's serous cavities and internal organs, providing a slippery, non-adhesive and protective surface. Genes with elevated expression in mesothelial cells included Upk1b, Upk3b, Cldn15, Fgf1, and Fgf9.
Four endothelial cell sub-clusters were identified with all clusters expressing the pan-endothelial markers Pecam1 and Emcn. The vascular endothelial markers Cd34 and Sox17 were expressed in every sub-cluster except EndoLym (lymphatic). The Endo1 cell cluster expressed Pcdh17, Rnf213, Sik3, as well as Gpihbp1, which is expressed exclusively in capillaries 61 . The Endo2 cell cluster showed elevated expression of Nrp2, Lepr, Slc16a9, Penk, and Fendrr, indicating venous cell identity 62 . As expected the lymphatic EndoLym cell cluster expressed Prox1 and Lyve1 as well as Ccl21a, Ccl21b, Cldn11, and Fxyd6, while the arterial Endo3 cells expressed Cxcl12, Gja4, Jag1, and Dkk2. Both Gja4 63 , and Jag1 64 expression have been linked to endothelial arterial cell specific differentiation.
Analysis of growth factor-receptor interactions. The scRNA-seq data provides a global view of the gene expression patterns of the multiple cell types of the developing uterus. It defines the transcription factor codes that drive the distinct lineage differentiation programs. Further, it describes the combinations of growth factors and growth factor receptors expressed by each cell type. We were particularly interested in the growth factor/receptor interactions taking place between the stromal cells and luminal epithelial cells that might play a role in promoting uterine endometrial gland formation.
Most of the Hox genes targeted in the ACD+/− mutants showed strongest expression in the stromal cells at PND12, with none showing significant expression in the epithelial compartment (Fig. 5). Indeed, for the HoxA cluster only the Hoxa9,10,11 genes showed significant expression in the developing PND12 uterus, with each expressed in the stromal, pericyte and myocyte compartments (Fig. 5). For the HoxB cluster none of the genes showed strong expression in stromal cells, but many showed significant expression in the epithelial and endothelial compartments. Of particular interest, none of the HoxC cluster genes showed strong expression in the developing uterus, although Hoxc9, Hoxc10 and Hoxc11 were upregulated in the ACD+/− mutant uterus (see below). There is however low level, but perhaps functionally important, expression of Hoxc11 in the stroma (Fig. 5). Also, www.nature.com/scientificreports www.nature.com/scientificreports/ the Hoxc10,11 genes are expressed in the adult human uterus, where uterine gland formation continues, giving evidence for normal uterine function for these genes 35 . For the HoxD cluster the Hoxd3, 8,9,10,11 genes showed expression in the PND12 uterine stromal, pericyte, myocyte and endothelial compartments (Fig. 5). The Hox codes for uterine stroma, pericytes and myocytes are remarkably similar.
This stromal expression of the ACD+/− mutated Hox genes, and absence of their expression in LE, suggests that the reduced invagination of LE cells to form glands in the ACD+/− mutants could be the result of altered cross talk between stroma and LE. Therefore, we restricted the analysis to growth factors expressed in stromal cells and their corresponding receptors in any cell type. Growth factors and receptors were filtered requiring detected expression in greater than 7.5% of cells in a cluster. This was to reduce the possible interactions to a more manageable number but resulted in a list that was not comprehensive. A Circos plot was generated that displays resulting potential growth factor/receptor interactions (Fig. 6) 65 . Of particular interest, stromal cells were found to express Bmp, Cxcl12, Gdf10, Igf1, Manf, Rspo3 and Tgfb2/3, all potentially able to bind to their corresponding receptors on epithelial cells.
The stromal growth factors and their corresponding epithelial receptors were then used to generate a network plot using GeneMania software 66 (Fig. 7). The functional GO enrichment analysis indicates that some of the genes within the network are involved in gland development (P = 3.98E-16), morphogenesis of branching epithelium (P = 3.54E-16), morphogenesis of a branching structure (P = 7.67E-16), branching morphogenesis of an epithelial tube (P = 1.46E-15), and epithelial tube morphogenesis (P = 3.21E-12), consistent with a role in the initiation and/or progression of uterine gland development.

scRNA-Seq comparison of WT with mutant ACD+/− and ACD+/−WTA11 uterine cells.
To better understand the altered gene expression patterns in each mutant cell type we also carried out scRNA-seq for ACD+/− and ACD+/−WTA11 PND12 uteri. All cells of the three genotypes were combined and clustered in an unsupervised manner resulting in 11 major cell types (Fig. 8, right panel). The clustered cells could then be color coded according to genotype (Fig. 8, left panel). Of interest, there were striking genotype specific differences for three clusters, the stromal, myocyte and epithelial cells. For example, within the stromal cluster the cells of the three genotypes clearly subclustered separately. The WT and ACD+/− cells were at the two ends of the cluster, with the ACD+/−WTA11 cells in between. For the myocyte cluster there was a similar result, with the WT and ACD+/− cells at the two ends of the cluster and once again the ACD+/−WTA11 cells in between. For the epithelial cell cluster there was also genotype dependent cell separation, although not as complete as observed in the stromal and myocyte clusters. It is interesting to note that the cell types with the greatest genotype specific differences (stromal and myocyte) were the cells with the strongest expression for the mutated Hox genes.
We then specifically examined the differential expression of WT versus ACD+/− stromal cells to find possible changes in gene expression that might result in perturbed stromal-LE signaling. Genes involved in signaling with reduced expression in the mutant stroma included Thbs1, Thy1, Igfbp6, Gas6, Cd44, Ndp, and Cxcl12. Thbs1 is a major activator of Tgf-beta 67 , so reduced Thbs1 expression in the mutant stroma suggests reduced Tgf-beta signaling. As noted earlier, Tgf-beta is expressed by stromal cells and its receptor is expressed by LE. Also of interest, Thy1 can interact with integrins to drive cell adhesion/migration responses 68 , which could play a role in www.nature.com/scientificreports www.nature.com/scientificreports/ GE formation. Of note, Gas6 mutant mice are viable and give normal size litters 69 , arguing against a key role for Gas6 in GE formation.
Of particular interest, the Ndp gene with reduced expression in mutant stroma encodes Norrin protein, which can interact with frizzled receptor to drive canonical Wnt signaling 70 , which in turn is critically important for GE formation. Consistent with reduced Wnt signaling in the mutant there was increased expression of Dkk2 (1.9 FC), which encodes a Wnt inhibitor. The mutants, however, also showed increased expression of Rspo1 and Rspo3, which activate canonical Wnt signaling. In addition, we observed decreased expression of Gli1 and Ptch2 in mutants, suggesting reduced Hedgehog signaling.

RT-PCR analysis of gene expression in mutant uteri.
We carried out RT-PCR to further validate selected gene expression differences in mutants. First, we were interested in the altered expression of the Hoxc9,10,11 genes themselves. As noted previously, in the wild type uterus Hoxc9 and Hoxc10 show very low expression, with Hoxc11 www.nature.com/scientificreports www.nature.com/scientificreports/ having low but significant expression in the uterine stroma (Fig. 5). The scRNA-seq data showed significant upregulation of all three genes in the ACD+/− uterus (data not shown). RT-PCR confirmed the upregulated expression of these genes (Fig. 9). These genes show upregulation even though each is heterozygous for a frameshift mutation in the first exon that would be expected to result in increased nonsense mediated decay. These results suggest that the Hoxc9,10,11 genes normally play a limited role in early uterus development but undergo increased compensatory expression when the dose of Hoxa9,10,11 and Hoxd9,10,11 genes is reduced. It is interesting to note that in the ACD+/−WTA11 uterus, with two wild type copies of Hoxa11, there was no significant upregulation of Hoxc9,10,11 (data not shown). As noted previously, they (Hoxc10,11) are also normally expressed in the adult uterus, where uterine gland formation continues 35 .  www.nature.com/scientificreports www.nature.com/scientificreports/ RT-PCR was used to confirm the reduced expression of Lef1, Sox9, Sp5, Gli1, Ptch2, and Cxcl15 in ACD+/− mutant uteri compared to WT. We examined both WT and ACD+/− PND10 and PND20 uteri. The Lef1 transcription factor is a mediator of Wnt signaling but is also a key target, with Wnt signaling reporter transgenes often including concatemerized Lef binding sites 74 . The observed significant reduction of Lef1 expression therefore gives evidence for reduced Wnt signaling in the ACD+/− mutant uterus (Fig. 10). The observed reduced expression of Sp5, a known Wnt target 75 , is also consistent with reduced Wnt signaling (Fig. 10).
The scRNA-seq and RT-PCR data also provide evidence for reduced hedgehog signaling. Gli1 and Ptch2 encode integral components of the hedgehog signaling pathway and are also targets of hedgehog signaling 76 . Expression of Gli and Ptch are increased by hedgehog signaling and thus they serve of markers of hedgehog signaling 77,78 . The RT-PCR validation of reduced expression of both Gli1 and Ptch2 in the mutant uterus therefore gives evidence for decreased hedgehog signaling (Fig. 10). In addition, RT-PCR showed lowered expression levels for Sox9 and Cxcl15. The Cxcl15 gene shows restricted expression in the GE of the uterus (Table S1) 79 , and its reduced expression in ACD+/− uteri is therefore likely a simple reflection of reduced GE representation in the mutants.
Immunhistochemistry. We further examined the reduced Wnt signaling in ACD+/− uteri using immunostaining with Lef1 antibody. At PND10 there was dramatic reduction of Lef1 expression in ACD+/− mutants (Fig. 11, top panels). We then asked if this reduced expression persisted in adults. We observed that during diestrus, estrus, and also in the E5 post vaginal plug uterus there was significant reduction of Lef1 expression in mutants (Fig. 11). In WT adult uteri Lef1 expression was restricted to stroma during diestrus and at E5, but during estrus Lef1 was expressed in both stroma and GE. Lef1 was also expressed in budding LE that was forming GE at PND10. In mutants the very low levels of Lef1 expression observed appeared restricted to stroma.

Discussion
In this study we identify novel female fertility function for the Hoxc9,10,11 genes, with ACD+/− mice showing a much more severe infertility phenotype than AD+/− mice. Histologic analysis showed an approximate ten fold reduction in the number of uterine glands in ACD+/− mice compared to wild type. Single cell RNA-seq was used to explore the gene expression programs driving distinct lineage trajectories in the wild type developing uterus, as well as to globally define the perturbed gene expression patterns in the multiple cell types of the ACD+/− and ACD+/−WTA11 mutant uteri. Unsupervised clustering of the WT cells identified eight distinct cell types, including epithelial, endothelial, stromal, myocytes, lymphocytes, myeloid, pericytes and mesothelial. Subclustering was carried out for five of these cell types, for example dividing endothelial cells into capillary, arterial, venous and lymphatic cell subtypes. Epithelial cells were divided into glandular and luminal subtypes. In addition, we identified a considerable heterogeneity of stromal and myocyte cells. The results define Figure 10. RT-PCR analysis of altered gene expression in ACD+/− mutants. RNA was isolated from PND10 and PND20 wild type and mutant whole uteri. There was reduced mutant expression of Lef1, consistent with reduced Wnt signaling in mutants. Sp5, and target of Wnt signaling, also showed reduced expression in mutants. Gli1 and Ptch2 also showed lower expression in mutants, suggesting perturbed hedgehog signaling. Cxcl15 and Sox9 have highest expression in the GE, and their reduced expression therefore likely reflects reduced GE representation in mutants. Values were scaled to Gapdh.
www.nature.com/scientificreports www.nature.com/scientificreports/ the differential gene expression of these various cell types and subtypes and characterize the transcription factor codes that drive lineage specific development. The scRNA-seq data further allows the analysis of potential growth factor/receptor interactions between the different cell types in the PND12 uterus. These results confirm and significantly extend the previous laser capture microdissection/microarray study 56 and scRNA-seq analysis of developing uterine epithelial cells 80 .
As might be expected, the ACD+/− and ACD+/−WTA11 mutant uteri show the greatest gene expression changes where the mutated Hox genes are most strongly expressed, in myocyte and stromal cells. For these cell types the ACD+/−WTA11 cells show a striking intermediate gene expression state, between the ACD+/− and WT cells, consistent with their intermediate severity genotype (Fig. 9). There were also pronounced differences in gene expression of the ACD+/− mutant epithelial cells, even though they show very little expression of the mutated Hox genes, suggesting that this might be the result of disrupted cross talk with the mutant stroma. www.nature.com/scientificreports www.nature.com/scientificreports/ Several growth factor/receptor pathways of interest were altered in the mutants. There was reduced Gli1 and Ptch2 expression in the ACD+/− uterus (Fig. 11), giving evidence for reduced hedgehog signaling. Both Gli1 and Ptch2 were previously shown to be expressed in neonatal uterine stroma 81 . There was also evidence of reduced Tgfbeta signaling in ACD+/− uteri. The mutant stromal cells show 2.5 FC reduced expression of thrombospondin (Thbs1), which functions to activate latent Tgfbeta 82 and the epithelial cells show reduced expression (2.5 FC) of Dcn, which binds Tgfbeta and enhances its activity 83 .
Wnt signaling was dramatically reduced in the ACD+/− mutant uteri. Lef1 is a member the Tcf/Lef family of transcription factors and interacts with β-catenin to up-regulate Wnt target genes. Lef1 showed expression in a subset of stromal cells of WT PND10 uteri, especially in the mesometrial region (Fig. 11). Reduced Lef1 expression in ACD+/− mutants was shown by scRNA-seq, RT-PCR (Fig. 10) and IHC (Fig. 11). The ACD+/− mutant uteri stroma also showed reduced expression of Ndp (2.7 FC), which encodes Norrin, a non-Wnt ligand for Frizzled receptor, capable of activating canonical Wnt signaling 84 . The mutants also showed lower expression of Sp5 (Fig. 10) a known Wnt target 75 . A majority of Lef1−/− mice die by two weeks of age with missing teeth, mammary glands, whiskers, hair 20 and most relevant, absence of uterine glands 85 . Previous studies have shown that Wnt signaling is essential for proper uterine gland formation. Wnt4 86 , Wnt7a 87 , Wnt5a 88 , Wnt11 1 and β-catenin [89][90][91] are all required for gland formation. The Wnt7a mutant is of particular interest. Wnt7a expression is in the luminal uterine epithelium postnatally 87,92 , with mutants showing near complete absence of uterine glands, loss of postnatal expression of Hoxa10 and Hoxa11, and partial homeotic transformation of uterus, surprisingly, towards vagina, as measured by the mutant stratified luminal epithelium, similar to the Hox9,10,11 mutants we previously described 25 . Further, Hoxa11 is required to maintain Wnt7a expression 88 . These combined observations suggest a positive feedback loop, with Hox genes activating Wnt signaling and Wnt signaling activating Hox genes. Of note, the ACD+/− uterus is not homozygous mutant for any Hox gene. The dramatic effect on Wnt signaling resulting from heterozygous mutation of nine related Abd-B Hox genes provides further evidence of their functional overlap. Taken together these results strongly suggest that reduced Wnt signaling in ACD+/− mutants contributes to the failure of formation of normal numbers of uterine glands.
Of interest, it should be noted that the PND12 scRNA-seq data detected minimal Hoxc9,10,11 expression in the developing wild type uterus, although there was significantly upregulation in the context of reduced Hoxa9,10,11/Hoxd9,10,11 expression. The redundant Hoxc9,10,11 may therefore be normally hidden. The data do not exclude, however, the possibility that Hoxc9,10,11 genes show expression in wild type uterus at other times. The strongly synergistic phenotype, with Hoxc9,10,11 female mice showing normal fertility, and AD+/− females showing approximate 50% reduced fertility, while ACD+/− mice show complete infertility, certainly argues in favor of redundancy.
The ACD+/− mutant uteri additionally showed reduced (3.8 FC) stromal expression of Cxcl12, also known as stromal cell-derived factor 1 (Sdf1), which encodes a ligand for Cxcr4, expressed by the LE. In multiple developing contexts Cxcr4/Cxcl12 interaction has been previously shown to drive cell migration and tubulogenesis [30][31][32][33][34] , key steps in uterine gland formation. Cxcl12 is a powerful chemoattractant for many cell types [93][94][95] , and inward movement of LE cells is a key first step in uterine gland formation. Indeed, Cxcl12 and its G protein-coupled receptor Cxcr4 have been proposed to control cell movement in many developing organs 96 . For example, in the developing pancreas, similar to the uterus, Cxcl12 is expressed in the mesenchyme while Cxcr4 is expressed in epithelial cells. Either genetic mutation or pharmacologic inhibition of Cxcl12 with AMD3100 resulted in a significant reduction of epithelial buds during pancreas branching morphogenesis 97 . Similar, in the developing kidney the inhibition of Cxcl12/Cxcr4 with AMD3100 causes a severe reduction of ureteric bud branching morphogenesis 33 . Taken together these results give evidence for a role for the reduced mutant expression of Cxcl12 in the observed formation of fewer uterine glands.

Conclusions
In summary this study identifies potentially redundant roles for the Hoxa9, 10,11,Hoxc9,10,11,Hoxd9,10,11 genes in uterine gland formation. We further used scRNA-seq to define the gene expression patterns of the multiple cell types of the WT developing uterus. The results define the gene expression patterns driving lineage specific development. In addition, scRNA-seq was used to characterize the perturbed gene expression levels of all developing uterus cell types in the ACD+/− and ACD+/−WTA11 mutants. Particularly striking was the observed reduced Wnt signaling and disruption of the Cxcl12/Cxcr4 axis in the mutants.

Materials and Methods
Animals. All animal procedures were approved by the Institutional Animal Care and Use Committee at Cincinnati Children's Hospital Medical Center in accordance to their guidelines for the care and use of laboratory animals (IACUC protocol 2015-0065). The Hox mutant mice used in this study were generated by recombineering as previously described 25,36 . The Hoxa9,10 mutant mice were made exactly as previously described for Hoxa9,10,11, using the same BAC targeting vector, but the recombination event inserted only the two Hoxa9,10 mutant genes instead of all three.
The female reproductive tracts were isolated immediately after euthanasia using IACUCC approved methods, washed in cold 1X PBS, and processed for immunohistochemistry or the isolation of RNA using standard methodologies. For blastocyst isolation, gestational day 4 (vaginal plug = Day 1) uterine horns were gently flushed with 1 ml of 1X PBS and blastocysts counted using a dissecting microscope. For all experiments a minimum of 3 animals were analyzed per genotype.
Hematoxylin and eosin (H&E) staining. Female reproductive tracts were isolated as described above, then fixed in 4% paraformaldehyde, rinsed in 1X PBS and dehydrated through ethanol washes. Uterine horns were cut into four pieces (ovary/oviduct and three uterine sections) and all pieces of tissue embedded into one www.nature.com/scientificreports www.nature.com/scientificreports/ paraffin block. Glands were counted in each transverse section from a minimum of three non-contiguous slides and the average number of glands per section determined.
Immunohistochemistry. Briefly, fixed uterine horns were sectioned (5 μm), sections were mounted on slides, baked, deparaffinized in xylene, and rehydrated in a graded series of ethanol. Antigen retrieval was conducted using 10 mM sodium citrate (pH 6.0) with boiling. Sections were blocked using 4% normal donkey serum in 1X PBST (pH 7.4) and incubated overnight at 4 °C with a rabbit monoclonal anti-Lef1 (C12A5) antibody (Cell Signaling Technology) at dilution of 1:2000. Sections were washed in 1X PBST and incubated with biotinylated goat anti-rabbit secondary antibody (Vector Laboratories) at a dilution of 1:200 for 30 minutes at room temperature. After washes in 1X PBST, the immunoreactive protein was visualized using the Vectastain ABC kit (Vector Laboratories) and 3,3′ diaminobenzidine tetrahydrochloride hydrate (Sigma) as the chromagen. Sections were counterstained with nuclear fast red before coverslips were applied using Permount.
RNA extraction and real-time PCR. RNA was isolated from uterine samples using the RNeasy Plus mini kit (Qiagen) after homogenizing the samples using Qiagen's Tissue Lyser II. The quantity and purity of each RNA was determined using an Agilent Bioanalyzer nanochip. Total RNA (1 μg) was reverse transcribed to cDNA using the SuperScript Vilo kit (Invitrogen). The cDNA was ethanol precipitated and resuspended in nuclease-free water. cDNA concentration was determined using a NanoDrop Spectrophotometer and samples were stored at −20 °C until use.
Real-time PCR was conducted using Applied Biosystems TaqMan Gene Expression Assays, TaqMan Gene Expression Master Mix, nuclease free water, and the Step One Plus Real Time PCR system (Thermo Fisher Scientific). Gapdh was used as the reference gene. All samples were analyzed in triplicate for each gene tested. Statistical analyses were performed using GraphPad Prism (Version 7).

Cold protease uterine dissociation protocol for Drop-seq analysis. For Bacillus licheniformis
enzyme mix 1 ml contained 890 μl of DPBS with no added Ca or Mg, 5 μl of 5 mM CaCl 2 , 5 μl of DNase (Applichem A3778, 125 U), and 100 μl Bacillus licheniformis enzyme, 100 mg/ml (P5380, Aldrich Sigma). We euthanized post-natal day 12 (PND12) mice (morning of birth = post natal day 0, PND0) and isolated uterine horns without ovary and oviduct. We incubated tissue in 1 ml of enzyme mix at 6 °C for 2-4 minutes, then began to dissociate uterus by drawing tissue through an 18 gauge needle attached to a 3 ml syringe. Repeated 6 °C incubation and passed tissue through progressively smaller gauge needles (23 and 21 gauge) until a single cell suspension was achieved. Progress was monitored with a microscope. When single cell status was achieved, we added an equal volume of ice cold 10% FBS in 1X PBS to inactivate protease. Filtered cell suspension through a 30 μm MACS filter. Rinsed filter with 3-5 mls of ice cold 0.01% BSA in 1X DPBS into a 50 ml centrifuge tube. Transferred cells to a 15 ml centrifuge tube and pelleted cells by centrifuging at 160 g for 5 min at 4 °C. Repeated washing and centrifugation steps twice more. Resuspended cells in 1-2 mls of 0.01% BSA in 1X DPBS and determine cell viability using Trypan blue. A hemocytometer was used to determine the final concentration of cells.
Cell-type clustering and marker genes were identified using Seurat 2.3.2 46 and R 3.4.3. Data were prefiltered at both the cell and gene level with the removal of cells with low library complexity (<500 expressed genes) as well as those with a high percentage (>20%) of unique molecular identifiers (UMIs) mapping to mitochondrial genes. In addition, cells with more than 8% histone gene expression and more than 0.4% hemoglobin gene expression were also excluded. Genes were selected that had >0 UMI in more than five cells. Following filtering, the number of cells, the average number of genes called expressed per cell, and average number of transcripts (UMI) per cell respectively were, WT (6343, 1222, 2283), ACD+/− (2134, 1258, 2865), and ACD+/−WTA11 (3378, 1297, 2531). Genes with the highest expression variability were used for the principal component analysis. Significant principal components were determined using the PCElbowPlot function. Data were visually represented after clustering using Seurat's t-SNE implementation. The FindAllMarkers function was used to determine the top marker genes for each cluster with a minimum of 10% of cells expressing the gene within the cluster and a minimum logFC threshold of 0.25. Sub clustering was performed by repeating the above analysis on individual clusters.
Differential expression between clusters was determined using the FindMarkers function in Seurat with pseudocount set to 0.01. Similar cell types from the WT experiment were compared to each of the two mutant mice.
The top differentially expressed genes were filtered on genes that had a log2 fold-change greater than 1 and were expressed in at least 10% of the cells for that cluster.
Data has been deposited in GEO under accession number GSE118180. All scRNA-seq data is freely available through GEO. All other data is presented in this paper.
Circos plot. The Circos plot was generated by identifying the growth factors and receptors that are expressed in greater than 7.5% of the cells within any given cluster. The possible interactions between the stromal growth factors and their corresponding receptors in the multiple cell types were plotted. The Circos plot was generated as previously described 65 . Growth factor and receptor analysis. The growth factors in the stromal groups and their corresponding receptors in the epithelial compartments were input into GeneMania to generate a network plot as previously described 66 .