Dbx2 regulation in limbs suggests interTAD sharing of enhancers

Abstract

Background: During tetrapod limb development, the HOXA13 and HOXD13 transcription factors are critical for the emergence and organization of the autopod, the most distal aspect where digits will develop. Since previous work had suggested that the Dbx2 gene is a target of these factors, we set up to analyze in detail this potential regulatory interaction.

Results: We show that HOX13 proteins bind to mammalian-specific sequences at the vicinity of the Dbx2 locus that have enhancer activity in developing digits. However, the functional inactivation of the DBX2 protein did not elicit any particular phenotype related to Hox genes inactivation in digits, suggesting either redundant or compensatory mechanisms. We report that the neighboring Nell2 and Ano6 genes are also expressed in distal limb buds and are in part controlled by the same Dbx2 enhancers despite being localized into two different topologically associating domains (TADs) flanking the Dbx2 locus.

Conclusions: We conclude that Hoxa13 and Hoxd genes cooperatively activate Dbx2 expression in developing digits through binding to mammalian specific regulatory sequences in the Dbx2 neighborhood. Furthermore, these enhancers can overcome TAD boundaries in either direction to co-regulate a set of genes located in distinct chromatin domains.

Keywords: Hox; TAD boundary; chromatin architecture; digit development; gene regulation.

FIGURE 1

Analysis of Dbx2 expression in developing digits. A‐D, WISH analysis of Dbx2, A, and Hoxa13/Hoxd13, B, Gdf5, C, and Mkx, D, in mouse embryonic forelimbs at different developmental stages. Scale bar: 250 μm

FIGURE 2

Single‐cell RNAseq analysis of Dbx2+ cell populations in developing hindlimbs. UMAP representations of the scRNAseq data from mouse E11, A, E13, B, and E15, C, mouse hindlimbs showing the expression of Dbx2, Hoxa13, and Hoxd13, as well as of different joint (Gdf5) and tendons/ligaments (Mkx, Scx) markers, ,

FIGURE 3

TAD organization around the Dbx2 locus. A, High resolution (5 kb bin size) Hi‐C map of the Dbx2 genomic region in mouse ES cells (top) and E14 embryonic cortex (bottom), and graphs showing the TAD‐separation score based on the HicFindTADs algorithm using different window size values (the curves calculated using standard parameters are displayed in gray and the average in blue). Data from Reference . B, 40 kb resolution (bin size) Hi‐C map of the Dbx2 genomic region in E12 mouse limb buds and graphs showing the TAD‐separation score (as in A). Data from Reference . On top, the gene loci are represented in blue (Dbx2) or gray boxes for other genes. C, ChIPseq profile showing the CTCF binding coverage in the Dbx2 genomic region in mouse E12 forelimbs. Arrowheads below the CTCF peaks indicate BS orientation, determined using the CTCFBS prediction tool (red: negative strand; blue: positive strand). Those BSs with a score < 5 or for which opposite orientations were predicted using different matrices are marked in gray. In the latter case, the orientation of the BS prediction with a higher score is indicated by the direction of the arrowhead

FIGURE 4

The Dbx2 regulatory landscape in mouse limb buds. A, 4C‐seq analysis of Dbx2 interactions using the Dbx2 promoter as a viewpoint in E12 mouse distal (light blue) and proximal (green) forelimbs. Profile overlap is in dark blue. Each curve marks the average profile of three independent biological replicates. Asterisks mark the region(s) displaying increased contact frequencies in the DFL vs. PFL. TADs in A and B are depicted with thick green lines. Protein coding loci are represented by blue (Dbx2) or gray boxes (all other genes) pointing toward the gene direction. B and C, 4C‐seq and ChIP‐seq analysis of H3K27ac and H3K27me3 marks in distal (light blue) and proximal (green) forelimbs, and HOXA13/HOXD13 binding profiles over the Dbx2 genomic region (below) and their respective peak calling. Profile overlap is in dark blue. Data from References , . The HOX13 bound putative regulatory elements (DLE1 to 3) are shown in red. The region framed by a dashed line in B is displayed at a higher resolution in C. The Vista enhancer mm1571 is represented by a blue rectangle. Zoom in view of the 4C‐seq and ChIP‐seq profiles from B

FIGURE 5

Putative Dbx2 enhancers are active in distal limb buds. A, Vista alignment of the DLE1‐3 regions and of the previously reported mm1571 regulatory element (depicted by purple and pink boxes, respectively). Evolutionarily conserved HOX binding sites are marked by vertical red lines at the bottom. B, X‐gal staining of embryo transgenic for the DLE1 and DLE2 regulatory sequences in E13 mouse forelimbs (left) and of the mm1571 Vista enhancer (image from Vista enhancer browser; https://enhancer.lbl.gov/ ⁵⁷). Distal limb expression of the mm1571 enhancers is marked by a red arrowhead. C and D, Quantitative PCR, C, and WISH, D, analysis of Dbx2 expression in the distal forelimb of wild‐type (wt) and DLE1^−/− littermates. Each point represents independent biological replicates; bars represent the mean ± SEM. Values are normalized to the Hmbs gene and to the wt. In B and D, Digits are numbered (I‐V) in the anterior to posterior order. Scale bar: 250 μm

FIGURE 6

HOXA13 and 5′ HOXD proteins cooperatively regulate Dbx2 expression in developing digits. A, Scheme of the different Hoxa and Hoxd paralogs expressed (in blue) in wt distal limbs and in those of mice carrying homozygote mutations disrupting or altering the expression of the mouse Hoxa13 and Hoxd paralogs. Silent genes are in gray. Red crosses represent inactivated genes. Dashed lines represent various deletions at the HoxD locus. The Del(Atf2‐Nsi) mice carry a large genomic deletion spanning the centromeric TAD flanking HoxD. They display virtually no expression of any Hoxd genes in digits. B, WISH analysis of Dbx2 expression in E12 mouse forelimbs of control and compound Hoxa13/Hoxd13 mutant mice. Scale bar: 250 μm. C and D, Quantitative PCR analysis of Dbx2 expression in the DFL of control embryos or in different Hoxa13, Hoxd13, and HoxD mutant alleles. Each point represents independent biological replicates; bars represent the mean replicate value ± SEM. P‐values are calculated based on t‐test comparison against wt values, E. Model explaining the cooperative role of Hoxa13 and Hoxd genes in Dbx2 regulation. Inactivated Hoxa/Hoxd paralogs in each mutant configuration are indicated in red. Arrow thickness represents the relative contribution of each HOX protein. Gray dashed arrows depict weak Dbx2 activation

FIGURE 7

Disruption of the DBX2 homeodomain. A, Dbx2 locus structure and predicted proteins in both control and mutated Dbx2 alleles. Blue scissors indicate the sgRNAs used for CRISPR/Cas9 genome editing. The homeodomain‐coding moiety is highlighted in pink. The first and last amino‐acid positions of the DBX2 homeodomain are indicated. Its three α‐helices (H1‐3) are depicted by light pink boxes. Gray lines represent the correspondence between the DBX2 H1‐H3 and its encoding sequence at the Dbx2 locus. B, Table showing the proportion of Dbx2 ^+/+, Dbx2 ^+/−, and Dbx2 ^−/− P14 offspring obtained from Dbx2 ^+/− X Dbx2 ^+/− crosses. C, Alcian blue and alizarin red staining of the forelimb of P7 wt or Dbx2 ^−/− littermates. Digits are numbered (I‐V) in the anterior to posterior order. Scale bar: 500 μm. D and E, Quantification of the length, D, and degree of ossification, E, of metacarpals (Mc) and phalanges (P1‐P3) of digits I to V in wt (green) or Dbx2 ^−/− (red) littermates. Bone length was calculated as the distance between the tips of the epiphysis. The degree of ossification was calculated as the ratio of the length of the alizarin red + domain and the total bone length. Each point represents a biological replicate. Blue lines depict the mean of all biological replicates

FIGURE 8

Nell2 and Ano6 expression in mouse limb buds. A, Synteny of the Dbx2 genomic region. Dbx2 is in blue and other genes in gray boxes. The red dashed rectangle depicts the deleted region reported in a human case. Gray and red dashed lines indicate syntenic relationship and orthologous gene loss, respectively. Gene insertions are in green. Scale bar: 100 kb. B and C, WISH of Dbx2, Nell2, and Ano6 in E13 mouse forelimbs, B, or in E11 to E13 whole embryos, C. D, UMAP representation of the scRNA‐seq data from E13 mouse hindlimbs for Dbx2, Nell2 and Ano6

FIGURE 9

Dbx2, Nell2, and Ano6 coregulation in mouse limb buds. A, 4C‐seq profiles showing the interactions of Dbx2, Nell2, Ano6, DLE1, and DLE2 in proximal (green) and/or distal (light blue) forelimbs (profile overlap is in dark blue). The gray contacts in the Ano6 viewpoint correspond to probably artefactual PCR product. Nell2, Ano6, DLE1, and DLE2 profiles are from a single experiment. The dashed rectangle represents the region analyzed in the right panel. C and D, Quantitative PCR, C, and WISH analysis, D, of Nell2 and Ano6 expression in E13 DFL of control and DLE1^−/− embryos. Each point represents independent biological replicates; bars represent the mean ± SEM. Values are normalized to the Hmbs gene and to the wt. Scale bar in B and F: 250 μm

FIGURE 10

TAD structure and expression of the chicken Nell2/Dbx2/Ano6 genes. A, Top, Hi‐C map of the Dbx2 genomic region in chicken HH20 wing buds and graphs showing the TAD‐separation score (bottom) using standard window size parameters (gray lines; average in blue). Protein‐coding gene loci are represented by blue (Dbx2) or gray boxes for all other genes. Data from Reference . TADs called by hicFindTADs are depicted with black boxes. Bottom, ChIPseq profile showing the CTCF binding coverage in the Dbx2 genomic region in chicken HH20 limbs. Arrowheads below the CTCF peaks indicate the orientation of binding sites as determined by using the CTCFBS prediction tool (red: negative strand; blue: positive strand). Those binding sites with a score < 5 or for which opposite orientations were predicted using different matrices are marked in gray. In the latter case, the orientation of the BS prediction with a higher score is indicated by the direction of the arrowhead. B, WISH analysis of Dbx2, Nell2, and Ano6 expression in chicken wings and in HH20/HH30 embryos. Scale bar 250 μm (limbs)/ 2 mm (embryos)

FIGURE 11

Dbx2 regulation in mouse and chicken. A, Scheme of the TAD architecture of the mouse Dbx2 genomic region (top) and its 3D organization (bottom). The Dbx2 and Nell2/Ano6 loci are depicted by blue and light purple boxes, respectively. All other genes are in gray boxes. DLE1 to 3 elements are shown in green. Green arrows indicate the regulation of DLE1 to 3 over Nell2 and Ano6 genes. B, Schemes of the TAD organization of the mouse and chicken Dbx2 genomic region and its regulation. DLE1 to 3 and Vista mm1571 (or its chicken orthologous) elements are depicted by green and brown round boxes, respectively. Because of the low resolution of the chicken Hi‐C, it was not possible to precisely resolve the location and extension of the Dbx2 interTAD domain (approximate TADs limits are depicted by red dashed lines). Green and brown arrows point to the DLE1 to 3 and neural tube enhancer regulatory activity, respectively. No expression of Dbx2, Nell2, or Ano6 was scored in the distal limb/wing of the chick embryo. In the mouse, the DLE1 to 3 sequences regulate the expression of Dbx2, Nell2, and Ano6 in the developing distal limb