A gene desert required for regulatory control of pleiotropic Shox2 expression and embryonic survival

Abstract

Approximately a quarter of the human genome consists of gene deserts, large regions devoid of genes often located adjacent to developmental genes and thought to contribute to their regulation. However, defining the regulatory functions embedded within these deserts is challenging due to their large size. Here, we explore the cis-regulatory architecture of a gene desert flanking the Shox2 gene, which encodes a transcription factor indispensable for proximal limb, craniofacial, and cardiac pacemaker development. We identify the gene desert as a regulatory hub containing more than 15 distinct enhancers recapitulating anatomical subdomains of Shox2 expression. Ablation of the gene desert leads to embryonic lethality due to Shox2 depletion in the cardiac sinus venosus, caused in part by the loss of a specific distal enhancer. The gene desert is also required for stylopod morphogenesis, mediated via distributed proximal limb enhancers. In summary, our study establishes a multi-layered role of the Shox2 gene desert in orchestrating pleiotropic developmental expression through modular arrangement and coordinated dynamics of tissue-specific enhancers.

Fig. 1. The Shox2 gene desert constitutes a hub for tissue-specific enhancers.

A Genomic interval containing the Shox2 TAD and previously identified Shox2-associated enhancer regions. Vista Enhancer Browser IDs (hs: human sequence, mm: mouse sequence) in bold mark enhancers with Shox2-overlapping and reproducible activities (arrowheads). The position of the human R4 enhancer driving reporter activity in the sinus venosus is indicated. B Heatmap showing H3K27 acetylation (ac) -predicted and ChromHMM-filtered putative enhancers and their temporal signatures in tissues with dominant Shox2 functions (see full Supplementary Fig. 1). Blue and red shades represent H3K27ac enrichment and mRNA expression levels, respectively. Distance to Shox2 TSS (+) is indicated in kb. Left: Shox2 expression pattern (Shox2-LacZ/+) at E11.5. C Transgenic LacZ reporter validation of predicted gene desert enhancers (DEs) in mouse embryos at E11.5. Arrowheads point to reproducible enhancer activity with (black) or without (white) Shox2 overlap. JGn, TGn, FGn: jugular, trigeminal, and facial ganglion, respectively. PA, pharyngeal arch. DRG, dorsal root ganglia. FL, Forelimb. HL, Hindlimb. TE, Telencephalon. DiE, Diencephalon. MB, Midbrain. HB, Hindbrain. MNP, medial nasal process. MXP, maxillary process. MDP, mandibular process. Reproducibility numbers are indicated on the bottom right of each representative embryo shown (reproducible tissue-specific staining vs. number of transgenic embryos with any LacZ signal). Corresponding Vista IDs of the elements tested are listed in Supplementary Table 1.

Fig. 2. 3D chromatin architecture across the Shox2 regulatory landscape in distinct tissues.

A C-HiC analysis of the genomic region containing the Shox2 TAD in wildtype mouse embryonic forelimb (FL), mandible (MD) and heart (HT) at E11.5 (see also Supplementary Fig. 2). The chr3:65977711-67631930 (mm10) interval is shown. Upper panels (for each tissue): Hi-C contact map revealing upstream (U-dom) and downstream (D-dom) domains flanking the Shox2 gene. Middle panels: Stronger (gray boxes, p < 0.01) and weaker (brown boxes, p > 0.01, <0.05) domain boundaries based on TAD separation score (Wilcoxon rank-sum test). A matrix showing normalized inter-domain insulation score (blue = weak insulation, red = strong insulation) is plotted below. Bottom panels: Virtual 4 C (v4C) using a Shox2-centered viewpoint shows Shox2 promoter interaction profiles in the different tissues. Shox2 contacting regions (q < 0.1, Supplementary Data 2) as determined by GOTHiC are shown on top. Red arrows point to chromatin domain anchors. Asterisk marks a high-density contact domain (HCD) observed only in heart tissue (chr3:66402500-66572500). Black arrow indicates reduction of internal D-dom contacts between elements inside the HCD and outside in the heart sample (see also Supplementary Fig. 2). B Top: CTCF enrichment in mESCs (gray) and newborn mouse hearts at P0 (orange). Bottom: CTCF motif orientation (red/blue) and strength (gradient). Protein coding genes (gene bodies) are indicated below. DEs, predicted gene desert enhancers validated in Fig. 1 (blue: tissue-specific activity). C C-HiC subtraction to visualize tissue-specific contacts for each tissue comparison (red/blue). Plots below display the corresponding subtracted inter-domain insulation scores. Dashed lines demarcate the HCD borders.

Fig. 3. Gene desert deletion reduces Shox2 in limb and craniofacial compartments.

A CRISPR/Cas9-mediated deletion of the intra-TAD Shox2 gene desert interval (GD^Δ) (mm10, chr3:66365062-66947168). Vista (hs) and newly identified gene desert enhancers (1-16, active in blue) are displayed along with TAD interval and CTCF peaks from mESCs. HCD, high-density contact domain (see Fig. 2). B, D ISH revealing spatial Shox2 expression in fore- and hindlimb (FL/HL), craniofacial compartments, and brain in GD^Δ/Δ embryos compared to wildtype (WT) controls at E10.5 and E11.5. Red arrowheads and red arrows point to regions with severely downregulated or reduced Shox2 expression, respectively. Red asterisk demarcates Shox2 loss in the anterior portion of the palatal shelves. White arrows indicate regions (diencephalon, DE and midbrain, MB) without overt changes in Shox2 expression. Scale bars, 500 μm (b) and 100 μm (d). C, E Quantitative mRNA analysis (qPCR) in limb and craniofacial tissues of WT and GD^Δ/Δ embryos. Box plots indicate interquartile range, median, maximum/minimum values (bars). Dots represent individual data points. ****P < 0.0001; *P < 0.05; n.s., not significant (two-tailed, unpaired t-test for qPCR). F DE9 and DE15 enhancer activities (Fig. 1C) overlap Shox2 expression in medial nasal process (MNP) and maxillary-mandibular (MXP-MDP) regions, respectively, in mouse embryos at E11.5. Asterisk marks anterior palatal shelf. “n” indicates number of embryos per genotype, or transgene analyzed, with similar results. Source data are provided in the Source Data file.

Fig. 4. Gene desert-mediated transcriptional control of cardiac Shox2 is essential for embryonic viability.

A ISH revealing Shox2 downregulation in the cardiac sinus venosus (SV) (red arrowhead) and nodose ganglion of the vagus nerve (red asterisk) in GD^Δ/Δ embryos at E10.5. White arrow indicates normal Shox2 expression in the dorsal root ganglia (DRG) of GD^Δ/Δ embryos. Scale bar, 100 μm. B Quantitative PCR (qPCR) revealing depletion of Shox2 in GD^Δ/Δ hearts compared to WT controls at E11.5. Box plots indicate interquartile range, median, maximum/minimum values (bars). Dots represent individual data points. ****P < 0.0001; n.s., not significant (two-tailed, unpaired t-test). C Co-localization of SHOX2 (green), HCN4 (red) and NKX2-5 (blue) in hearts of GD^Δ/Δ and WT control embryos at E11.5. SHOX2 is lost in the HCN4-marked SAN pacemaker myocardium in absence of the gene desert (dashed outline). Nuclei are shown in gray. Scale bars, 50 μm. D SAN enhancer candidate regions in the gene desert interval (VS-250) essential for Shox2 expression in the SV. Top: Virtual 4 C (v4C) Shox2 promoter interaction signatures in embryonic hearts (HT) and limbs (FL) (gray) at E11.5 overlapped with HT (red) and FL (green) -specific subtraction profiles. Below: ATAC-seq tracks from embryonic hearts at E11.5 and SAN cells from sorted Hcn4-GFP mouse hearts at P0 (Supplementary Data 3)^,. Desert enhancers (DEs) (black) and putative SAN enhancer elements with distance to the Shox2 TSS in kb (+) are indicated. Cons, vertebrate conservation track by PhyloP. E Transgenic LacZ reporter validation in mouse embryos at E11.5. Left: the +325 element drives transgenic LacZ reporter expression exclusively in the SV. Right: the +325A subregion drives Shox2-overlapping SV activity, similar to +325B (Supplementary Fig. 5). The interval shared between +325A/B subregions contains a conserved core element (marked gray) that interacts with TBX5 in embryonic hearts at E12.5. “n” denotes fraction of biological replicates with reproducible results. Single numbers represent the total of transgenic embryos obtained, including those without staining. RA, right atrium. RV, right ventricle. OFT, outflow tract. Corresponding Vista Enhancer IDs (mm) are listed in Supplementary Table 5. Source data are provided in the Source Data file.

Fig. 5. Enhancer-mediated transcriptional robustness safeguards Shox2 in the heart.

A H3K27 acetylation ChIP-seq (H3K27ac) and RNA-seq profiles from human fetal heart compartments at post conception week 17 (pcw17) across the human orthologous sequence of the +325-mouse sinus venosus (SV) enhancer and the SHOX2 interval. The left ventricle (LV) dataset has been previously published. +268, distance to SHOX2 TSS. Cons, mammalian conservation by PhyloP. B Top: Generation of a + 325 SV enhancer deletion (4.4 kb) allele in mice (SV-Enh^Δ). Below: Shox2 mRNA distribution (ISH) in SV-Enh^Δ/Δ compared to WT mouse embryos at E10.5. Arrowhead points to downregulated Shox2 in the SV. Asterisk and arrow mark Shox2 expression in the nodose ganglion of the vagus nerve and dorsal root ganglia (DRG), respectively. C qPCR analysis of Shox2 and Rsrc1 mRNA levels in SV-Enh^Δ/+ and SV-Enh^Δ/Δ embryonic hearts at E10.5 compared to WT controls. Box plot indicates interquartile range, median, maximum/minimum values (bars) and individual biological replicates (n). P-values are shown, with ****P < 0.0001 (two-tailed, unpaired t-test). Three outliers, two datapoints of Shox2 Δ/+ replicates and one for Rsrc1 (Δ/Δ), are outside of the scale shown. N.s., not significant. “n” indicates number of biological replicates analyzed, with similar results. LA, left atrium. RA, right atrium. RV, right ventricle. Source data are provided in the Source Data file.

Fig. 6. The gene desert encodes a series of distributed proximal limb enhancers (PLEs) with subregional specificities.

A Re-processed ChIP-seq datasets from mouse embryonic limbs at E10.5 (CTCF, ATAC-seq, H3K27ac) and E12 (H3K27ac, H3K27me3) showing epigenomic profiles at the Shox2 locus^,,. Bars above each track represent peak calls across replicates. DFL, distal forelimb. PFL, proximal forelimb. On top: TAD extension in mESCs (black bars) with desert enhancers identified in Fig. 1 (DEs 1-16; blue indicates validated activity). The extension of the gene desert deletion (GD^Δ) and SV control region (VS-250) is indicated. B 4C-seq interaction profiles from two independent biological replicates (R1, R2) of proximal forelimbs at E12.5 (red outline). Black arrowhead indicates the 4C-seq viewpoint at the Shox2 promoter. Gray arrowheads point to CTCF-boundaries of the Shox2-TAD. Green lines (in A) indicate Shox2-interacting elements with putative proximal limb activities. Gray lines mark 4C-seq peaks overlapping previously validated DEs without such activities (Fig. 1C). C Identification of proximal limb enhancers (PLEs) through transgenic LacZ reporter assays in mouse embryos at E12.5. Embryos shown are representatives from stable transgenic LacZ reporter lines (Supplementary Fig. 6A). Reproducibility numbers from original transgenic founders are listed for each element (bottom right).

Fig. 7. Limb-specific loss of gene desert function leads to Shox2 downregulation and defective stylopod morphogenesis.

A Schematics illustrating gene desert inactivation (GD^Δ) in the presence of reduced limb Shox2 dosage based on Prx1-Cre-mediated Shox2 deletion (Shox2^Δc). B Shox2 transcript levels determined by qPCR in fore- (FL) and hindlimbs (HL) of wildtype (WT), Shox2^Δc/+ and GD^Δ/Shox2^Δc embryos at E11.5. One outlier (FL WT datapoint) is outside of the scale shown. C Micro-CT scans of fore (FL)- and hindlimb (HL) skeletons of GD^Δ/Shox2^Δc and Shox2^Δc/+ control mice at postnatal day 42 (P42). Red arrowheads point to severely reduced stylopods in GD^Δ/Shox2^Δc individuals compared to controls (black arrowheads). “n”, number of biological replicates with reproducible results. Scale bar, 5 mm. All images at same scale. D Micro-CT stylopod quantification at P42 reveals significant reduction of stylopod (humerus/femur) length in GD^Δ/Shox2^Δc mice compared to WT (P = 9.6 × 10^-14/P = 9.6 × 10^-14), Shox2^Δc/+ (P = 1.27 × 10^-13/P = 9.6 × 10^-14) and GD^Δ/+ (P = 1.27*10^-13/P = 9.6 × 10^-14) controls. Box plot indicates interquartile range, median, maximum/minimum values (bars) with dots representing individual biological replicates (n). ****P < 0.0001; ***P < 0.001; **P < 0.01; n.s., non-significant (two-tailed, unpaired t-test for qPCR; ANOVA for micro-CT). Source data are provided in the Source Data file.

Fig. 8. Graphical summary.

A Identification of the Shox2-flanking gene desert as a reservoir for distributed transcriptional enhancers with Shox2-overlapping activities in limb, craniofacial, cardiac, and cranial nerve/neuronal cell populations. Light green indicates limb enhancers with highly subregional or limb type restricted activities. B The Shox2 gene desert encodes distributed tissue-specific enhancers that are englobed in a dynamic chromatin domain (D-dom) with tissue-invariant loop anchors and a cardiac-specific high-density contact domain (HCD) that may influence the activity of contained enhancers. Additional gene desert enhancers are likely to participate in the regulation of cardiac Shox2 in SAN progenitors. Dashed line demarcates the VS-250 interval essential for cardiac Shox2 expression. C Cumulative functions of gene desert enhancers orchestrate pleiotropic Shox2 expression essential for proximal limb morphogenesis, craniofacial patterning, and cardiac pacemaker development.