Recombineering-based dissection of flanking and paralogous Hox gene functions in mouse reproductive tracts

Hox genes are key regulators of anterior-to-posterior patterning in the developing embryo. They are organized in four clusters (HoxA,B,C,D) in mammals and are expressed in a spatiotemporal order in the embryo that matches their 3′-to-5′ arrangement on the chromosomes.

The study of the 39 mammalian Hox genes is greatly confounded by their overlapping functions. Paralogous genes on different Hox clusters are evolutionarily derived from one ancestral Hox gene. They display similar expression patterns during development, and at least in some cases their coding sequences are functionally interchangeable (Greer et al., 2000). As a result, they show a high degree of functional redundancy. For example, mutation of either Hoxa11 or Hoxd11 results in relatively minor developmental defects (Small and Potter, 1993; Davis and Capecchi, 1994), whereas mutation of both confers dramatic patterning defects in kidneys and forelimbs (Davis et al., 1995).

There is also evidence for functional redundancy of flanking Hox genes. This is particularly relevant to the AbdB class of genes at the 5′ end of the clusters, where a single ancestral gene has given rise to five Hox paralogous groups (Hox9-13). In the developing kidney (Patterson and Potter, 2004), skeleton (Haack and Gruss, 1993) and uterus (Gendron et al., 1997; Ma et al., 1998), flanking Hox genes often show similar expression patterns. In addition, homeobox swap experiments show that the homeoboxes of flanking HoxA genes can be functionally equivalent, whereas those of Hox genes that are further separated are functionally distinct (Zhao and Potter, 2002).

Functional redundancy of flanking Hox genes was also demonstrated by nonallelic noncomplementation. For example, Hoxb5/Hoxb6 transheterozygotes show distinct skeletal malformations, even though each single heterozygote shows no phenotype (Rancourt et al., 1995). Mice that are transheterozygous for Hoxa10 and Hoxa11 also display synergistic defects in their limbs and reproductive tracts (Branford et al., 2000). These experiments, however, just begin to define flanking Hox gene redundancies. A preferred genetic approach is to mutate multiple flanking Hox genes to more completely remove overlapping function.

Hoxa10 and Hoxa11 are crucial for the appropriate development of both the male and female reproductive tracts. Morphological differentiation of the female reproductive tract from the Mullerian ducts takes place from embryonic day (E) 15 to postnatal day 14, concurrent with expression of the AbdB HoxA genes (Taylor et al., 1997). Hoxa9 is expressed strongly in developing oviducts, with weaker expression in the uterus. Both Hoxa10 and Hoxa11 show high expression in the uterus, with slight differences in their anterior expression boundaries. Development of the male reproductive tract from the Wolffian duct closely parallels the process in the female (Joseph et al., 2009). At E17 Hoxa10 is strongly expressed in the vas deferens, with an anterior boundary reaching to the future junction with the caudal epididymis (Podlasek et al., 1999). Previous studies with Hoxa10^-/- and Hoxa11^-/- mice have shown partial anterior homeotic transformations of the uterus and vas deferens that are similar but not identical. The difference in phenotypes is most apparent in the female, with the top 25% of the Hoxa10^-/- uterus displaying an anteriorization to a more oviduct-like structure (Benson et al., 1996), whereas the entire Hoxa11^-/- uterus is thinner and smaller than in wild type (WT), resembling the oviduct (Gendron et al., 1997).

Hoxd9, Hoxd10 and Hoxd11 are also expressed in the developing male and female reproductive tracts (Dollé et al., 1991), although less well characterized. In studies of mice individually mutant for Hoxd9, Hoxd10 or Hoxd11, anatomical transformations of reproductive tracts have not been reported. Hoxd9^-/- mice are fertile and give birth to normal sized litters (Fromental-Ramain et al., 1996). Interestingly, studies of Hoxd10 and Hoxd11 mutants reveal normal fertility in females but subfertility in males that is associated with an inability to produce a vaginal plug, which is likely to result from hindlimb locomotor defects (Carpenter et al., 1997; Davis and Capecchi, 1994).

In this report, we describe the use of recombineering to study Hox gene functional redundancies in the reproductive tracts of male and female mice. With a modified recombineering strategy we were able to generate BAC targeting constructs over 100 kb in length that facilitated the simultaneous frameshift mutation of multiple flanking Hox genes, while leaving regional enhancers intact. We generated mice with flanking mutations in Hoxa9, Hoxa10 and Hoxa11 (Hoxa9,10,11) and the paralogous genes Hoxd9, Hoxd10 and Hoxd11 (Hoxd9,10,11). The Hoxa9,10,11 mutant mice provided a genetic setting that revealed previously undetected redundant fertility and reproductive tract patterning functions for Hoxd9,10,11. Dramatic anterior homeotic transformations of uterus and vas deferens were observed. Surprisingly, however, the epithelial lining of the mutant uterus underwent a conspicuous change in gene expression and structure that caused it to resemble that of the more posterior vagina. In addition, the results demonstrate that these Hox genes play an essential role in the regulation of the uterine immune system and that small non-coding RNAs are among their key targets.

BAC clones used were RP23-20F21 (HoxA cluster) and RP23-101K1 (HoxD cluster) from CHORI BACPAC Resources (http://bacpac.chori.org/home.htm). Plasmid constructs for recombineering of BACs were made by subcloning PCR-amplified DNA into a modified pl451 (Pentao et al., 2003) (supplementary material Table S1), with kan/neo flanked by ‘once only’ Lox66 and Lox71. Sequences were chosen so as to introduce a small deletion/frameshift near the beginning of the first exon. Recombineering modification of BACs was carried out as previously described (Warming et al., 2005) (precise modifications are listed in supplementary material Table S1).

The modified BACs were isolated using NucleoBond PC 500 (Clontech), linearized using PI-SCEI (New England BioLabs), and electroporated into embryonic stem cells (ESCs) (SE2 cells, strain 129 made in the S.S.P. lab.). ESC colonies surviving G418 selection were screened by quantitative PCR (qPCR) to count the number of remaining wild-type alleles, with proper recombinations leaving only one. A primer-probe set (Integrated DNA Technologies) was designed in which the binding sequence of the signaling probe was homologous to the deleted region of the Hox gene. ESCs were then further tested by PCR to check that each Hox gene was properly targeted (primer sequences are listed in supplementary material Table S1). Mice were maintained on a mixed 129/CD1 background (CD1 from Charles River).

Males and diestrus females at 3-7 months of age were used for gross analysis and histology. Tissue sections from adult WT and Hox mutant mice were processed onto the same slide for immunofluorescence (Hirota et al., 2010). Antibodies used for immunofluorescent staining were specific to SMAA (1:500; LMAB-1A4, Seven Hills Bioreagents), OVGP1 (1:250; sc-48754, Santa Cruz Biotechnology), RPS12 (1:250; 16490-1-AP, Proteintech) and KRT5 (1:500; gift from Geraldine Guasch, Cincinnati Children’s Hospital Medical Center).

Whole-mount in situ hybridization for Hoxd9,10,11 was performed as described (Benson et al., 1996). Primer sequences are listed in supplementary material Table S1.

RNA was isolated from male and diestrus female (n=3/genotype) reproductive tracts using the RNeasy kit (Qiagen). Total RNA (1 μg) was amplified using the WT Expression Kit (Ambion) and hybridized to Affymetrix GeneChip Mouse Gene 1.0 ST arrays; analysis was performed using GeneSpring GX (Agilent Technologies). Data were deposited in the GEO database with accession number GSE41993. Hierarchical clustering was performed using the Chebyshev distance metric.

RNA-Seq using Illumina TruSeq was performed on pooled total RNA (1.5 μg) from the same three samples/genotype that were analyzed by microarray. Data were analyzed in Avadis NGS as previously described (Brunskill and Potter, 2012) and deposited in the GEO database with accession number GSE42339.

Uterine horns were dissected from WT and Hoxd9,10,11^-/- diestrus mice, minced into small pieces, placed in 3 ml HBSS (Life Technologies), 0.83 Wünsch units of Liberase TM (Roche) was added and the samples incubated at 37°C for 45 minutes with intermittent pipetting. Samples were centrifuged at 1200 rpm (250 g) for 5 minutes, washed twice with 2 ml 5 mM EDTA in PBS, strained through a 70 μm cell strainer, pelleted and resuspended in FACS buffer (5% FBS, 0.1% NaN₃). Live cells present in the resulting cell suspensions were counted on a hemocytometer using a 1:2 dilution with Trypan Blue, and similar numbers of cells (<10⁶) from WT and mutant uteri were stained with fluorochrome-conjugated antibodies (1:50) for myeloid-lineage markers: CD45 (PTPRC), FITC (553079, BD Biosciences); GR1 (LY6G), PE-Cy7 (552985, BD Biosciences); F4/80 (EMR1), PE (122616, BioLegend); CD11b (ITGAM), PacBlue (101224, BioLegend). Single-stained compensation controls were prepared using a CD1 mouse uterus.

Samples were analyzed on an LSRII cytometer (BD Biosciences). Single-stained controls prepared from uterine tissue were used to calculate compensation for spectral overlap of the fluorochromes used. We collected 100,000-200,000 events from each sample that fell in the initial gating on side scatter (SSC) versus forward scatter (FSC) to include cells that are likely to be alive. Using FlowJo software (TreeStar) we gated on live, single cells using SSC versus FSC, then on the CD45⁺ population. Counts for the CD45⁺, CD11b⁺, F4/80⁺ and GR1⁺ cell populations were obtained and expressed as a percentage of the live or CD45⁺ populations. Total CD45⁺ cell counts in WT and Hoxd9,10,11^-/- uteri were determined by multiplying the percentage of CD45⁺ cells by the total number of cells harvested from the uteri as determined by Trypan exclusion using a hemocytometer. Similarly, total CD11b⁺ cell counts were determined by multiplying the percentage of CD45⁺ CD11b⁺ cells by the total CD45⁺ cell count.

To better understand the overlapping functions of flanking Hox genes it is necessary to create mice in which sets of flanking Hox genes have been mutated. One approach is to use Cre-Lox technology to generate deletions that encompass multiple Hox genes. This strategy, however, also removes interspersed shared enhancers, thereby disrupting the expression patterns of the remaining Hox genes (Zákány et al., 2004). It is therefore difficult to attribute the resulting phenotype to the loss or misexpression of particular Hox genes. A preferred approach would be to introduce frameshift mutations into multiple flanking Hox genes. We reasoned that recombineering methods could be adapted to achieve this objective.

We modified recombineering to allow the simultaneous targeting of nearby genes. The BAC targeting constructs were engineered to carry multiple frameshift mutations (Fig. 1A). A DNA segment with a kan/neo selectable marker flanked by two blocks of sequence homology to the Hox gene of interest was recombined into the first exon and, following identification of cells carrying the desired modification, the marker was removed with inducible Cre. We used the ‘once only’ LoxP sequences, Lox66 and Lox71. Each carries one mutation, and the single remaining LoxP following recombination is an inactive double mutant that does not interfere with subsequent serial modifications of additional Hox genes. The resulting BAC targeting constructs were used to produce mice with mutations in Hoxa9,10,11, the three genes of the HoxA cluster with strong uterus expression, as well as in their paralogs on the HoxD cluster, Hoxd9,10,11. The remaining kan/neo sequences were removed by breeding to mice with germ line Cre expression (EIIa-cre from Jackson Labs).

The expression of Hoxd9,10,11 in the developing reproductive tracts has not been defined in detail (Dollé et al., 1991). We found that Hoxd9,10,11 are expressed during development (Fig. 1B) and in the adult uterus (supplementary material Fig. S1) with expression patterns very similar to Hoxa9,10,11 (Ma et al., 1998; Satokata et al., 1995), placing these genes in the right place and time to impact patterning and function of the reproductive tract.

Mice homozygous mutant for Hoxa9,10,11 or Hoxd9,10,11, as well as mice deficient for nine alleles Hoxa9,10,11^-/-d9,10,11^+/- or Hoxa9,10,11^+/-d9,10,11^-/-, were viable and produced in the predicted Mendelian ratios. Hoxa11^-/-d11^-/- mice usually suffer perinatal death due to severely hypoplastic kidneys (Davis et al., 1995). As expected, therefore, Hoxa9,10,11^-/-d9,10,11^-/- mice died shortly after birth.

Both male and female Hoxa9,10,11^-/- mice were infertile, confirming and extending previous studies with Hoxa10^-/- and Hoxa11^-/- mutants (Gendron et al., 1997; Satokata et al., 1995). Male Hoxa9,10,11^-/- mice showed completely penetrant bilateral cryptorchidism, with intra-abdominal testes located just below the inferior pole of the kidney. Cryptorchidism results in defective spermatogenesis (supplementary material Fig. S2). Studies of Hoxa10- or Hoxa11-deficient males previously reported cryptorchidism linked to abnormal development of the gubernaculum, but with incomplete penetrance and variable location of undescended testes (Lewis et al., 2003; Satokata et al., 1995; Rijli et al., 1995). As expected, both male and female Hoxa9,10,11^-/-d9,10,11^+/- mice were infertile.

We also observed that Hoxd9,10,11^-/- male mice were infertile and female mice showed evidence of subfertility. Five matings of Hoxd9,10,11^-/- males with WT females failed to produce vaginal plugs over 6 months, consistent with previous reports of copulation defects in Hoxd10 mutants and subfertility of Hoxd11 mutants, even though male reproductive tracts appear normal in both (Favier et al., 1995; Carpenter et al., 1997). In crosses between Hoxd9,10,11^-/- females and WT males, the litters were somewhat smaller (8.7±1.3 pups/litter, average ± s.e.m.) than for either Hoxd9,10,11^+/- (11.0±1.0 pups/litter) or WT (11.5±1.8 pups/litter) females. In addition, litters born to Hoxd9,10,11^-/- mothers frequently died postnatally with no milk spots or were severely runted (5/12 Hoxd9,10,11^-/- litters).

We found more striking evidence for Hoxd9,10,11 function in female fertility in compound Hox mutants. About 26% (13/49) of double heterozygous (Hoxa9,10,11^+/-d9,10,11^+/-) females with a vaginal plug did not produce pups. Furthermore, females that did become pregnant produced small litters (average of 5.3 pups/litter, n=13) compared with Hoxa9,10,11^+/- (average 13 pups/litter, n=4) and Hoxd9,10,11^+/- (average 11 pups/litter, n=8) single heterozygous mice. Importantly, matings of Hoxa9,10,11^+/-d9,10,11^-/- females (n=4) with WT males failed to produce any pups despite the presence of vaginal plugs.

Hoxa10 and Hoxa11 are known to play a role in developmental patterning of the reproductive tracts (Benson et al., 1996; Gendron et al., 1997), but Hoxd9,10,11 have not previously been implicated in this process. We analyzed Hoxa9,10,11^-/-, Hoxd9,10,11^-/- and compound Hox mutant uteri, reasoning that the presence of Hoxa9,10,11 in previous studies masked the patterning function of Hoxd9,10,11.

Hoxa9,10,11^-/- mice displayed severe anterior homeotic transformations of reproductive structures; uteri were extremely thin along their entire length, resembling oviducts (Fig. 2B). The uterotubal junction, which connects the uterus and oviduct, was poorly defined, elongated and shifted posteriorly compared with WT (Fig. 2A). Histological analysis of Hoxa9,10,11^-/- uteri showed a significant reduction of both myometrial and stromal layers and an increase in luminal branching (Fig. 2H), more similar to WT oviduct (Fig. 2J) than WT uterus (Fig. 2G). Immunostaining for α-smooth muscle actin (SMAA) showed disruption of a thinned uterine muscle layer (Fig. 2N) compared with WT (Fig. 2M). In addition, the mutant uterus displayed a dramatic reduction in uterine gland tissue compared with WT. These changes reflect an apparent transformation of the mutant uterus to a more oviduct-like structure. The surprising strength of the phenotype is likely to be the result of overlapping functions of the flanking Hoxa9,10,11 genes. By contrast, the Hoxd9,10,11^-/- uterus appeared grossly normal without any obvious anteriorization (Fig. 2C). The histology and SMAA staining patterns of Hoxd9,10,11^-/- uteri (Fig. 2I,O) were indistinguishable from those of WT (Fig. 2G,M). However, deletion of one set of Hoxd9,10,11 alleles from the Hoxa9,10,11 homozygous mutant (Hoxa9,10,11^-/-d9,10,11^+/-) significantly exacerbated the anteriorization of the uterus observed in Hoxa9,10,11^-/- mice (Fig. 2E). These uteri were extremely thin and coiled for nearly half their length. Histological analysis showed that Hoxa9,10,11^-/-d9,10,11^+/- uteri closely resembled oviductal architecture, with dramatically decreased diameter, a branched luminal epithelium and severely diminished stromal and myometrial layers (Fig. 2K).

There was also evidence of functional redundancy of patterning for HoxA and HoxD paralogs in Hoxa9,10,11^+/-d9,10,11^-/- mice. Whereas Hoxa9,10,11^+/- and Hoxd9,10,11^-/- mutants showed uteri similar to WT, Hoxa9,10,11^+/-d9,10,11^-/- compound mutants exhibited a thin uterus with occasional slight coiling at the anterior end (Fig. 2F), thinner myometrial and stromal layers, as well as a more branched luminal epithelium (Fig. 2L). Even Hoxa9,10,11^+/-d9,10,11^+/- uteri often displayed an elongated uterotubal junction (Fig. 2D). These results demonstrate that Hoxd9,10,11 function in female reproductive tract patterning.

In Hoxa9,10,11^-/- males the vas deferens showed a partial anterior homeotic transformation towards an epididymis. Gross and histological analyses of the reproductive tract revealed elongation of the junction between the caudal epididymis and vas deferens, with posterior shifting and coiling of the vas deferens near the epididymal junction (Fig. 3A,B). The mutant vas deferens showed dramatically reduced muscle and stromal layers (Fig. 3G,H). Immunostaining for SMAA showed a decrease in the thickness, as well as a disorganization, of the muscle layer (Fig. 3M,N).

By contrast, the Hoxd9,10,11^-/- vas deferens appeared grossly normal, but with occasional minor coiling near the epididymal junction (Fig. 3C). Histological analysis and SMAA immunostaining showed no apparent abnormalities (Fig. 3I,O) compared with WT (Fig. 3G,M). However, removal of one set of Hoxd9,10,11 alleles from Hoxa9,10,11 homozygous mutants (Hoxa9,10,11^-/-d9,10,11^+/-) resulted in a significant enhancement of the Hoxa9,10,11^-/- vas deferens anteriorization (Fig. 3E), with extreme coiling from the epididymal junction extending posterior approximately one-third of its length. Histology and SMAA staining (Fig. 3K; supplementary material Fig. S2) revealed a vas deferens structure similar to the WT caudal epididymis (Fig. 3J; supplementary material Fig. S2), with thinner muscle and stromal layers. Removal of one set of Hoxa9,10,11 alleles from Hoxd9,10,11^-/- (Hoxa9,10,11^+/-d9,10,11^-/-) gave rise to a vas deferens that typically displayed coiling of the proximal end along with a smaller caudal epididymis (Fig. 3F). These mice rarely presented with unilateral cryptorchidism and more commonly with canalicular testes located just above the scrotum. Histology and SMAA staining of the vas deferens showed a reduced muscle layer, particularly for the anterior portion (Fig. 3L; supplementary material Fig. S2). Even the Hoxa9,10,11^+/-d9,10,11^+/- double heterozygotes usually showed some gross coiling of the anterior vas deferens (Fig. 3D). These compound mutant phenotypes further demonstrate the functional redundancy of Hoxa9,10,11 and Hoxd9,10,11 in the development of the vas deferens.

Interestingly, the size of the caudal epididymis was reduced in Hoxa9,10,11^-/-, Hoxa9,10,11^-/-d9,10,11^+/- and Hoxa9,10,11^+/-d9,10,11^-/- mice. Staining for SMAA in the Hoxa9,10,11^-/- caudal epididymis showed a thinning of the muscle layer, suggesting that this structure now resembled the more anterior caput segment of the epididymis (supplementary material Fig. S2). The anteriorization of the Hoxa9,10,11^-/-d9,10,11^+/- epididymis was even more obvious, with extreme thinning of the muscle layer in the caudal epididymis, further implicating Hoxd9,10,11 in male reproductive tract patterning.

We used microarrays to better understand the molecular basis of the observed anterior transformations in the Hoxa9,10,11^-/- uterus. In comparing mutant and WT, 128 genes were differentially expressed [fold change (FC)>2, P<0.05]. Fifty-three percent (39/74) of upregulated genes were also upregulated in WT oviducts versus uteri (FC>2, P<0.05), consistent with the observed gross anteriorization. Hoxa9,10,11 therefore normally function, in part, to repress the expression of these oviduct-related genes in the uterus. For example, Ovgp1, which encodes oviductal glycoprotein 1, normally shows much stronger expression in oviducts than uteri, but was dramatically upregulated in mutant uteri (supplementary material Table S2).

Other genes showing increased expression in mutant uteri, and associated with anteriorization, included the growth factor/cytokine genes Fgf1, Fgf16, Fgf18, Ccl8, Wnt2b, Slit2, Cxcl13 and Il1f8. Downregulated genes included the growth factor Bmp7 and the transcription factors Eya1, Lef1 and Hand2. Two microRNA (miRNA) genes, Mir10b and Mir181b-1, were also downregulated in the mutants.

The microarray data for Hoxd9,10,11^-/- uteri presented a different landscape. Consistent with the normal morphology and histology, there was very little molecular evidence of anteriorization towards oviduct. Nevertheless, there were some striking differences in the mutants. Surprisingly, Hoxd9,10,11^-/- uteri showed more genes with altered expression (189) than Hoxa9,10,11^-/- mutant uteri. The upregulated genes in Hoxd9,10,11^-/- uteri exhibited a clear inflammatory signature, with cytokine genes (Ccl6, Ccl8, Ccl9, Ccl11), antigen recognition genes (Cd209f, Clec4a1, Clec4a3) and major histocompatibility complex class II components (H2-DMa, H2-Aa), suggesting increased populations of antigen-presenting cells.

The downregulated genes were also very informative. Strikingly, nearly all encoded short non-coding RNAs (ncRNAs), including small nucleolar RNAs (snoRNAs), miRNAs and small Cajal body-specific RNAs (scaRNAs). Interestingly, one upregulated gene in the Hoxd9,10,11^-/- uterus was Rbm3, which when overexpressed can cause a global depletion of miRNAs (Dresios et al., 2005). These results show that Hoxd9,10,11 are positive regulators of small ncRNAs.

Hox paralogs are thought to have predominantly overlapping sets of downstream targets. It is therefore interesting to consider the overlap of perturbed gene expression profiles of the Hoxa9,10,11^-/- and Hoxd9,10,11^-/- uteri. Surprisingly few genes (<10%) were present in the upregulated gene lists of both mutants. This is consistent with the observed anteriorization of Hoxa9,10,11^-/- in contrast to Hoxd9,10,11^-/- uteri. There was greater overlap among the genes showing reduced expression in the mutants. Of the 58 genes with lower expression in Hoxa9,10,11 mutants, 25 also showed reduced expression in Hoxd9,10,11 mutant uteri. Furthermore, over 90% of these shared genes encoded small regulatory ncRNAs.

As more Hox genes were removed, the homeotic transformation of the uterus towards oviduct became more extreme. Comparison of significantly differentially expressed genes (P<0.05) in Hoxa9,10,11^-/-, Hoxa9,10,11^-/-d9,10,11^+/-, Hoxa9,10,11^+/-d9,10,11^-/- and WT uteri revealed a subset of 110 co-regulated genes that were changed at least 2-fold in one genotype (supplementary material Table S2). Hierarchical clustering produced a heatmap connecting molecular anteriorization and genotype (Fig. 4A). Fold changes in gene expression generally correlated with the severity of anteriorization, such that Hoxa9,10,11^-/-d9,10,11^+/- had the greatest fold change followed by Hoxa9,10,11^-/- and Hoxa9,10,11^+/-d9,10,11^-/- (supplementary material Table S2). Ovgp1 provides an example of this expression/morphology relationship, showing increases of 9.0-, 6.5- and 5.1-fold in Hoxa9,10,11^-/-d9,10,11^+/-, Hoxa9,10,11^-/- and Hoxa9,10,11^+/-d9,10,11^-/- uteri, respectively (supplementary material Table S2; Fig. 4B). Ovgp1 expression was unchanged in Hoxd9,10,11^-/- uteri, in agreement with the absence of anteriorization in this genotype.

Similar to the upregulation of anterior genes in Hoxa9,10,11^-/- uteri, many of the genes upregulated in the mutant vas deferens were those that are normally highly expressed in the more anterior epididymis. Microarray analysis of Hoxa9,10,11^-/- versus WT vas deferens found 474 differentially expressed genes, 75% (221/295) of which are normally expressed at much higher levels in the epididymis than in the vas deferens. Rps12, which is highly expressed in the epididymis and upregulated in the mutant vas deferens, provides one example (Fig. 4C). Very few genes showed changes in expression common to both uterus and vas deferens in Hoxa9,10,11^-/- mutants. These included the downregulated transcription factor Eya1 and placental Plac8 and the upregulated retinoic acid-binding protein Crabp2 and pro-inflammatory Retnla.

Compared with Hoxa9,10,11^-/-, a much smaller set of genes was differentially expressed in the Hoxd9,10,11^-/- vas deferens (91 genes with FC>2, P<0.05) (supplementary material Table S3). Much like Hoxd9,10,11^-/- uteri, there was little molecular evidence of anterior transformation, as only five of these genes are normally expressed at higher levels in the epididymis. Comparison of this dataset with paralogous Hoxa9,10,11 mutants revealed that only three of 36 upregulated genes were shared.

There was, however, a striking prevalence of small ncRNAs in the 55 downregulated genes, similar to in Hoxd9,10,11^-/- uteri. Thirty percent of the downregulated genes in the Hoxd9,10,11^-/- vas deferens were also downregulated in the Hoxd9,10,11^-/- uterus and all of them were ncRNAs. The downregulation of small ncRNAs in both male and female Hoxd9,10,11^-/- reproductive tracts underscores the potential role of these Hox genes in driving ncRNA expression.

Microarray analyses of the Hoxa9,10,11^+/-d9,10,11^-/- and Hoxa9,10,11^-/-d9,10,11^+/- vas deferens revealed expression changes that corresponded to an increased severity of anterior transformation, with 139 and 682 differentially expressed genes (P<0.05 and FC>2), respectively. Comparison of the significantly altered genes (P<0.05) in these datasets along with the Hoxa9,10,11^-/- dataset yielded 201 co-regulated genes that were changed at least 2-fold in one genotype (supplementary material Table S3). These genes were used to generate a heatmap, which clearly shows that the expression signature in the Hox mutant vas deferens was dramatically shifted to resemble that of the WT epididymis (Fig. 4D).

We then used RNA-Seq to examine the altered gene expression patterns of the Hox mutant reproductive tracts. Genes that were defined as differently expressed by microarrays were mostly confirmed as differently expressed by RNA-Seq. For example, 84% of genes identified by microarrays as differently expressed in the Hoxa9,10,11^-/- uterus were also found by RNA-Seq (supplementary material Table S4). The major reason for disagreement appeared to result from the different methods used for cDNA synthesis. The TruSeq RNA-Seq technology used an initial poly(A)⁺ selection step, whereas the Ambion microarray technology included a random primer for cDNA synthesis. Therefore, the microarrays detected gene expression differences in non-polyadenylated RNAs not seen by RNA-Seq. This difference was particularly pronounced for the Hoxd9,10,11^-/- uterus, where most of the gene expression differences detected by microarrays involved small ncRNAs that are often non-polyadenylated or only transiently polyadenylated. In this case, only 36% of the gene expression differences seen with arrays were also detected by RNA-Seq (supplementary material Table S4). The two methods were therefore complementary.

RNA-Seq analysis typically resulted in a 5- to 15-fold increase in the number of significantly differentially expressed genes identified (FC>2, P<0.05) compared with microarray analysis. Whereas the microarrays generally distinguished hundreds of gene expression differences, RNA-Seq found thousands. There appeared to be multiple reasons for this discrepancy, including lower background with RNA-Seq, the different amplification technologies employed and the limitations of microarrays in terms of gene representation.

RNA-Seq detects alternative splicing events. We observed that the Hox clusters displayed a number of interesting noncanonical splicing patterns, many of which we have previously described in the developing kidney (Brunskill and Potter, 2012). For example, in the WT uterus the last exon of Hoxd11 was frequently spliced to the last exon of Hoxd10, resulting in a frameshifted non-coding transcript (Fig. 5A). Interestingly, this splicing was absent in the Hoxa9,10,11^-/- uterus, and the same pattern held for the WT and Hoxa9,10,11^-/- vas deferens. This suggests splicing cross-regulation between Hox clusters. The Hox mutants also showed distinct RNA processing patterns for many other genes (supplementary material Fig. S3).

The predominant expression of Hoxa9,10,11 in the uterus was clearly shown in RNA-Seq data from the entire HoxA cluster (Fig. 5B). Hoxa9,10,11 were highly expressed in the WT uterus, with very little expression of other HoxA genes. The HoxD cluster also showed high levels of expression for Hoxd9,10,11, with moderate to high expression for Hoxd3,4,8 (Fig. 5C). Hoxb6,7,8 were strongly expressed in the uterus with no difference among genotypes (Fig. 5D). The HoxC cluster exhibited very low expression in the WT uterus (Fig. 5E), but paralogs of the genes targeted in this study (Hoxc9,10,11) were significantly upregulated in all Hoxa9,10,11 and Hoxd9,10,11 mutant uteri (Fig. 5F). It is also interesting to note that, although Hoxa9,10,11 and Hoxd9,10,11 clearly have overlapping functions, there was no compensatory expression of Hoxd9,10,11 in Hoxa9,10,11 mutants or vice versa.

We also observed cross-regulatory interactions with non-paralogous Hox genes in the uterus. All other HoxA genes (Hoxa1-7 and Hoxa13) and the 3′ HoxC genes (Hoxc4-8) were upregulated to varying degrees in the Hoxa9,10,11^-/-, Hoxa9,10,11^-/-d9,10,11^+/- and Hoxa9,10,11^+/-d9,10,11^-/- uteri (Fig. 5F). Hoxc4 was also upregulated in the Hoxd9,10,11^-/- uterus. In addition, Hox genes were upregulated in the mutant vas deferens, although less extensive than for the uterus.

RNA-Seq data provided a more comprehensive definition of the upregulated immune processes in Hox mutant uteri. Gene ontology analyses revealed that immune processes were significantly enriched (false discovery rate correction, P<0.05) in all Hox mutant uteri datasets, with a particularly strong signature in Hoxd9,10,11^-/- and Hoxa9,10,11^+/-d9,10,11^-/- mice (supplementary material Table S5). To define altered immune cell populations in the Hox mutant uteri we performed flow cytometric analysis. We found a dramatic 3-fold increase in the number of CD45⁺ hematopoietic cells in Hoxd9,10,11^-/- versus WT uteri (Fig. 6A). The majority of CD45⁺ cells were CD11b⁺ myeloid lineage cells in both WT and Hoxd9,10,11^-/-, but the percentage of CD11b⁺ cells in the CD45⁺ population was significantly increased in mutant versus WT uteri (88% versus 72%; Fig. 6B), and the total number of CD11b⁺ cells was nearly 4-fold higher in the mutants (Fig. 6C). Further analysis of myeloid lineage cells in mutant uteri showed an influx of either F4/80⁺ macrophages or GR1⁺ granulocytes (representative flow panels are shown in Fig. 6D). More specifically, two of seven Hoxd9,10,11^-/- uteri had an increased percentage of F4/80⁺ cells among CD45⁺ cells (72% versus a WT average of 52%) and five of seven mutant uteri had an increased percentage of GR1⁺ cells (30% versus a WT average of 7.6%). The reasons for this variability are unknown. Immunofluorescent staining of uterine sections for F4/80 also showed increased expression throughout the Hoxd9,10,11^-/- uterus (Fig. 6E). To determine whether these changes were uterus specific or reflective of alterations in the global hematopoietic system we also performed hemavet analyses on peripheral blood from the mice used for flow cytometric analyses. We found significant increases in total white blood cells and lymphocytes with trends toward increases in neutrophil and monocyte populations in Hoxd9,10,11^-/- blood (supplementary material Fig. S4).

Uterine-specific natural killer cells (uNKs) play key roles in embryo implantation and placental development as well as providing innate immune protection (Lash and Bulmer, 2011). A previous study of Hoxa10^-/- mice showed inhibited differentiation of uNKs in the decidua basalis, as detected by the absence of granzyme A expression (Rahman et al., 2006). Confirming and extending this observation, RNA-Seq analysis showed seven granzyme genes with expression downregulated 3- to 100-fold in the Hoxa9,10,11^-/- mutants, with smaller fold changes for these genes in the Hoxd9,10,11^-/- mutants (see supplementary material Table S5).

RNA-Seq analysis uncovered a squamous epithelium gene expression signature suggesting a novel posteriorization of the Hoxa9,10,11^-/- uterus. For reasons that are not clear, these changes in gene expression were not detected with microarrays. Several keratin genes typically associated with stratified, squamous epithelium were highly upregulated in Hoxa9,10,11^-/- uteri (supplementary material Table S5). This was surprising given that both the uterus and oviduct are lined with a columnar epithelium. A transition to a stratified, squamous epithelium normally occurs further posterior at the junction with the cervix. Immunofluorescent staining for KRT5 was performed on sections from the anterior, middle and posterior ends of WT and Hoxa9,10,11^-/- uteri. Strikingly, the epithelium of the Hoxa9,10,11^-/- posterior uterus displayed pseudostratification and strong KRT5 expression in the most basal epithelial cells, similar to the staining observed in WT vaginal tissue (Fig. 7A). We also found dramatically increased KRT5 expression in the middle and anterior portions of the mutant uterus, although these sections displayed normal columnar epithelium.

To better define the extent of stratification and KRT5 expression in the mutant uterus, we analyzed longitudinal sections and stitched together composite images to provide a full-length view of the Hoxa9,10,11^-/- and WT uterus (Fig. 7B). In the WT there is a very clear transition from the low KRT5-expressing columnar epithelium of the lower uterus to the brightly staining KRT5⁺ stratified epithelium of the cervix. However, the mutant displays an extended transition zone, with bright KRT5⁺ basal cells appearing in the lower uterus and KRT5⁺ columnar cells extending rostral to near the oviducts. Hoxa9,10,11^-/- uteri also showed a small but significant increase in the expression of Hoxa13 (Fig. 5F), which is typically expressed only in the cervix and vagina.

These results strongly suggest that Hoxa9,10,11 function to repress the expression of both the anterior oviductal and posterior cervix/vaginal genes in the uterus, and are therefore essential for specifying proper uterine architecture.

We have developed a novel recombineering strategy that allows the generation of BAC-based DNA targeting constructs designed to simultaneously introduce frameshift mutations into multiple flanking genes. This strategy circumvented the deletion of shared regulatory elements located between Hox genes, a complication that confounded previous Cre-LoxP-mediated Hox deletion studies. We generated mice with simultaneous mutations in Hoxa9, Hoxa10 and Hoxa11, the only three genes of the HoxA cluster with abundant expression in the uterus. We reasoned that the triple mutations might reveal reproductive tract functions masked in single mutants by flanking gene functional redundancies. We also generated triple mutant mice with frameshifts in Hoxd9, Hoxd10 and Hoxd11. Although these genes had not been implicated in reproductive tract patterning, we suspected that previously hidden roles might be uncovered by combining them with the paralogous Hoxa9,10,11 mutants.

Hoxa9,10,11 mutants showed synergistic reproductive tract anterior homeotic transformations, as defined by morphology, histology and gene expression profile. This extends previous studies that provided evidence of Hoxa10 and Hoxa11 redundancy in the female reproductive tract (Benson et al., 1996; Gendron et al., 1997; Branford et al., 2000). These results strongly suggest overlapping functions for the Hoxa9,10,11 flanking genes. In addition, we show a previously undetected functional overlap of Hoxa9,10,11 and Hoxd9,10,11 in patterning the reproductive tracts. Although Hoxd9,10,11^-/- mice showed few, if any, signs of reproductive tract anteriorization, the removal of one Hoxd9,10,11 copy from Hoxa9,10,11^-/- mice significantly exacerbated anterior homeotic transformations of the uterus and vas deferens. In addition, Hoxa9,10,11^+/-d9,10,11^-/- and Hoxa9,10,11^+/-d9,10,11^+/- mice displayed synergistic morphological changes in the reproductive tracts. The removal of at least some Hoxa9,10,11 gene activity, therefore, was necessary to unmask Hoxd9,10,11 reproductive tract patterning functions.

It is important to note that this study did not remove all Hox function during reproductive tract development. Owing to early postnatal lethality we were not able to examine adult reproductive tracts of Hoxa9,10,11^-/-d9,10,11^-/- homozygous mutant mice. In addition, Hoxb9 and Hoxc9,10,11 remained intact, and although these genes normally show very low expression levels, for example in the uterus, we did observe significant compensatory upregulation in the mutants. It remains likely, therefore, that removal of all Hox function would result in even more dramatic phenotypes.

To better understand the molecular nature of the mutant phenotypes we used microarray and RNA-Seq gene expression profiling technologies that were both cross-validating and complementary. The results identified sets of more anterior oviduct and epididymis genes that were ectopically expressed in the Hoxa9,10,11^-/- uterus and vas deferens, respectively. Although Hoxd9,10,11^-/- mutants did not show significant molecular anteriorization, when combined with paralogous HoxA mutations, as in Hoxa9,10,11^-/-d9,10,11^+/-, there were increases in anterior gene expression correlating with more complete homeotic transformations. Clearly, Hoxa9,10,11, and to a lesser extent Hoxd9,10,11, normally function to repress the expression of these oviduct and epididymis genes in more posterior structures. These effects could be either direct or indirect.

We also found that Hoxd9,10,11 play a role in female fertility. Hoxd9,10,11^-/-, Hoxa9,10,11^+/-d9,10,11^+/- and Hoxa9,1011^+/-d9,10,11^-/- mice display varying levels of reduced fertility not seen in Hoxd9^-/-, Hoxd10^-/- or Hoxd11^-/- mutants. This strongly suggests overlapping functions for these flanking and paralogous Hox genes.

Surprisingly, most of the downregulated genes in both the mutant uterus and vas deferens were non-polyadenylated small regulatory ncRNAs. The majority of these ncRNAs in the Hoxd9,10,11^-/- reproductive tracts were functionally diverse snoRNAs. A number of miRNAs were also downregulated in both the Hoxa9,10,11^-/- and Hoxd9,10,11^-/- uteri. Two of these miRNAs, Mir181 and Mir10b, were previously shown to downregulate Hox genes (Naguibneva et al., 2006; Sun et al., 2011), and Mir7 was reported to be a target of Hox regulation (Reddy et al., 2008). Our results suggest reciprocal regulation of miRNAs by Hox genes, resulting in a negative-feedback loop. For example, Hox genes normally drive increased expression of Mir181, which in turn targets Hoxa11 mRNA for degradation. During the terminal differentiation of myoblasts, Hoxa11 normally inhibits MyoD expression. Increased Mir181 expression represses Hoxa11, thereby allowing elevated MyoD expression and promoting myogenesis (Naguibneva et al., 2006). The dysregulation of Mir181b-1 and MyoD might be responsible, in part, for the disruption of muscle and aberrant muscle-related gene expression observed in the Hox mutant reproductive tracts. Importantly, myocardin, a key regulator of smooth muscle differentiation, was downregulated in the Hoxa9,10,11^-/- and both nine allele mutant (Hoxa9,10,11^-/-d9,10,11^+/- and Hoxa9,10,11^+/-d9,10,11^-/-) vas deferens.

One of the most striking results of the gene expression profiling was the dramatic upregulation of inflammatory markers in Hox mutant uteri. This signature was most pronounced in Hoxd9,10,11^-/- uteri, where nearly all upregulated genes were related to immune processes. The inflammatory response was validated by immunostaining and flow cytometry showing an influx of CD45⁺ CD11b⁺ myeloid lineage cells in mutant uteri. The increased CD11b⁺ cells were not uniform in identity among mutants. The most prevalent granulocytes in circulating blood are neutrophils, which are exclusively phagocytic cells that respond to acute infection. Macrophages are also phagocytes that typically respond to chronic infection, but have additional roles in antigen presentation and signaling. The increased numbers of these phagocytes in the Hoxd9,10,11^-/- uterus could suggest infection or inflammation in this tissue or the improper recruitment of these cells to the uterus. It is well known that infiltration of leukocytes into the uterus occurs in response to estrogen (Zheng et al., 1988), so it is possible that the increases we observed result from dysregulated levels of, or response to, estrogen in the mutants. Increased uterine leukocytes could also reflect a global increase in leukocytes. Analysis of peripheral blood did show significantly increased total white blood cells, with trends toward increased neutrophils and monocytes in Hoxd9,10,11 mutants. Interestingly, although HoxA, HoxB and HoxC genes are expressed in hematopoietic cells (Moretti et al., 1994), HoxD gene expression has not been detected in these cells. This suggests that any hematopoietic effects of Hoxd9,10,11 mutation are probably indirect.

Surprisingly, we observed an apparent posteriorization of the Hoxa9,10,11^-/- uterine epithelium. As determined by gene expression profile, histology and immunochemistry, the normal columnar epithelial layer of the uterus was partially posteriorized to a more vaginal-type stratified, squamous epithelium. This finding is anti-dogmatic, as Hox mutations almost always produce anteriorization. Interestingly, vitamin A deficiency (VAD) has been shown to produce a columnar-to-squamous transition in the mouse uterus (Ponnamperuma et al., 1999). Taken together with our results and the fact that Hox genes are targets of retinoic acid signaling (Marshall et al., 1996), this suggests that VAD could result in reduced expression of Hoxa9,10,11 in the uterus leading to squamous metaplasia. In addition, Wnt7a mutants show decreased expression of Hoxa10 and Hoxa11 and uterine posteriorization, as evidenced by stratification of the uterine epithelium and altered gene expression (Dunlap et al., 2011). These results suggest that Hoxa9,10,11 positively stamp the uterus identity, actively repelling influences from both anterior and posterior.

We thank Adam Kushner for help with maintaining the mouse colony and Shawn Smith and Hung-Chi Liang for help with the microarray and RNA-Seq experiments.