The regulatory landscape of the Dlx gene system in branchial arches: Shared characteristics among Dlx bigene clusters and evolution

The mammalian Dlx genes encode homeobox‐type transcription factors and are physically organized as convergent bigene clusters. The paired Dlx genes share tissue specificity in the expression profile. Genetic regulatory mechanisms, such as intergenic enhancer sharing between paired Dlx genes, have been proposed to explain this conservation of bigene structure. All mammalian Dlx genes have expression and function in developing craniofacial structures, especially in the first and second pharyngeal arches (branchial arches). Each Dlx cluster (Dlx1/2, Dlx3/4, and Dlx5/6) has overlapping, nested expression in the branchial arches which is called the "Dlx code" and plays a key role in organizing craniofacial structure and evolution. Here we summarize cis‐regulatory studies on branchial arch expression of the three Dlx bigene clusters and show some shared characteristics among the clusters, including cis‐regulatory motifs, TAD (Topologically Associating Domain) boundaries, CTCF loops, and distal enhancer landscapes, together with a molecular condensate model for activation of the Dlx bigene cluster.

Keywords: branchial arches; cis‐regulatory elements; distal enhancers; Dlx genes; in vivo enhancer assay

The Dlx genes are key regulator for craniofacial development and evolution. Multiple redundant enhancers located in both close and distant locations to the Dlx genes are summarized, and possible Dlx bigene cluster activating mechanism and evolution are introduced.

The Dlx genes are homeobox‐type transcriptional regulatory genes homologous to the Drosophila Distal‐less gene (Cohen & Jurgens, [4]). Mammals have six genes, Dlx1, Dlx2, Dlx3, Dlx4, Dlx5, and Dlx6, which are paired by combinations of Dlx1/Dlx2, Dlx3/Dlx4, and Dlx5/Dlx6 to form three clusters (McGuinness et al., [21]; Nakamura et al., [24]; Ozcelik, Porteus, Rubenstein, & Francke, [26]; Simeone et al., [36]; Stock et al., [38]). In the common ancestor of vertebrates before the two rounds of whole genome duplication (2R‐WGD), a tandem duplication of the Dlx gene took place, then followed by 2R‐WGD which generated the current three Dlx bigene clusters, with one cluster lost during evolution (Ellies et al., [8]). Therefore Dlx1/Dlx4/Dlx6 and Dlx2/Dlx3/Dlx5 are paralogous groups, respectively (Stock et al., [38]; Sumiyama, Irvine, & Ruddle, [39]). The expression of the paired Dlx genes is similar, suggesting that they are controlled by a common expression mechanism. Each Dlx gene is involved in early vertebrate development of the forebrain, limbs, placenta, branchial arches and many other tissues (Panganiban & Rubenstein, [27]).

The Dlx clusters maintain physical linkage to the Hox clusters. In humans, the Dlx1/2 cluster is on chromosome 2 and linked to the HoxD cluster, the Dlx3/4 cluster is on chromosome 17 and linked to the HoxB cluster, and the Dlx5/6 cluster is on chromosome 7 and linked to the HoxA cluster, respectively (McGuinness et al., [21]; Morasso, Yonescu, Griffin, & Sargent, [23]; Nakamura et al., [24]; Quinn, Johnson, Nicholl, Sutherland, & Kalionis, [31]; Scherer et al., [34]; Simeone et al., [36]; Sumiyama et al., [39]). The Dlx gene group corresponding to the HoxC cluster on chromosome 12 was lost during evolution. While the Hox genes play a role in determining positional information along the anterior‐posterior body axis (Hox‐code), the nested expression pattern of the Dlx genes (Figure 1) provide positional identity to the neural crest cells along the proximodistal axis in the branchial arches (Dlx code) (Depew, Lufkin, & Rubenstein, [5]).

Knockout of the Dlx genes is known to alter the positional information in the branchial arches and causes homeotic transformations that modify the upper and lower jaw identities (Depew et al., [5]). Whether this positional information depends on the unique protein properties of each Dlx gene (qualitative) or on the total Dlx dosage (quantitative) is yet to be determined (Depew, Simpson, Morasso, & Rubenstein, [6]). The regulatory mechanism of the nested expression pattern of Dlx genes should be important in morphological development and evolution. Here we focus on the cis‐regulatory elements in the Dlx clusters involved in branchial arch development and describe their shared characteristics among the functionally somewhat redundant Dlx clusters.

Because of the conservation of structure and sequence among various vertebrates, the intergenic region between the paired Dlx genes has been a major target in the search for distal enhancer elements (Ghanem et al., [11]; Park et al., [28]; Sumiyama et al., [40]; Zerucha et al., [47]). Several conserved elements have been tested by in vivo enhancer assays by combining with the heat‐shock protein promoter (hsp68 promoter) and the lacZ reporter gene (Kvon, [18]; Sumiyama & Ruddle, [42]).

In the Dlx1/2 intergenic region, I12a was identified as a distal enhancer element with branchial arch activity (Ghanem et al., [11]; Park et al., [28]) (Figure 2a). The I12a enhancer drives expression that partially overlaps Dlx1/2 mRNA expression in the maxillary, mandibular and hyoid processes, but does not fully recapitulate the endogenous expression pattern (Figure 3a). I12a is the only enhancer that has been identified as having branchial arch activity in the Dlx1/2 intergenic region when tested by in vivo enhancer assay, but histone modification mark H3K27ac (assumed as activation signal) and open chromatin signature (the ATAC‐seq signal, often correlates with CTCF occupancy, promoter, and distal enhancer) of developing branchial arches at E10.5 (Minoux et al., [22]) extends almost in the entire intergenic region (Figure 3a). The intergenic region also includes I12b, which is primarily a forebrain enhancer without branchial arch activity (Ghanem et al., [11]; Figures 2a and 3a).

In the Dlx5/6 intergenic region, the forebrain and branchial arch enhancer I56i was identified (Zerucha et al., [47]; Figure 2a). It drives expression in both the mandibular arch of PA1 and hyoid arch (PA2), which partially mimics the endogenous Dlx5/6 expression pattern (Ghanem, Yu, Poitras, Rubenstein, & Ekker, [13]; Zerucha et al., [47]). Almost the entire intergenic region of Dlx5/6 is marked by H3K27ac activation signal in mandibular and hyoid arches (Figure 3b), as observed in I12a. This intergenic region also includes I56ii, which is not a branchial arch enhancer but is active in the developing forebrain and in the specific GABAergic interneurons (Ghanem et al., [11]; Poitras et al., [29]; Figures 2a and 3b).

In the Dlx3/4 intergenic region, the I37‐2 element (Figure 3c) is a branchial arch enhancer showing expression at the distal tip of mandibular and hyoid arches which recapitulates the endogenous expression pattern at E10.5 (Sumiyama & Ruddle, [42]). Deletion of this I37‐2 element from the 79kb PAC transgene construct showed reduced mandibular and hyoid expression, suggesting that I37‐2 is a core cis‐element of branchial arch expression (Sumiyama & Ruddle, [42]). The intergenic region of Dlx3/4 also includes the limb enhancer I37‐1 (Sumiyama et al., [40], [39]; Figure 2a) and also exhibits H3K27ac signal in E10.5 mandibular and second pharyngeal arches in almost the entire region (Figure 3c). This indicates some contribution of the whole intergenic sequence to the branchial arch expression.

As expected from overlapping branchial arch expression, there may be some shared trans factor input among all Dlx intergenic cis‐regulatory elements. In fact, homeodomain‐factor related core binding motifs are found in all elements, some of which are conserved between I56i and I12a and between I56i and I37‐2 (Ghanem et al., [11]; Sumiyama & Ruddle, [42]). Interestingly some of the putative DLX binding motifs are shared among all three branchial arch cis‐regulatory elements (Figure 2b), indicating a possible common mechanism for establishing/maintaining branchial arch enhancer activity. It may be possible that these enhancers with the common motif were generated by 2R‐WGD, but it is not conclusive due to the limitation of aligned sequence length. So far little is known about transcription factors which directly bind and control these branchial arch enhancers. What kind of combination of transcription factors is necessary to control the branchial arch expression via these enhancers should be determined in the future, in order to clarify commonalities and differences between the intergenic enhancers.

In addition to the intergenic similarities, the Dlx clusters have some structural similarities. The 5' upstream sequences of the Dlx2/3/5 paralogous group have not been reported to have branchial arch related enhancer activities. Transgenic analysis of the Dlx3 PAC deletion clone shows that the Dlx3 upstream region has neither branchial arch nor limb bud enhancer activities (Sumiyama et al., [40]; Sumiyama & Ruddle, [42]). Interestingly there are CTCF binding site clusters immediately upstream of the Dlx2/3/5 paralogs (Figure 3a–c). These CTCF binding site clusters are bound by CTCF and cohesin and correspond to the TAD boundary. CTCF and Smc1 ChIA‐PET data in ESCs (Dowen et al., [7]; Weintraub et al., [45]) show that many long‐range interactions occur at these sites toward the 5' upstream region of the Dlx1/4/6 paralogs, and Hi‐C (Micro‐C) data (Hansen et al., [15]) also show that "flame" (Fudenberg, Abdennur, Imakaev, Goloborodko, & Mirny, [10]) or "stripe" (Vian et al., [44]) arise from this point (Figure 3a–c). These observations suggest that the Dlx clusters share a similar long‐range interaction machinery in close proximity to the promoter of the Dlx2/3/5 paralogs to recruit distal enhancers. It is probable that the Dlx promoters likely contribute little to tissue specificity, but the recruited distal enhancers are involved in tissue‐specific expression.

The intergenic branchial arch enhancers do not fully recapitulate endogenous expression and do not show very strong activity, as compared to the expression of the complete PAC transgene (Sumiyama & Ruddle, [42]). Several branchial arch enhancers have been recently identified in the Dlx clusters by extensive in vivo enhancer screening based on both sequence conservation and epigenetic signals. These branchial arch enhancers are found in the 5' upstream region of the Dlx1/4/6 paralogs at various distances from the Dlx promoters.

In the Dlx1/2 cluster, the URE2 enhancer element was identified approximately 10 kb away from the Dlx1 promoter (Ghanem et al., [12]). URE2 shows forebrain and branchial arch expression in E11.5 embryos. Expression is strong in the hyoid arch and relatively weak in the maxillary and mandibular arches (Ghanem et al., [12], [13]).

In the Dlx5/6 cluster, multiple branchial arch enhancers were identified (Birnbaum et al., [1]). eDlx#18 and eDlx#19 are located approximately 500kb away from the Dlx6 promoter. Human eDlx#18 was tested in transgenic mice and was shown to have medial branchial arch expression in E11.5 embryos (Birnbaum et al., [1]). Human eDlx#19 shows strong, endogenous Dlx5/6‐like expression in the maxillary and mandibular arches of E11.5 embryos (Birnbaum et al., [1]). In vivo enhancer assay of the mouse counterpart (meDlx#19) shows expression very similar to human expression in E11.5 embryos. This expression occurs around E10.5 and weakens after E12.5 (Sumiyama, K., unpublished data). Another branchial arch enhancer Mef2‐BA (not shown in Figure 3b) has been reported being active in E9.5 embryo (Verzi et al., [43]). Note that it locates in close proximity to the Dlx6 promoter and completely overlaps the promoter of the non‐coding RNA Dlx6os1. Smc1 ChIA‐PET and Hi‐C signals indicate that eDlx#18/eDlx#19 and the Dlx5/Dlx6 promoters are linked by CTCF‐loop and physically located close to each other in a 3D genome conformation (Figure 3b).

In the Dlx3/4 cluster, the branchial arch enhancer was identified approximately 70 kb away from the Dlx4 promoter (Ruf et al., [32]; Sumiyama, K., unpublished data). This enhancer (mU4‐Tad3) well recapitulates the Dlx3/4 expression pattern of E10.5 embryos. Note that this enhancer is not included in the 79 kb PAC transgene construct (Sumiyama et al., [40]), which also mimics the expression of the mandibular and hyoid arches, suggesting enhancer redundancy. This enhancer is also located near the CTCF peak, and the Smc1 ChIA‐PET and Hi‐C data indicate that it is physically close to the Dlx3/4 promoters in 3D genome conformation (Figure 3c).

Such an enhancer landscape suggests that similar TAD structure and CTCF loops may exist among all Dlx clusters. In the Dlx1/2 cluster, no long‐range distal enhancer has been identified so far, but it is probable that branchial arch enhancers locate somewhere far upstream within the TAD structure. The Dlx1/2 TAD appears to extend beyond 400kb in which there are a few branchial arch enhancer candidates located in the intron of the Metap1d and Slc25a12 genes (indicated by gray ellipses in Figure 3a), estimated by ATAC‐seq and H3K27ac signals. These candidates will need to be assayed for in vivo enhancer activity in the future.

As shown in the previous section, every Dlx cluster has multiple branchial arch enhancer elements which can drive similar expression, which suggests enhancer redundancy like that observed in other gene systems (Osterwalder et al., [25]). Adding to the in vivo tested enhancers, there are many other cis‐elements with ATAC‐seq/H3K27ac positive signals (Figure 3a–c). Understanding how such redundant multiple enhancers function simultaneously to activate the Dlx cluster is of great interest. Some interesting facts have been shown in previous experiments. Firstly, the ATAC‐seq and H3K27ac signals are considered good indicators of tissue‐specific enhancers, but when tested in an in vivo enhancer assay, only a fraction of the ATAC‐seq/H3K27ac positive candidates can drive lacZ expression (Sumiyama, K., unpublished data). This suggests that only some special kind of cis‐elements can initiate enhancer activity in vivo. Secondly, even in vivo test qualified enhancers can only initiate transient reporter expression but fail to maintain expression at later stages. Endogenous Dlx branchial arch expression begins around E9.5 and continues after E12.5, but the branchial arch enhancers listed here show transient expression around E9.5 to E10.5 and are not capable of maintaining later expression. PAC transgenesis (Tg) experiments of the Dlx3/4 cluster exhibited interesting features of such "transient" type of elements (Sumiyama et al., [40]; Sumiyama & Ruddle, [42]). When used in intact form, the PAC 79 kb Tg construct can drive both early phase expression around E9.5 and later stage expression after E13.5. When the I37‐2 element is deleted from the PAC, no early or later expression can be driven (Sumiyama & Ruddle, [42]). This means that PAC sequences other than the I37‐2 element include the cis‐elements needed to support and maintain proper expression, but without I37‐2 they cannot function as enhancers.

Based on these observations, we propose a putative model explaining the activation of the Dlx bigene clusters by multiple enhancers (Figure 4). To simplify explanation, we will refer to the special class of enhancers that can initiate activation as "leading enhancers" and all other cis‐regulatory elements that cannot initiate enhancer activity as "following enhancers". The intergenic leading enhancers (I12a, I56i, and I37‐2) and the distal leading enhancers (URE2, eDlx#18, eDlx#19, and mU4‐Tad3) are physically located relatively close by CTCF‐cohesin loop (Figure 4). Upon activation, transcription factors (TFs) accumulate on the leading enhancers and could form small seeds of phase‐separated molecular condensates consisting of TFs and co‐activators (Boija et al.., [2]; Hnisz, Shrinivas, Young, Chakraborty, & Sharp, [17]; Lesne, Baudement, Rebouissou, & Forne, [19]; Sabari et al., [33]; Shin et al., [35]). The seed condensates on the leading enhancers eventually fuse with each other, resulting in a larger molecular condensate and it physically excludes surrounding non‐activated chromatin (Figure 4). At this point, the TF bound following enhancers may join the condensate and enhance/maintain expression activity. Note that a single leading enhancer may be able to drive expression (Sumiyama & Ruddle, [42]), but multiple leading enhancers possibly help to form a larger and more stable condensate. The impact on Dlx expression caused by deletion of a redundant leading enhancer needs to be evaluated quantitatively in the endogenous genome context in the future. The large distance between leading enhancers helps to hold larger condensates (possibly super‐enhancers) for stable and vigorous expression. This is probably consistent with the observations that potent tissue‐specific developmental enhancers are often located far from the target promoter (Gonen et al., [14]; Lettice et al., [20]).

Mutations in the Dlx protein coding region often cause craniofacial dysplasia in humans. The Dlx3 gene has been identified as a candidate gene for tricho‐dento‐osseous syndrome (TDO) (Hart et al., [16]; Price, Bowden, Wright, Pettenati, & Hart, [30]), which cause craniofacial malformations. The Dlx4 gene is associated with orofacial clefting and abnormal jaw development (CL/P) (Wu et al., [46]). The Dlx5/6 gene is associated with split‐hand‐foot‐malformation 1 (SHFM1), which includes craniofacial abnormalities (Elliott & Evans, [9]).

Because the Dlx gene function is pleiotropic, mutation in the protein coding region can result in a severe phenotype, and therefore tissue‐specific cis‐regulatory mutations are expected to play a major role in phenotypic evolution. So far only a little is known about the evolutionary role of cis‐regulatory mutations in Dlx clusters. Regarding craniofacial variations in human, a SNP (rs10238953, A > G mutation, G is minor allele, minor allele frequency is about 0.14) in eDlx#18 cis‐regulatory element is associated with non‐pathogenic human morphological variation of chin (Claes et al., [3]). Because eDlx#18 has an enhancer activity in the medial part of the mandibular process (Birnbaum et al., [1]), it is highly probable that the eDlx#18 enhancer function is somehow modified by this mutation. Functional analysis will be needed for this SNP in the future. Because the eDlx#18 enhancer is highly conserved and other tetrapod species have A at this site, it is probable that the rs10238953 G allele has emerged and evolved in human lineage. The eDlx#18 enhancer is an interesting example of allelic variation with morphological phenotype which can be selected in the evolutionary process. The I37‐2 element shows accelerated evolution at the common ancestral branch of eutherians and metatherians, which may modify branchial arch enhancer activity (Sumiyama et al., [41]). It may be involved in theria craniofacial evolution since Dlx4 may be involved in secondly palate closure (Wu et al., [46]). As Dlx3/4 holds multiple branchial enhancers, it is not clear whether a single enhancer mutation causes the phenotypic variation. Phenotypic effect of cis‐regulation should not be concluded with a single in vivo enhancer‐reporter assay, but it also must be tested in the original 3D genomic context, because the presence of multiple enhancers may buffer a single mutation effect (Osterwalder et al., [25]).

The Dlx gene cluster is an ideal model system to study long‐range cis‐regulatory mechanisms in clustered genes, because it is relatively simple compared to a more complex cluster like the Hox genes. The Dlx clusters share some regulatory mechanisms, therefore paralogous comparative study may help further understanding. Recent advances in genome‐editing technology definitely promise to help in understanding the large scale cis‐regulation by introducing large scale deletions or replacing/inserting mutated cis‐elements into the original 3D genomic context. The leading enhancers are a particularly interesting target of genome‐editing to study their role in development and evolution. This class of enhancer tends to show deep conservation and may be related to the enhanceosome model in which the TF binding motifs keep composition and positioning (motif grammar) (Spitz & Furlong, [37]). It would be of great interest to see if a synthetic enhancer that is modified without deviating from the motif grammar can change the Dlx expression profile and associated phenotype. Both correct enhancer grammar and proper enhancer location within TAD structure should be clarified by in vivo experiments with endogenous sequence modification in the future to comprehend the cis‐regulatory element driven morphological evolution.

This work was supported by Grants‐in‐Aid for Scientific Research (18H02490), Japan. We would like to thank two anonymous reviewers and Dr. Steven Q Irvine for many helpful comments and suggestions.