Conservation and Variations of Bimodal HoxD Gene Regulation During Tetrapod Limb Development

In all tetrapods examined thus far, the development and patterning of limbs require the activation of gene members of the HoxD cluster. In mammals, they are controlled by a complex bimodal regulation, which controls first the proximal patterning, then the distal structure, allowing at the same time the formation of the wrist and ankle articulations. We analyzed the implementation of this regulatory mechanism in chicken, i.e. in an animal where large morphological differences exist between fore-and hindlimbs. We report that while this bimodal regulation is globally conserved between mammals and avian, some important modifications evolved at least between these two model systems, in particular regarding the activity of specific enhancers, the width of the TAD boundary separating the two regulations and the comparison between the forelimb versus hindlimb regulatory controls. Some aspects of these regulations seem to be more conserved between chick and bats than with the mouse situation, which may relate to the extent to which forelimbs and hindlimbs of these various animals differ in their functions. AUTHOR SUMMARY The morphologies of limbs largely vary either amongst tetrapod species, or even between the fore-and hindlimbs of the same animal species. In order to try and evaluate whether variations in the complex regulation of Hoxd genes during limb development may contribute to these differences, we compared their transcriptional controls during both fore-and hindlimb buds development in either the mouse, or the chicken embryos. We combined transcriptome analyses with 3D genome conformation, histone modification profiles and mouse genetics and found that the regulatory mechanism underlying Hoxd gene expression was highly conserved in all contexts, though with some clear differences. For instance, we observed a variation in the TAD boundary interval between the mouse and the chick, as well as differences in the activity of a conserved enhancer element (CS93) situated within the T-DOM regulatory landscape. In contrast to the mouse, the chicken enhancer indeed displayed a stronger activity in fore-than in hindlimb buds, coinciding with the observed striking differences in the mRNA levels. Altogether, differences in both the timing and duration of TAD activities and in the width of their boundary may parallel the important decrease in Hoxd gene transcription in chick hindlimb versus forelimb buds. These differences may also account for the slightly distinct regulatory strategies implemented by mammals and birds at this locus, potentially leading to substantial morphological variations.

In all tetrapods examined thus far, the development and patterning of limbs 23 require the activation of gene members of the HoxD cluster. In mammals, they are 24 controlled by a complex bimodal regulation, which controls first the proximal patterning, 25 then the distal structure, allowing at the same time the formation of the wrist and ankle 26 articulations. We analyzed the implementation of this regulatory mechanism in chicken, 27 i.e. in an animal where large morphological differences exist between fore-and hindlimbs. 28 We report that while this bimodal regulation is globally conserved between mammals and 29 avian, some important modifications evolved at least between these two model systems, 30 in particular regarding the activity of specific enhancers, the width of the TAD boundary 31 separating the two regulations and the comparison between the forelimb versus hindlimb 32 regulatory controls. Some aspects of these regulations seem to be more conserved 33 between chick and bats than with the mouse situation, which may relate to the extent to 34 which forelimbs and hindlimbs of these various animals differ in their functions. 35

AUTHOR SUMMARY 37
The morphologies of limbs largely vary either amongst tetrapod species, or even 38 between the fore-and hindlimbs of the same animal species. In order to try and evaluate 39 whether variations in the complex regulation of Hoxd genes during limb development 40 may contribute to these differences, we compared their transcriptional controls during 41 both fore-and hindlimb buds development in either the mouse, or the chicken embryos. 42 We combined transcriptome analyses with 3D genome conformation, histone 43

INTRODUCTION 57
Tetrapod limbs are organized into three parts bearing skeletal elements; the 58 stylopodium (humerus/femur), the zeugopodium (radius/fibula, ulna/tibia) and the 59 autopodium, the latter including the acropod (phalanges, metacarpals/metatarsals) and the 60 mesopodium (carpals and tarsals) [1]. Limbs can display large variations in their 61 morphologies, either between tetrapod species or within the same species, between fore-62 and hindlimbs, as a result of their adaptation to different functions and ecological niches. 63 For example, frogs display particular shapes of carpal and tarsal elements, with an 64 elongated proximal tarsal whenever detectable[2], whereas geckos forelimb skeletal 65 elements somehow resemble those of their hindlimbs [3]. Another example of this 66 morphological flexibility are the forelimbs of bats, which do have digits early on similar 67 to those of other mammals but that subsequently elongate to make flight possible [4]. 68 In this context, birds are a fascinating animal group as they evolved forelimbs 69 (wings) and hindlimbs (legs) specialized for flying or for terrestrial locomotion, 70 respectively [5]. Recent studies using comparative genomics approaches either amongst 71 birds or between bats and mice, have revealed that some DNA sequences potentially 72 involved in limb development and which are highly conserved can display differential 73 enhancer activities as compared to their mouse orthologous sequences [6,7]. Furthermore, 74 the analysis of several domestic pigeons displaying variations in foot feathering within 75 the same species, suggested that changes in cis-regulatory elements in the genes encoding 76 forelimb-or hindlimb-specific transcription factors may contribute to a partial 77 transformation from hindlimb to forelimb identity [8]. Taken together, these observations 78 suggest that both the gain of species-specific enhancers, the different activities of the 79 same regulatory sequences and alterations in DNA sequences amongst various species 80 and/or within the same species, contributed to generate these large morphological 81

differences. 82
In addition to their essential role during axial patterning and organogenesis in 83 vertebrates[9,10], Hox genes are required for proper growth and skeletal patterning of 84 tetrapod limbs. In particular, genes belonging to the HoxA and HoxD clusters for both 85 fore-and hindlimbs, as well as some genes of the HoxC cluster during hindlimb 86 development [11,12]. Regarding both Hoxd and Hoxa genes, recent development in 87 In this context, it was recently reported that a bat regulatory sequence located within 123 T-DOM and controlling Hoxd genes displays differential enhancer activity in the limbs, 124 when compared to its mouse orthologous sequence [6]. These findings suggest that some 125 variations in limb morphology may be associated with the mechanism of bimodal gene 126 regulation described at the HoxD locus. However, this mechanism was reported to be 127 implemented during the development of forelimb buds only and hence it remained unclear 128 as to how much the important differences between fore-and hindlimbs either amongst 129 various species or within the same species, may be related to variations of this mechanism. 130 To tackle this issue, we used a comparative regulatory approach involving chick and 131 mouse embryonic fore-and hindlimbs, mostly for two reasons. First chicken embryos, 132 unlike mice, display striking differences between the morphologies of their adult wings 133 (forelimbs) and legs (hindlimbs) (Fig. 1A, B, left). Secondly, it was reported that Hoxd 134 gene expression domains during chick wing-and leg buds development showed important 135 deviations when compared to their mouse counterparts [23,31]. These features suggested 136 that the bimodal regulatory system at work at the mouse HoxD locus may be operating 137 slightly differently during the development of the avian appendicular skeletons. 138 Here we combine analyses of transcriptome, 3D genome conformation, histone 139 modification and mouse genetics to show that this bimodal regulatory mechanism is 140 highly conserved in birds. However, in chicken leg buds, the duration of the T-DOM 141 regulation is much reduced, as suggested by dynamic changes in histone modifications, 142 along with a slight difference in chromatin interaction profiles at specific enhancer 143 regions, which accounts for the concurrent reduction in Hoxd genes expression in the 144 zeugopod. Finally, by using mutant mouse embryos lacking a large part of T-DOM, we 145 uncovered differences in the effect of T-DOM deletion for Hoxd genes between fore-and 146 hindlimbs. Therefore, while these regulatory mechanisms are virtually the same amongst 147 tetrapod species and even within the same species between the fore-and hindlimbs, slight 148 differences are scored, which may partly contribute to the observed morphological 149 differences. 150

Transcription of Hoxd genes in mouse and chick limb buds 153
We first compared the expression patterns of Hoxd genes in mouse fore-and 154 hindlimbs at E12.5 (Fig. 1A) with those observed in chick at either HH28 (equivalent to 155 E12.25 to E12.5) (Fig. 1B) or HH30 (equivalent to E13 to E13.5) (Fig. 1D). In mouse 156 fore-and hindlimbs, expression of Hoxd13 and Hoxd12 was strong in the prospective 157 acropod region (hereafter termed "distal"), whereas the domains of Hoxd11 and Hoxd10 158 transcripts were split into the distal and zeugopod (hereafter termed "proximal") regions, 159 except the future mesopod, which was labelled by a Col2a1-expressing domain (Fig. 1A,  160 arrowheads). These expression patterns were similar in both fore-and hindlimbs, except 161 for a clearly weaker expression level in the hindlimb proximal domain. 162 When compared with Hoxd gene expression in the mouse counterparts, at least two 163 dramatic differences were confirmed in chick limbs. First, unlike the Hoxd12 expression 164 pattern observed in murine limbs, the chick Hoxd12 gene was strongly expressed in chick 165 proximal forelimb (wing) (Fig. 1B). Secondly, the expression of all Hoxd genes was 166 significantly reduced in the chick proximal hindlimb (leg) by stage HH28, when 167 compared to chick proximal wing and mouse proximal limbs[23,31]. As a result, the 168 expression domains of the chick Hoxd12 in wing buds was much like that of Hoxd11 or 169 Hoxd10, in contrast with the mouse where Hoxd12 was never strongly expressed in the 170 proximal domain. However, the transition between the two Hoxd-expressing domains 171 also labelled the future wing mesopod (Fig.1B, arrowheads). Of note, expression of all 172 Hoxd genes was weak in proximal leg buds, again in contrast to what was observed during 173 mouse limb bud development (Fig. 1). 174 To characterize these differences in more detail, we performed RNA-seq analyses 175 by using HH30 limb buds in order to more easily micro-dissect the various domains and 176 thus exclude a potential contamination of the future mesopod from the distal domain. 177 RNA profiles confirmed the differences detected by WISH. First, Hoxd11 to Hoxd8 were 178 expressed at lower levels in the mouse proximal hindlimb when compared to forelimb, a 179 situation re-enforced in chick proximal leg where only a weak Hoxd12 signal was 180 detected (Fig. 1C, D, upper tracks). In contrast, more reads were scored for Hoxa10 to 181 Hoxa11 in both mouse and chick proximal hindlimb when compared to forelimb (S1 Fig). In the distal domains, transcription patterns and profiles from mouse and chick were 183 similar between fore-and hindlimbs, at both HoxA and HoxD clusters (Fig. 1C, D, lower 184 tracks, S1 Fig). However, the chick profile revealed a consistently higher transcription of 185 Hoxd12, which was as strong as Hoxd13, whereas the amount of transcripts in the mouse 186 counterpart was much lower than that for Hoxd13 (Fig. 1C, D, lower tracks). Taken 187 together, these initial results indicated that both the expression status and the transcript 188 domains of Hoxd genes displayed significant differences, either amongst species or within 189 the developing fore-and hindlimb buds, in particular in chicken. 190

Bimodal regulation in both fore-and hindlimb buds 192
To determine to what extent these differences could result from variations in the 193 implementation of the bimodal regulatory mechanism, we performed comparative 4C-194 seq analyses. We used a variety of viewpoints located at comparable positions to reveal 195 potential interactions in both mouse and chicken limb buds. To do this, we cross-196 annotated Hoxd genes regulatory sequences identified in the mouse genome onto the 197 chick genome by using the LiftOver tool in UCSC. These annotations were then used for 198 all following experiments. In both fore-and hindlimbs, interactions were scored between 199 Hoxd genes and the regulatory sequences island III and Prox, as hallmarks of C-DOM 200 transcriptional activity, or with the CS39 sequence as a proxy for T-DOM activity in either 201 the distal or the proximal region, respectively [17,18]. As seen in mouse forelimbs, 202 Hoxd11 mainly contacted T-DOM sequences, especially CS39, in mouse proximal 203 hindlimb cells, i.e., in cells where T-DOM was fully active and where C-DOM was silent 204 ( Fig. 2A, top). In contrast, in mouse distal hindlimb cells, it preferentially interacted with 205 C-DOM sequences such as island III and Prox ( Fig. 2A, bottom). Quantification of 206 contacts indicated 75% of telomeric contacts in proximal forelimb cells and 50% in distal 207 forelimb cells, showing that Hoxd11 had reallocated 25% of its global interactions toward 208 the C-DOM TAD in distal cells. Likewise, mouse hindlimb cells showed the same 209 interaction profiles, with 70% of telomeric contacts in proximal hindlimb cells and 45% 210 in distal hindlimb cells ( Fig. 2A). This comparison indicated that the bimodal regulation 211 is virtually identical between fore-and hindlimbs in mouse. 212 We then examined the related interaction patterns in chick wing and leg cells by 213 using a region between Hoxd11 and Hoxd10 as a viewpoint (termed 'Hoxd10-11'), i.e. a 214 sequence located as close as possible to the bait used in the mouse experiments. In chick 215 proximal wing cells, Hoxd10-11 mostly interacted with the CS39 and CS93 regions 216 located in T-DOM, as well as with a region near the Hnrnpa3 gene, where the distal TAD 217 border is observed in the murine locus (Fig. 2B, black arrowhead). These predominant 218 contacts with T-DOM were reduced in chick distal wing cells, as in the mouse, and 25% 219 common feature of bats and chicken, suggesting that the chick CS93 sequence may have 242 a limb enhancer activity similar to that reported for the bat BAR116. We examined the 243 enhancer activity of chick CS93 using a transgenic mouse lacZ reporter system and 244 compared it to the activity of full-length mouse CS93 sequence by using lentivector-245 mediated transgenesis [32,33]. We initially cloned a 2kb large sequence containing chick 246 CS93 and those surrounding sequences (Fig. 3B), which showed higher interactions with 247 Hoxd10 to Hoxd11 in the 4C profiles obtained from proximal wing cells (Fig. 2B, track 248 1). We noted that the surrounding sequences are not particularly conserved among species, 249 whereas the CS93 region of the chick genome is more conserved with the bat than with 250 the mouse counterpart (Fig. 3B). By using the BLAT search tool in UCSC, we also found 251 that most of conserved regions from the bat BAR116 and the chick CS93 sequences can 252 be aligned onto the mouse CS9 region, respectively (Fig. 3A). 253 We assessed their enhancer activities and, unlike for the mouse BAR116, the full-254 length mouse CS93 triggered lacZ transcription in E10.5 limb buds with an expression 255 localized to the prospective stylopod and zeugopod at E12.5 (Fig. 3C). The staining was 256 weaker in hindlimb than in forelimb buds, possibly due to the delay in limb 257 development [34,35]. Accordingly, the 380bp localized in 5' of the mouse CS93 seem to 258 be essential for expression. On the other hand, we found that the 2kb large sequence 259 containing the chick CS93 displayed limb enhancer activity in mouse limb buds at E12.5 260 (Fig. 3D). The transgene expression driven by chick CS93 came under two different 261 patterns. The first one displayed lacZ staining throughout the forelimbs (n=2/5), which 262 was similar to the staining observed when the bat BAR116 sequence was assessed in 263 mouse forelimb bud (Figs 3 and 4 in[6]). The second pattern concerned mostly the 264 proximal forelimb buds (n=3/5), as seen when the mouse CS93 was used (Fig. 3C). In 265 both cases, a weaker expression was observed in hindlimb bud, as in the case of bat 266 BAR116. These results suggest that the downregulation of Hoxd genes in chick leg bud 267 involves a weaker activity of-and less interactions with the CS93 sequence. Taken 268 together with previous research [6], these results indicate that, unlike for the bat BAR116 269 and the chick CS93, the mouse CS9 is devoid of limb enhancer activity. Accordingly, the 270 380bp localized in 5' of the mouse CS93 seem to be essential for expression. 271

Implementation of the regulatory switch between TADs in mouse and chicken 273
We observed that the differences in Hoxd12 expression were more substantial 274 between chick fore-and hindlimb buds than in the murine counterparts (Fig. 1C, D).  The murine Hoxd9 to Hoxd11 genes, but not Hoxd12, are located within the TAD 286 boundary and interact both with T-DOM and with C-DOM. In contrast, in chicken limb 287 buds, Hoxd12 was able to switch contacts from T-DOM to C-DOM, suggesting that the 288 TAD boundary in chick could be located between Hoxd12 and Hoxd13 (see also S2B and 289 These results showed that the bimodal regulatory mechanism and the sequential 297 transition from the proximal to the distal global controls are implemented during chick 298 limb development similarly to what was described in mice. Therefore, the differences in 299 gene expression observed both between mice and chicken and between chick wing-and 300 leg buds cannot be solely explained by visible variations in the respective interaction 301 profiles. Instead, they ought to involve the distinct use of enhancers (or group thereof-) 302 within an otherwise globally conserved chromatin architecture.  At stage HH20 (approximatively E10 in mouse), the H3K27ac enrichment in T-322 DOM was still substantial in wing buds. In marked contrast, however, this level appeared 323 dramatically reduced in leg buds ( except that after its initial onset, T-DOM activity was terminated much more rapidly than 338 in the wing bud, followed by a decrease in accumulation of H3K27ac at the target HoxD 339 cluster itself. wing-and leg bud cells at HH20 (Fig. 5A). This extended the conclusion reached after 393 the 4C analyses, i.e. that the TAD boundary region in chick was displaced towards the 5' 394 part of the gene cluster when compared to mouse limb bud cells [19]. 395 In mouse limb cells, this TAD boundary falls within a region where multiple CTCF 396 sites are occupied [43][44][45]. CTCF is an architectural protein that both helps defining 397 constitutive domains of interaction and triggers enhancer-promoter contacts [46]. We thus 398 examined the presence of bound CTCF at the chick HoxD locus and surrounding TADs 399 and found that the profiles were comparable between wing-and leg buds at HH20 (S5 400 Noteworthy, the orientation of the CTCF motives were conserved between mouse and 404 chick. However, we found less bound CTCF in the chicken HoxD cluster than in the 405 mouse counterpart, which could affect the strength and/or stability of the TAD boundary 406 in chick (S5C Fig, bottom). 407 Using CHi-C at 5kb resolution, the distribution of contacts was relatively similar 408 between wing-and leg bud cell populations, despite the reduced level of H3K27ac in the 409 T-DOM and near the TAD border in leg bud cells described above. The contacts between 410 the HoxD cluster (galGal5, chr7:16,315,000-16,395,000) and the T-DOM region 'a' 411 (galGal5, chr7: 16,045,000-16,255,000) were indeed comparable between wing-and leg 412 bud cells (Wilcoxon sum test, W=222,690, p-value=0.663) (Fig. 5B). However, when the 413 comparison was further limited to a 30kb large region including CS93 (galGal5, 414 chr7:16,090,000-16,120,000) rather than to the entire region 'a', the contact intensity in 415 leg bud cells was clearly below that scored in wing bud cells (W=5,495, p-value=0.0213). 416 These reduced contacts between Hoxd genes and the surroundings of region CS93 417 confirmed the observation made by looking at the 4C profiles obtained using the later 418 stage (HH30) (Fig. 2). 419 The fact that bound CTCF was not detected around the CS93 region suggests that suggesting that these stable contacts are associated with CTCF, as described in mouse 430 developing limb buds [48]. 431 432

Regulation of T-DOM by HOXA13 433
We asked what may cause both the robust reduction in H3K27ac marks in chicken 434 T-DOM and the decrease in contacts between Hoxd genes and the CS93 region in leg bud 435 cells at HH20. As it was reported that HOX13 proteins bind T-DOM-located sequences 436 concomitantly to its inactivation and that the absence of these proteins prevented T-DOM 437 to terminate operating[20,49], we re-assessed the expression dynamics of Hoxa13. We RT-qPCR experiments, using chick and mouse entire fore-and hindlimb buds from HH20 443 to HH22 and E10.5 to E10.75, respectively (Fig. 6). While these developmental stages 444 are not strictly equivalent between chick and mouse[50], they were selected because the 445 size difference between chick wing-and leg buds is not yet too large between HH20 and 446 HH22[51]. Also, Hoxa13 starts to be expressed in mouse forelimb buds at around 447 E10.5 [52]. While the onset of Hoxa13 expression was detected in chick wing bud at 448 HH22, Hoxa13 transcripts were already well present in chick leg bud at HH20-21 ( Fig.  449 6A, left). Also, the expression level of this gene in leg buds was markedly stronger than 450 in wing buds (Fig. 6A, right). Hoxa11 expression was also higher in chick leg buds than 451 in forelimb buds, as was also observed in the RNA-seq dataset (S4C and S6A Fig),  452 suggesting that the entire chicken HoxA cluster was activated in leg buds well before it 453 was switched on in wing buds. This was nevertheless not a general phenomenon for Hox 454 genes, for the expression of Hoxd13 was comparable between wing-and leg buds (Fig.  455 6A, bottom and S6B, S6C Fig, left). 456 In the mouse, the development of the forelimb bud precedes that of hindlimb buds 457 by about half a day. In contrast, the initiation of both wing-and leg bud in chicken is 458 almost concomitant and the growth of the leg bud exceeds that of the wing bud [34,51]. 459 However, even when considering these developmental differences, the dramatic 460 variations we scored between both the timing of Hoxa13 activation and its transcript 461 levels between the chick wing-and leg buds were different from the situation observed 462 in murine fore-and hindlimb bud (Fig. 6B, top). Also, expression of other Hoxd genes 463 increased in both mouse fore-and hindlimb buds, as development proceeded (Fig. 6B   RT-qPCR (Fig. 7A right, S7A Fig. left). In such mutant proximal forelimb buds, Hoxd11 493 to Hoxd8 transcripts were depleted for more than 90% when compared to control 494 proximal forelimbs (Fig. 7A right, S7A Fig. left). Surprisingly however, Hoxd11 495 transcripts were not as dramatically affected in proximal mutant hindlimbs, and 496 transcripts produced by the Hoxd10 to Hoxd8 genes were decreased in amount by 50 to 497 60% when compared to control animals ( Fig. 7A right, S7A Fig. left). The reduced level 498 of Hoxd gene expression elicited by the mouse T-DOM deletion in the forelimb bud thus 499 mimics the situation observed in chick proximal leg bud (S7A Fig. right) and these 500 significant differences also exist in the way T-DOM operates in murine forelimb versus 501 hindlimb buds. 502 Interestingly, the sustained expression of Hoxd genes in T-DOM deletion mutant 503 proximal hindlimb buds completely disappeared when a larger deletion was engineered 504 between the Mtx and Titin (Ttn) genes (Fig. 7A), indicating that the genomic regions 505 between SB3 and Ttn, i.e. telomeric to the T-DOM TAD, contributes to the difference in 506 gene expression observed between the mouse fore-and hindlimb buds when T-DOM is 507 removed. 508 To identify potential forelimb to hindlimb differences in chromatin re-509 organization after the deletion of T-DOM, we generated 4C profiles from the mutant allele 510 using the Hoxd11 promoter as a viewpoint (Fig. 7B). In control proximal fore-and 511 hindlimb cells, Hoxd11 mostly contacted the intact T-DOM, with a particularly strong 512 interaction with and around the CS39 region (Fig. 7B, tracks 1 and 3). In proximal cells 513 deleted for T-DOM, interactions within the HoxD cluster were increased and ectopic 514 contacts were established (or strongly re-enforced) with the newly fused neighboring 515 telomeric TAD (Fig. 7B, tracks 2 and 4). We used this 4C-seq dataset to determine three 516 candidate regions of potential enhancer activity, referred to as hindlimb enhancer (HE) 1 517 to 3 (Fig. 7B, C, red arrows). We cross-checked this selection either with DNaseI 518 hypersensitive sites (HS) data from E11.5 hindlimb buds (GSM1014179)[53], with 519 potential enhancer regions as defined by the Limb-Enhancer Genie tool[54] and with 520 H3K4me1 ChIP-seq datasets obtained from control and mutant hindlimb proximal 521 domains ( Fig. 7C and S7B Fig). According to these criteria, HE1 turned out to be the 522 most promising region and we thus assessed its enhancer potential in developing 523 transgenic limb buds. 524 By using a transgenic enhancer reporter system, the HE1 region reproducibly drove 525 lacZ expression in proximal fore-and hindlimb buds, lateral plate and somatic mesoderm 526 at E12.5 (Fig. 7D). This result indicated that the HE1 enhancer activity is not specific for 527 the proximal hindlimb, even though it was potentially active in a hindlimb-specific 528 manner after deletion of T-DOM. To confirm this observation, we performed a 4C-seq 529 analysis by using the HE1 sequence as a bait (S7C and S7D Fig) transcripts is decreased when compared to forelimb buds (Fig. 1C). Also, the deletion of 539 TADs flanking the HoxD cluster, we wondered whether the structures of these TADs were 561 somehow modified in birds, or at least whether they would show some variations either 562 between the two species, or between the bird wing and leg buds. 563 A global analysis of 4C and capture Hi-C datasets did not reveal any salient 564 differences between mammals and birds regarding the way they implement this complex 565 bimodal limb regulation at the HoxD locus. The TADs appeared well conserved between 566 the two species, as well as the presence in chick of most -if not all-regulatory elements 567 that had been described in the mouse counterparts, on both sides on the gene 568 cluster[17,18], even though the chick TADs were reduced in size. We thus conclude that 569 the bimodal regulatory strategy described in mammals (see [59]) is implemented in a 570 similar manner during bird development. This re-enforces the idea that the function of 571 Hox genes at these early steps of limb development is mostly to set up and organize the 572 basic plan of the future appendages, rather than elaborating or fine-tuning on a pre-573 patterned structure. 574 575

Interspecies comparison of the TAD boundary at HoxD 576
The bimodal regulatory mechanism and the transition between the two global 577 controls at HoxD are thus well conserved amongst tetrapods. However, the distinct 578 expression of Hoxd12 in proximal limbs between mouse and chick suggests that the width 579 of the TAD boundary at the HoxD locus may slightly vary between the two species. By 580 using transcriptome, 4C and Hi-C datasets, we previously observed slightly different 581 positions of this boundary in mouse distal versus proximal limb cells, due to the fact that 582 Hoxd10 and Hoxd11 respond to both T-DOM and C-DOM regulations. We thus proposed 583 that the boundary was located between Hoxd11 and Hoxd12 in proximal cells, while 584 between Hoxd9 and Hoxd10 in distal cells [17,19] (Fig. 8A). 585 In contrast, the chick Hoxd12 is strongly expressed in proximal wing buds, 586 suggesting that the TAD boundary, in this context, expands towards the 5' part of the gene 587 cluster, close to the Hoxd13 locus (Fig. 8B) observed. However, when the full length CS93 sequence was injected, a robust enhancer 617 activity was scored in a proximal limb region (Fig. 3). This discrepancy between two 618 experiments involving almost the same sequences may be caused by the positions of the 619 regions used for the mouse transgenic enhancer assays, the mouse sequence extending a 620 bit more towards one of its extremities. Either the enhancer activity was provided by this 621 sub-fragment, or this fragment may be required for the expression of a more widespread 622 activity of the full DNA sequence. It remains that the BAR116 sequence may not be a bat 623 specific enhancer. were obtained with various distal to proximal distributions of the lacZ staining, a clear 631 imbalance was scored between forelimb and hindlimb cells, with a much stronger 632 expression in the former than in the latter. Therefore, the chick enhancer sequence 633 behaved more like the bat sequence than like their murine counterparts, an effect that was 634 re-enforced by the sequence alignments, which revealed more similarities between chick 635 and bats than between the two mammalian species. This similarity correlates with Hoxd 636 gene expression and may relate to the salient distinctions in morphologies between fore-637

Premature termination of T-DOM regulation in chick hindlimb buds 640
The  (Fig. 4). Instead, 647 H3K27me3 marks are still present over C-DOM at this early stage, unlike in the early 648 mouse limb buds ( Fig. 4 and S3 Fig)[17], suggesting that Hoxd13 early activation in chick 649 may be driven by the T-DOM regulation until the C-DOM regulation is implemented and 650 takes it over. This idea is well supported by our CHi-C analysis showing that the TAD 651 boundary is expanded towards the Hoxd13 locus in the early chick limb buds. 652 In addition, a major difference was observed in the activation of Hoxa13 between 653 chick and mouse hindlimb buds, with an earlier and stronger activation of Hoxa13 in 654 chick leg buds at HH20, when compared to both mouse hindlimb buds and chick wing 655 buds. This indicates that the T-DOM activity may be readily terminated by the early 656 presence of HOXA13 protein and that the C-DOM regulation is ready to be active earlier 657 in chick leg buds than in wing buds. The potential causes for both this early activation of 658 Our genetic approach however makes it difficult to assess whether this sequence is 677 used for Hoxd regulation under normal circumstances, or whether the interactions 678 observed are mostly due to its new proximity to the target genes induced by the deletion 679 of T-DOM. In the former case, this may indicate that as in chick and bats, the global C-680 DOM regulation may be more active in forelimb than in hindlimb buds and hence the 681 HE1 enhancer may not be necessary. In the case of the mouse, this deficit of regulation 682 during proximal hindlimb development could have been compensated for by evolving 683 additional enhancers outside the TAD (Fig. 8C)  Total RNA was extracted using the RNeasy Micro Kit (QIAGEN), following the 733 manufacturer instructions. 1 µg total RNA was used for cDNA synthesis with SuperScript 734 VILO (Invitrogen). RT-qPCR was performed on a CFX96 real-time system (BIORAD) 735 using the GoTaq qPCR Master Mix (Promega). Each RT-qPCR was carried out with at 736 least two biological replicates and each experimental information is described in S1 Table. 737 Primer sequences for qPCR are listed in S2 Table. 738

4C-seq and Data analysis 740
Each mouse and chick limb tissue were fixed separately with 2% formaldehyde, lysed 741 and stored at -80°C. Samples were digested with NlaIII and DpnII as primary and 742 secondary restriction enzymes, respectively, and ligation steps were performed using high 743  Table). PCR products were multiplexed and 745 sequenced with 100bp single-end, followed by post-processing (de-multiplexing, 746 mapping, and 4C analysis) using the HTS station (http://htsstation.epfl.ch) [70]. Fragment 747 scores were normalized to the mean score of fragments falling into a region defined as 748 the bait coordinated ± 1 Mb, and the data were smoothed by using a running mean with 749 a window size of 11 fragments. The information regarding fragments excluded during the 750 procedure is provided in S3 Table. A   Conservation plot of mouse CS93 and bat BAR116 using the 2kb region of chick CS93 921 as a reference. The peaks represent a conservation higher than 50%. Pink regions are 922 conserved non-coding sequences. The phylogenetic tree shows the highest conservation 923 of the chick CS93 with the bat BAR116 sequences. (C) Enhancer activities of mouse 924 CS93 and chick CS93 in mouse fore-and hindlimb buds from E10.5 to E12.5. The lacZ 925 expression pattern (left) showed that mouse CS93 has an enhancer activity in the proximal 926 region of developing limb buds at E10.5 to E12.5. In contrast to the mouse, the chick 927 CS93 (right) showed differential enhancer activity between fore-and hindlimb buds at 928 E12.5, as was also reported for the bat BAR116 sequence. The numbers of lacZ positive 929 embryos over total transgene integrated are indicated.  Mouse HE1 is active in the proximal fore-and hindlimb buds and in the trunk at E12.5. 1001 The number indicate stained embryos over total number of integrations. PFL, proximal 1002 forelimb; PHL, proximal hindlimb; HL, hindlimb. 1003  Table  1037 DNA sequences of primers used for both RT-qPCR analyses and genotyping. 1038 S3 Table  1039 DNA sequences of primers used for 4C analyses. Custom barcodes (4bp shown by 1040 NNNN) were introduced in between the Illumina adapter sequence and the specific 1041 viewpoint sequences to do multiplexing using different samples with the same viewpoint. 1042

S4 Table 1043
Public datasets used in this research. 1044 LEGENDS TO SUPPLEMENTAL FIGURES 21 S1 Fig (related to Fig 1). Hoxa genes expression in mouse and chick limb buds. proximal or distal wing-and leg buds at HH20 and HH28. Stronger enrichments were 50 observed in both whole leg buds at HH20 and proximal leg buds at HH28, when compared 51 to the corresponding samples from wing buds. The Y axis represents the strand-specific 52 RNA-seq read counts, normalized by the total number of million mapped reads. H3K27ac around the 5' Hoxa genes were observed in leg buds at both HH19 and HH20, 61 whereas less marks were scored at HH20, in the region covering the e10 to e16 enhancers, 62 when compared to wing-and leg buds at HH19. At HH28, profiles established from 63 proximal or distal region were comparable between wing-and leg buds. (B) 3' Hoxa 64 S2