Control of Hoxd gene transcription in the mammary bud by hijacking a preexisting regulatory landscape

Characterization of MB-Enriched Transcripts by RNA-Sequencing.

We complemented previous transcriptome studies of the adult MG (39–42) by a transcriptional analysis of microdissected E13.5 embryonic MBs. MB pairs 2 and 3 were thus obtained, and their expression profiles were established by RNA-sequencing (RNA-seq). To distinguish MB-specific expression from the background transcription found in the presumptive dermis and epidermis, we dissected from the same animals an equivalent area of tissue immediately adjacent to the MB (hereafter referred to as “skin”), including the nonmammary epithelium and its underlying mesenchyme (Fig. 1A). We identified 409 genes specifically up-regulated (SI Materials and Methods) and 204 down-regulated genes in the MBs compared with control skin (false discovery rate < 10%, minimum expression fold change = 1.5; Fig. S1 A–C and Datasets S1 and S2). Up-regulated genes included the mammary epithelium markers Krt8, Bmp2, Fgf17, Gata3, Pthlh, and Tbx3 (30, 43–47), as well as the mammary mesenchyme-specific genes Tbx18, Esr1, and Meox2 (30), thus validating our experimental approach (Fig. 1B and Datasets S1 and S2). Among the differentially expressed genes, 28 long noncoding RNAs (lncRNAs) were up-regulated and four were down-regulated in the E13.5 MBs, compared with embryonic skin (Fig. S1 A and B and Datasets S1 and S2).

To classify the MB-enriched genes functionally, we performed gene ontology (GO) term enrichment analysis (SI Materials and Methods) using Gorilla (48) (P < 10⁻⁴). Enriched GO terms were further clustered using Revigo (49) (similarity threshold = 0.7). Among the 20 most enriched GO term categories, we found several terms related to epithelial structures and MG development, as well as BMP signaling, previously described to be essential for MG development (e.g., refs. 43, 50) (Fig. 1C and Dataset S3).

To evaluate spatial expression patterns within the MB further, we analyzed the overlap between our list of MB up-regulated genes and an existing estimate of both ectoderm- and mesenchyme-specific genes derived from a microarray-based analysis of the posterior MB at E12.5 (30). We found that of our 409 identified MB up-regulated genes, 65 had been reported as ectoderm-specific in this previous work, whereas 19 of them were found to be enriched in mammary mesenchyme (Fig. 1D). However, the majority of the MB up-regulated genes in our dataset were expressed at similar levels in the E12.5 mammary mesenchyme and ectoderm (Fig. S1D). We think that differences in the embryonic stage and tissue [E12.5 posterior MB analyzed by Wansbury et al. (30) versus E13.5 anterior MB analyzed here], as well as the increased sensitivity of RNA-seq compared with microarrays (51), may account for the moderate overlap between the two datasets. In support of the latter, most of our MB-enriched genes were expressed at moderate to low levels in E12.5 mammary primordia (Fig. S1E), and were thus likely not detectable as being differentially expressed with microarrays. The distribution of these genes was nevertheless clearly shifted compared with the distribution of genes classified as down-regulated in our analysis, which tended to be expressed at lower levels in the E12.5 mammary primordium (Fig. S1E).

The Embryonic MG Hoxome.

Due to their potential role in MG development (discussed above), we particularly looked at the Hox gene family (Fig. S2A) and noticed that members of the HoxB and HoxD clusters were expressed at high levels in the MB compared with the nearby skin control, whereas HoxA and HoxC paralogs were either not differentially expressed or even down-regulated in these structures (Fig. 2A and Fig. S2A). Among HoxB and HoxD genes, Hoxb3, Hoxb6, Hoxb9, and Hoxd9 were significantly up-regulated in the MBs compared with the skin tissue (Fig. 2A and Fig. S2A). Hoxd3 was also up-regulated in the MBs, although to a lower level than the former genes. Although the absolute mRNA levels of Hoxd10 were also significantly up-regulated in the MB, the levels of Hoxd10 remained very low (Fig. 2A and Fig. S1F). We performed whole-mount in situ hybridization (WISH) at different stages of MB development (Fig. 2B and Fig. S2B) and detected strong Hoxd9 expression in all five MB pairs at E12.5 and E13.5, whereas neither Hoxd10 nor Hoxd11 transcripts were scored. Hoxd8 expression was also observed, but only in the first to third pairs of MBs at E12.5. Its expression was strongly down-regulated to become virtually nondetectable 1 day later (Fig. 2B).

In support of our RNA-seq data, WISH also revealed expression of Hoxb6 and Hoxb9 in all MBs, whereas Hoxa9 transcripts were only scored in the fourth and fifth MB pairs (Fig. S2B). We confirmed that Hoxc9 was not transcribed in the MBs at any of the developmental times analyzed (36), and we could not detect any Hoxd3 transcripts by WISH, despite a slight, yet significant increase, compared with skin (Fig. S2A), likely due to the different sensitivities of the two approaches. To determine the precise localization of the Hoxd8 and Hoxd9 transcripts, we cryostat sectioned WISHed embryos, and both Hoxd8 and Hoxd9 were detected in the mammary mesenchyme below the mammary epithelium (Fig. 2C). Although Hoxd9 was detected at a high level in all mammary mesenchymal cells, Hoxd8 transcripts were found in the most lateral portion of this mesenchyme and were barely detectable in its most central part at E12.5. This expression pattern likely paralleled the rapid down-regulation of Hoxd8 transcription occurring in the MB between E12.5 and E13.5. This dataset, along with previous contributions (36–38), indicates that Hox genes are differentially expressed both among the various pairs of MBs and during the development of these structures, suggesting that temporal and spatial Hox combinations may contribute to MG development and specification.

Hoxd8 and Hoxd9 Transcription Depends on Distinct Regulations.

To search and identify regulatory elements controlling Hoxd8 and Hoxd9 expression in MBs, we initially analyzed a BAC transgenic line containing the entire HoxD cluster (TgBAC^HoxD; Fig. 3A). To discriminate transgenic Hox genes from their endogenous counterparts, TgBAC^HoxD mice were crossed with the HoxD^{Del(1–13)d11Lac} line [also known as Del(1–13)d11lac] where all Hoxd genes are deleted (Fig. 3A). Backcrosses generated Del(1–13)d11lac^−/−;TgBAC^HoxD embryos where only transgenic Hoxd genes were present. Neither Hoxd8 nor Hoxd9 was expressed in the MB in these embryos (Fig. 3B), suggesting that DNA elements located outside the region covered by the BAC^HoxD were required.

We further analyzed a series of targeted deletions within the HoxD cluster spanning different intervals around the Hoxd8-to-Hoxd9 locus (Fig. 3C) and obtained contrasting results. Although Hoxd9 expression was not affected in any of these deletions (Fig. 3D, Upper), Hoxd8 transcripts were strongly down-regulated in all homozygous mutant embryos where the Hoxd4-to-Hoxd8 intergenic region was absent [HoxD^Del(1–4i) or HoxD^Del(i); Fig. 3D, Bottom]. Hoxd8 transcription was scored in other deletions maintaining this DNA region, however, such as the HoxD^Del(1–4) or HoxD^Del(10–13). Although phylogenetic foot-printing analysis of the Hoxd4-Hoxd8 intergenic region revealed the presence of sequences conserved in mammals (Fig. S3A), a transgenic construct carrying this region upstream of a LacZ reporter cassette did not display any β-gal activity in the E12.5 MBs (Fig. S3B), thus corroborating the results obtained with the BAC^HoxD transgenic fetuses. Altogether, these results indicated that although enhancers outside of the HoxD cluster drive Hoxd9 expression in the MBs, Hoxd8 transcription depends on the combined activity of regulatory elements located in its immediate 3′ vicinity and in the flanking gene deserts.

The selective expression of Hoxd8 and Hoxd9 in the MBs was not due to the incapacity of other Hoxd genes to be transcribed there, because mutant lines carrying deletions encompassing the Hoxd8 and/or Hoxd9 locus ectopically transcribed Hoxd10, Hoxd11, or Hoxd12 (Fig. S4). This ectopic expression was at least partially driven by regulatory elements located outside the cluster because it was also scored in embryos homozygous for the HoxD^Del(1–10) deletion, which removes the putative local Hoxd8 enhancer(s) [Fig. S4C; Del^(1–10)]. In all of these cases of ectopic expression, the concerned Hoxd genes were positioned closer to the T-DOM as a result of the deletion, suggesting that the T-DOM may contain embryonic MB enhancers specifically interacting with this central part of the gene cluster (Fig. S4E). Ectopic expression of Hoxd13 was nevertheless never scored in any of the mutant lines analyzed in this study, suggesting that not all promoters are responsive to the MB enhancers (Fig. S4).

An MB Enhancer Region Is Located in the T-DOM.

To verify that distal MB enhancer(s) were located in the T-DOM, we tested two large inversions that separate the HoxD cluster from either one or the other gene desert (Fig. 4A). In the HoxD^{Inv(Itga6-TgHd11LacNsi)} inversion, a 3-Mb region containing the entire C-DOM was inverted, repositioning potential enhancers far away from the gene cluster [Fig. 4A; Inv(Itga6-Nsi)]. On the other hand, the HoxD^{Inv(Itga6-Attp)} inversion disconnected the HoxD cluster from the T-DOM, still keeping its relative distance with potential C-DOM enhancers [Fig. 4A; Inv(Itga6-Attp) and SI Materials and Methods]. Although Hoxd9 expression in the MBs remained unaffected in homozygous Inv(Itga6-Nsi) embryos, it was fully abrogated in the MBs of Inv(Attp-Itga6) homozygous mice (Fig. 4B, Upper). To rule out a potential negative effect of the novel neighboring DNA sequence associated with the HoxD cluster following the latter inversion, we also analyzed a large deletion removing almost the entire T-DOM [Fig. S5A; Del(Attp-SB3)]. Embryonic lethality was prevented by balancing this allele with a deletion of the HoxD cluster [Del(1–13)d11Lac; SI Materials and Methods], and Hoxd8 and Hoxd9 expression in the mammary mesenchyme was expectedly lost in Del(Attp-TpSB3)^+/−/Del(1–13)d11Lac^+/− embryos, compared with the Del(1–13)d11Lac^+/− control littermates (Fig. S5 B and C).

To locate the MB regulatory region more precisely within the T-DOM, we used targeted deletions covering this area, including the HoxD^{Del(Attp-TpSB2)}, HoxD^{Del(TpSB2-TpSB3)}, HoxD^{Del(65-TpSB3)}, and HoxD^{Del(TpSB2-65)} alleles (Fig. 4A and Fig. S5A), which were all analyzed over the Del(1–13)d11Lac balancing allele. The expression of Hoxd8 and Hoxd9 in the MB was not significantly affected in either Del(Attp-SB2) or Del(CS65-SB3) embryos (Fig. 4B and Fig. S5B). In contrast, expression was no longer scored in Del(SB2-SB3) or Del(SB2-CS65) embryos (Fig. 4B, Bottom). It is noteworthy that this region was previously shown to contain regulatory elements controlling Hoxd gene expression both in the proximal limb buds and the intestinal cecum (17, 18). This DNA segment also contains the transcriptional start sites of Hog and Tog, two lncRNAs coregulated with selected Hoxd genes in the developing cecum (18), and whose transcription was also slightly enriched in E13.5 MBs compared with the surrounding control tissue (Fig. 5B).

This particular DNA region maps at the boundary between two sub-TADs within the T-DOM (17) (Fig. 5A), and chromosome conformation capture approaches have revealed its strong contacts with the HoxD gene cluster in either transcriptionally active or inactive cellular contexts (17, 18, 23, 52). Of note, the interactions observed between this DNA region and Hoxd9 are usually stronger than with the other Hoxd genes (Fig. 5C and Fig. S6 A and B). This privileged interaction involving Hoxd9 was not scored with any other sequence within the T-DOM, except for the other limb-specific enhancer CS65 (17), ruling out a technical bias associated with the Hoxd9 probes used for chromosome conformation capture (4C) studies. The same preferential interaction was observed when Hoxd9 and Hoxd11 4C profiles were compared from either E10.5 brain cells (where Hoxd genes are not expressed) or anterior and posterior trunk cells (expressing different combinations of Hoxd genes) (53) (Fig. S6 B and C), supporting the idea that Hoxd9 specifically contacts this region, regardless of its transcriptional activity.

A high density of DNA sequences conserved among mammals were found in the 60-kb surrounding this Hoxd9 interacting region, with six of them specifically conserved in all eutherian mammals and marsupials but not in monotremes (Fig. 5D, blue arrows). A 20-kb-large deletion was engineered, which removed four of the six eutherian-specific elements, as well as the previously reported conserved sequences CS38 to CS40 (17). Embryos homozygous for this deletion [HoxD^{Del(CS38–CS40)}, also known as Del(CS38–CS40)] were viable and displayed no major alteration in Hoxd gene expression in either the proximal limb or the cecum at E13.5. Likewise, the expression of Hoxd8 in the E12.5 MBs remained unchanged (Fig. S5D). In contrast, Hoxd9 transcripts were no longer observed in these structures (Fig. 5E), confirming that expression of Hoxd8 and Hoxd9 in the MBs depends on different regulatory mechanisms. Despite the loss of Hoxd9 expression in the MB of Del(CS38–CS40)^−/− embryos, no major alterations were observed either in the milking behavior or in the survival of the offspring of knockout females, in agreement with the apparent lack of phenotype observed in the Hoxd9^−/− mice (36).

To assess whether additional regulatory elements controlling Hoxd8/Hoxd9 in the developing MB were located outside of this deleted region, we tested two stable transgenic lines carrying random insertions of BAC clones flanking the CS38–CS40 region (Fig. 5B; BAC^T1 and BAC^T2) carrying an integration of a LacZ reporter under the control of the β-globin minimal promoter. Neither of the lines displayed any β-gal activity in the MB at any stage analyzed (Fig. 5E). Therefore, although the CS38–CS40 region is required for Hoxd9 transcription, the control of Hoxd8 transcription depends on the combined activities of local and distal regulatory elements.

A Eutherian-Specific MB Element.

This 20-kb-large DNA region necessary to drive Hoxd9 expression in the developing MB interacts strongly with the HoxD cluster in all tissues analyzed thus far (17, 18, 53). A close examination of 4C interaction profiles revealed that the strongest contacts coincide with the CS38 and CS39 conserved elements (Fig. S6D). The latter had been previously described as an enhancer driving Hoxd expression in developing limb buds (17, 52). Immediately 5′ of CS39, we identified a ca. 200-bp sequence conserved among eutherian mammals. A stable transgenic line carrying this eutherian-specific element, the CS39 sequence, and a LacZ reporter cassette revealed β-gal activity in the developing MBs at E12.5 and E13.5, but not in the E11.5 MPs (Fig. 6A), thus matching the expression of endogenous Hoxd9. Instead, the CS38 region did not display any enhancer activity when tested in a classical transgenic assay.

Phylogenetic footprinting of the sequence contained in the TgNCS39 transgenic construct (Fig. 6B) revealed four regions of high conservation (Peak1–Peak4). Whereas Peak1 is specific to therian mammals, Peak2–Peak4 matched the vertebrate conserved CS39 region (17). We generated transgenic constructs carrying different combinations of these conserved regions (Fig. 6C) to define which sequence displayed the MB enhancer capacity. All constructs carrying the Peak1 sequence activated LacZ transcription in the MBs (Fig. 6 D and E), in agreement with the conservation of Peak1 in therian mammals. In contrast, the TgPeak3-Peak4/LacZ reporter displayed β-gal staining in the limb buds but not in the E12.5 MB (Fig. 6 D and E). We next generated lentiviral reporter constructs carrying either Peak1 or Peak2 individually, and only the TgPeak1/LacZ construct was able to drive LacZ reporter expression in the developing MBs, in addition to a high background activity in the rest of the embryo (Fig. 6 D and E). Both constructs were consistently active in the limb bud. However, the TgPeak1/LacZ construct showed significantly less expression intensity and robustness than observed with the TgPeak1-Peak2/LacZ reporter. In addition, β-gal activity was detected in the vibrissae and in the incipient hair follicles, as well as in the presumptive skin dermis, although at low levels. These results suggest that, albeit Peak1 is the minimal region required for MB enhancer activity, other sequences might contribute to drive robust and specific transcription into these structures.

The Peak1 sequence [hereafter referred to as mammary bud regulatory element (MBRE)] was not detected in the monotreme platypus or in other nonmammalian vertebrates (Fig. 6B). To evaluate the regulatory potential of this region across eutherian and noneutherian mammals, we cloned a DNA segment orthologous to the murine Peak1-Peak2 sequence from both the platypus and opossum, and tested it in lentiviral transgenic assays. In agreement with its lack of conservation, the platypus transgene failed to drive LacZ activity in the developing MBs (one of nine), although it had similar expression in the developing limb (five of nine) compared with the murine sequence. The opossum sequence also displayed reproducible LacZ staining in the proximal limb (five of seven) but showed only sporadic expression in the developing MB (two of seven). Nevertheless, it could control expression in the hair follicle placodes, thus partially mimicking the reporter activity of the murine TgPeak1/LacZ construct.

Mouse Mutant Lines.

The Inv(Itga6-Attp) mutant allele was generated using the STRING method (76). The parental alleles contained a LoxP site at the Itga6 locus located 3 Mb 5′ of the cluster (77) and at an Attp located 3.6 kb 3′ of Hoxd1 (17). This allele was maintained on a B6/CBA F1 hybrid background. The deletion Del(SB2-SB3) was obtained by targeted meiotic recombination (TAMERE) (78) by recombining the LoxP sites at the SB2 and SB3 loci insertion sites (17) located 151 kb and 1,000 kb downstream of the HoxD cluster, respectively. The Del(CS38–CS40) line was produced via homologous recombination of a loxP-PGKneo-loxP cassette. The cells containing the insertion were selected by neomycin resistance. Mice were crossed with a Hprt-Cre strain to obtain Del(CS38–CS40)^+/− animals. A cassette containing the LacZ coding sequence under the control of the SV40 promoter was introduced at the HindIII site of the Hoxd1 second exon of the BAC RP23-400H17. A stable transgenic line carrying a random integration of the modified BAC was generated by pronuclear injection and maintained on a B6/CBA F1 hybrid background.

RNA-Seq Data Analysis.

We aligned raw RNA-seq reads on the mouse mm10 genome assembly using TopHat 2.0.9 (79), implemented in Galaxy (80). We computed gene expression levels based on genomic annotations from Ensembl release 82 (81). We extracted uniquely mapping reads from TopHat alignments using Galaxy tools, and we counted the number of reads attributable to each gene using HTSeq (82). We further computed reads per kilobase of exon per million mapped reads gene expression levels using Cufflinks (83).

We performed differential expression analyses using DESeq2 (84). Specifically, we contrasted a generalized linear model that explains the variation in read counts for each gene, as a function of the tissue (MB or adjacent skin), to a null model that assumes no effect of the tissue. We ran the Wald test and the likelihood ratio test (LRT) implemented in DESeq2 independently. The P values were corrected for multiple testing with the Benjamini–Hochberg approach. We posed that genes are differentially expressed if the false discovery rate (FDR) was below 0.1, if the absolute fold expression change estimated with the Wald test was higher than 1.5, and if the absolute fold change estimated with the LRT test was higher than 2.

A GO term search was performed with Gorilla (48) by comparing MB enriched gene (Wald test fold change > 1.5, FDR < 0.1) versus all expressed genes (mean base levels > 5), using a P value threshold <10⁻⁴. The identified GO terms were further analyzed using the Revigo algorithm to extract relevant GO categories by clustering closely related GO terms using a similarity threshold of 0.7 (49).

Phylogenetic Footprinting Analysis.

Sequences corresponding to the genomic intervals containing the MBRE-CS39 enhancers and the intergenic region i from different species were downloaded from the University of California Santa Cruz (UCSC) genome browser (Table S4). The sequences were aligned using the Vista Shuffle Lagan algorithm (85). Sequences were considered as evolutionary conserved when showing at least 70% base pair identity over a 100-bp window.

Chromatin Conformation Capture (4C) Sequencing Dataset Analysis.

The global quantifications of Hoxd gene interactions with the telomeric and centromeric gene deserts, as well as with the CS38–CS40 region, were calculated as in the study by Andrey et al. (17). The coordinates [in Mus musculus 9 (mm9)] used for these quantifications correspond to the following intervals: chr2:74,619,573–75,443,573, chr2:73,752,211–74,474,511, and chr2:74,958,108–75,003,828, respectively. The counts over the CS38–CS40 region were normalized to the sum of those counts scored in both the C-DOM and T-DOM.

Hox Genes in the MGs.

Our comparative transcriptome analysis identified 409 genes enriched in the E13.5 MBs, among which were several Hox genes. The mammary mesenchyme is instrumental in the development of the MGs and can induce ectopic expression of mammary ectoderm markers in the nonmammary skin epithelium of different tetrapod species (54, 55). Also, the mesenchyme underlying the ectodermal precursor of different skin appendages determines the type of structure to be formed (56, 57). Therefore, the molecular identities of these different mesenchymes must be distinctly identified. In a study by Wansbury et al. (30), the authors dissected E12.5 posterior MBs and separated the ectoderm from the underlying mammary mesenchyme, making it difficult to identify which genes are specifically expressed in this structure. Indeed, an important number of the genes identified in that study are also expressed in the adjacent epithelium and mesenchyme (30). Our data, thus, represent a substantial complement in the description of the gene network operating during MB development. We confirmed the strong expression of the Hoxa9, Hoxb9, and Hoxd9 paralogous genes (36) in the MBs and report the detectable transcription of additional Hox genes, such as Hoxd8, Hoxb3, and Hoxb6.

Hoxd8 transcription is rapidly down-regulated, however, whereas other Hox gene mRNAs remain at high steady-state levels. Also, Hoxa9 and Hoxd8 expression is restricted to a subset of MBs, with the latter expressed only in MB1 to MB3 and the former in the posteriorly located MB4 and MB5. These observations, along with the reported role of Hoxc8 in the early formation of the MPs from the milk line (36), suggest that spatial and temporal Hox combinations contribute to the patterning and development of the different pairs of MGs (Fig. 7A). Hairs and feathers are other skin-derived appendages whose early development shares similarities with MP induction and invagination (58–61). Of note, they also express different Hox gene combinations (32–34, 62). Furthermore, within the same skin structure, a considerable level of interspecies heterogeneity in the expression of these genes can be observed (31), suggesting that various Hox codes were differentially coopted during evolution for the specification of epidermal appendages. In this context of multiple Hox gene expression, it is not surprising that our deletion of the entire CS38–CS40 region, including the Hoxd9 MBRE element, did not elicit a strong abnormal phenotype. Functional redundancy between group 9 HOX proteins in these and other developing structures was indeed previously reported (9, 36, 63).

Distinct Regulations in the MB.

We show that the telomeric gene desert is required for the expression of both Hoxd8 and Hoxd9 in the embryonic MB. These two genes are nevertheless regulated by different mechanisms. Hoxd9 transcription depends on a 20-kb-large DNA segment (region CS38–CS40) containing the eutherian conserved element MBRE, whereas sequences inside the HoxD cluster have little impact upon its regulation. On the other hand, the MBRE-containing region only weakly, if at all, contributes to Hoxd8 expression in the MBs, which instead requires a 13-kb-large region located 3′ of the Hoxd8 locus within the cluster (Fig. 7B). Surprisingly, although the Del(SB2-CS65) deficiency abrogated Hoxd8 expression in the MBs, none of the BAC clones covering this region displayed enhancer activity in the MBs. One explanation is that regulatory elements driving Hoxd8 expression may be located within the SB2–SB3 region, yet outside of the transgenic BAC tested in this study. Also, the concomitant activity of various Hoxd8 enhancers located both in the 3′ vicinity of the gene and in the telomeric gene desert may be required such that neither one of these sequences alone would be sufficient. Therefore, although the precise mechanisms behind these regulations remain to be elucidated, our results reveal that Hoxd8 and Hoxd9 are differentially controlled in the developing MBs. This observation is intriguing, because long-range regulations involving the HoxD flanking gene deserts were previously reported to involve series of contiguous genes systematically, thus allowing for a coordinated function.

A Regulatory Hub.

The particular topology of chromatin domains at Hox loci could have triggered the emergence of novel enhancers, through preexisting structural and regulatory contacts (16). Noteworthy, region CS38–CS40 is located at the boundary between the two sub-TADs of the T-DOM and was reported to establish robust contacts with the HoxD cluster in all tissues analyzed thus far, regardless of the transcriptional state of the gene cluster (17, 53, 64). These constitutive contacts may be associated with the presence of several CTCF binding sites clustered in this region. Here, we show that among Hoxd genes, the strongest interaction with this region, and particularly with the CS39 enhancer, is established by Hoxd9, suggesting that this locus may help secure a robust constitutive interaction with this region, which might be used by various regulatory regions located around the CS3–CS40 region to interact with their sets of target genes in the vicinity of Hoxd9 (Fig. 7C). In this view, these constitutive contacts may serve as a regulatory hub by dragging various tissue-specific enhancers always at the same position within the HoxD cluster, centered on Hoxd9, as previously observed for the proximal limbs and the cecum (17, 18) (Fig. 7C). Also, this strong constitutive interaction with Hoxd9 could have been instrumental in the emergence of MB enhancers in mammals by providing the necessary chromatin architecture, and thus optimal conditions, for a novel regulatory sequence to evolve, a mechanism proposed for the de novo appearance of enhancers located in the centromeric gene desert (19).

Evolution of an MB Enhancer.

MGs likely originated from an apocrine gland already present before the divergence of the sauropsid and synapsid lineages, which further evolved in the synapsid lineage (the precursors of mammals). It is believed that the ancestral MG was a unit composed of a hair follicle, a sebaceous gland, and an apocrine gland, thus termed the apo-pilo-sebaceus unit (APSU) (65, 66). Monotremes have ∼200 repetitions of APSUs per MG, which secrete directly to the body surface; hence, they lack mammary papilla (or nipples). In contrast, marsupial and eutherian MGs release their secretion into mammary ducts that converge toward a larger duct opening into nipple-like structures (65). Therian MGs differ in their number and positions, as well as in the shape and function of the nipples, which have adapted according to the nursing behavior and number of the progenies (29, 67). It is therefore possible that changes in the transcription profiles of the mammary mesenchyme, which is required for the development of both structures, accompanied the evolution and diversification of MG morphology and functions across the mammalian lineages.

Our investigations on the evolutionary origin of the enhancer activity driving Hoxd9 expression in the mouse embryonic mammary mesenchyme identified a sequence (MBRE) capable of directing reporter gene expression in the mesenchyme of the developing MBs. The MBRE, however, requires the presence of sequences located nearby (within Peak2 conserved throughout vertebrates) to ensure specific activity into the MBs. We were not able to identify the monotreme MBRE, and the orthologous sequence from platypus did not display enhancer activity in the mammary mesenchyme. Of note, the marsupial orthologous sequence (from opossum) did trigger LacZ reporter expression in the mesenchyme underlying the forming hair follicle placodes, whereas expression in the developing MBs was only sporadic.

The lack of evidence about Hox gene expression in the MBs of noneutherian mammals makes it difficult to interpret such changes in the regulatory activity of the MBRE in these lineages. Also, interspecies divergence within enhancer sequences does not necessarily mean that their regulatory output would be fundamentally modified (68), and enhancers associated with genes under positive selection in mammals were shown to evolve rapidly (69). As a possible scenario however, this sequence may have acquired a nonspecific enhancer activity in APSU-associated mesenchyme in the common ancestor of metatherians and eutherians, thus preluding the actual murine MB regulatory element. Alternatively, this sequence exerted its MB activity in the therian common ancestor and further diversified toward a hair follicle-specific activity in the marsupial lineage. Our observation that a transgenic construct carrying only MBRE retains its broad APSU enhancer activity, thus resembling the opossum sequence, nevertheless suggests that this activity was an ancestral property. In this view, the MB specificity of this sequence evolved through the acquisition of responsiveness to upstream mammary mesenchyme determinants. The nature of these upstream factors is currently under study.

It was recently shown that transposable elements can acquire a regulatory potential during evolution, and thus have an impact on gene expression (e.g., refs. 70–72). On the other hand, vertebrate Hox clusters are largely devoid of transposable elements, unlike their flanking regulatory regions. Interestingly, the CS38–CS40 region is also relatively poor in transposable element insertions compared with the surrounding HoxD telomeric gene desert. This low number of transposons could be due either to the high density of regulatory elements found within this region or to the importance of this DNA segment in the establishment of a regulatory conformation at the HoxD locus.

Mouse Strains.

Mice were handled according to the Swiss law on animal protection (LPA), with the requested authorization (GE/81/14 to D.D.). Mice were raised and killed according to good laboratory practice standards. Genetically modified mice were maintained and crossed in heterozygosis. The mouse mutant lines used in this work and the primers used to genotype them are described in Table S1. The description of the Inv(Attp-Itga6) and Del(CS38–CS40) alleles can be found in SI Materials and Methods.

RNA Extraction and RNA-Seq.

Pairs of embryonic MBs 2 and 3 from E13.5 embryos were dissected in cold PBS and immediately processed for RNA extraction. To evaluate MB-specific gene expression, we dissected an equivalent portion of ectoderm and its underlying mesoderm tissue immediately adjacent to the MB. RNA extraction was performed using the RNAeasy Micro Kit (Qiagen) following the manufacturer’s instructions. The RNA-seq libraries were prepared from 100 ng of pure total RNA using the TruSeq Stranded mRNA protocol from Illumina with polyA selection. Libraries were sequenced on a HiSeq 2500 machine as single-end, 100-bp reads. The RNA-seq datasets produced in this study are publically available in the Gene Expression Omnibus (GEO) database (accession no. GSE84943).

Cloning.

The probe sequences of Hoxa9, Hoxb6, Hoxb9, and Hoxc9 were amplified using Toptaq (Qiagen) and specific primers (Table S2), and cloned into pGEMT-easy vector following the manufacturer’s instructions. For the cloning of lentiviral constructs, the target DNA was amplified by PCR using the Expand High Fidelity PCR system (Sigma) and specific primers (Table S3), and ligated in the PCR8/GW/TOPO (Thermo Fisher Scientific). Evolutionary conserved sequences were identified using the Vista alignment algorithm (see SI Materials and Methods). Coordinates used for the alignments are listed in Table S4. The opossum and platypus Peak1 and Peak2 regions were identified based on sequence conservation corresponding to the mouse orthologous CS39 Peak2. The cloned sequences correspond to the coordinates chr4:188,167,036–188,167,486 (monDom5) and Ultra514:15,095,67–15,096,077 (ornAna1) of the opossum and platypus genomes, respectively. These sequences were synthesized in vitro and inserted into the pENTR vector (GeneArt; Thermo Fisher Scientific). All of the PCR8/GW/TOPO or pENTR clones were verified by standard Sanger sequencing using the T7 primer located in the vector backbone, and the clones carrying the correct insert were recombined into the pRRLbLacGW (18) by LR reaction (Thermo Fisher Scientific). The final recombined constructs were verified by Sanger sequencing.

Lentiviral Transgenesis and β-Gal Staining.

Lentiviral transgenesis was performed as in the study by Friedli et al. (73), and E12.5 embryos were dissected and stained for β-gal activity following a standard protocol.

WISH.

Probes used in WISH were synthesized using T3, T7, or SP6 RNA polymerase (Table S2). They were subsequently purified using the RNA Easy Mini Kit (Qiagen). The WISH protocol was the protocol used by Woltering et al. (74).

Chromosome Conformation Capture (4C) Sequencing Dataset Analysis.

The 4C-sequencing data from E11.5 whole embryo (75) were downloaded from the GEO database (accession no. GSE79048), whereas the 4C-sequencing data from E10.5 brain and anterior and posterior trunk were obtained from Noordermeer et al. (53). The pipeline used to analyze the 4C-seq data is described in SI Materials and Methods.