An evolutionarily conserved cis-regulatory element of Nkx3.2 contributes to early jaw joint morphology in zebrafish

The acquisition of movable jaws was a major event during vertebrate evolution. The role of NK3 homeobox 2 (Nkx3.2) transcription factor in patterning the primary jaw joint of gnathostomes (jawed vertebrates) is well known, however knowledge about its regulatory mechanism is lacking. In this study, we report a proximal enhancer element of Nkx3.2 that is deeply conserved in gnathostomes but undetectable in the jawless hagfish. This enhancer is active in the developing jaw joint region of the zebrafish Danio rerio, and was thus designated as jaw joint regulatory sequence 1 (JRS1). We further show that JRS1 enhancer sequences from a range of gnathostome species, including a chondrichthyan and mammals, have the same activity in the jaw joint as the native zebrafish enhancer, indicating a high degree of functional conservation despite the divergence of cartilaginous and bony fish lineages or the transition of the primary jaw joint into the middle ear of mammals. Finally, we show that deletion of JRS1 from the zebrafish genome using CRISPR/Cas9 leads to a transient jaw joint deformation and partial fusion. Emergence of this Nkx3.2 enhancer in early gnathostomes may have contributed to the origin and shaping of the articulating surfaces of vertebrate jaws.

mechanism is lacking. In this study, we report a proximal enhancer element of Nkx3.2 that is 23 deeply conserved in gnathostomes but undetectable in the jawless hagfish. This enhancer is 24 active in the developing jaw joint region of the zebrafish Danio rerio, and was thus 25 designated as jaw joint regulatory sequence 1 (JRS1). We further show that JRS1 enhancer 26 sequences from a range of gnathostome species, including a chondrichthyan and mammals, 27 have the same activity in the jaw joint as the native zebrafish enhancer, indicating a high 28 degree of functional conservation despite the divergence of cartilaginous and bony fish 29 lineages or the transition of the primary jaw joint into the middle ear of mammals. Finally, 30 we show that deletion of JRS1 from the zebrafish genome using CRISPR/Cas9 leads to a 31 transient jaw joint deformation and partial fusion. Emergence of this Nkx3.2 enhancer in 32 Introduction 36 The establishment of jaw joints was one of the major events that enabled the evolutionary 37 transition from jawless to jawed vertebrates. The earliest articulated jaws are found in fossil 38 placoderms from the Silurian period 423 MYA (million years ago) (Zhu et al., 2013). The 39 primary jaw joint of non-mammalian gnathostomes, including actinopterygians, amphibians, 40 reptiles and birds is located within the first pharyngeal arch and articulates Meckel's 41 cartilage and the palatoquadrate. These cartilages are derived from cranial neural crest cells 42 that migrate into the first pharyngeal arch and later ossify into anguloarticular and quadrate 43 bones, respectively (Schilling and Kimmel, 1994;Tucker et al., 2004). 44

45
The transition from jawless to jawed vertebrates and the underlying gene regulatory 46 network changes are not yet fully understood. However, the transcription factor Nkx3.2 47 (Bapx1), which acts as a chondrocyte maturation inhibitor (Provot et al., 2006), is thought to 48 have played a major role in the evolution of the primary jaw joint (Cerny et al., 2010). Nkx3.2 49 displays focal expression in the first pharyngeal arch, between Meckel's 50 cartilage/anguloarticular and the palatoquadrate/quadrate in non-mammalian vertebrates, 51 as shown in zebrafish (Miller et al., 2003), Xenopus (Square et al., 2015), and python and 52 chicken (Anthwal et al., 2013). In the lamprey, a jawless vertebrate, Nkx3.2 is expressed in 53 the ectomesenchyme surrounding the pharyngeal arches (Miyashita, 2018). cartilage (Piotrowski et al., 2003). In addition to Tbx1, the same binding region matches 151 many other Tbx factors as they all possess highly similar binding motifs, and some of these 152 Tbx factors are also known to function in craniofacial development (Papaioannou, 2014  As JRS1 appears functionally conserved and active specifically in the jaw joint, we next 305 investigated the function more closely by using CRISPR/Cas9 genome editing to generate a 306 zebrafish line with the entire JRS1 enhancer sequence deleted from the genome. Two 307 sgRNAs were targeted to sequences either side of JRS1 resulting in a 445bp deletion 308 spanning the conserved core of the JRS1 and ~100 bp flanking sequences ( Figure 5A, 309 Supplementary Figure 3). The resulting enhancer deletion allele is termed nkx3.2 Δ JRS1 . We 310 generated homozygous mutant zebrafish (nkx3.2 Δ JRS1/ Δ JRS1 ) and used Alcian blue and alizarin 311 red staining to characterise the craniofacial morphology in comparison to heterozygotes and 312 wild-type siblings at 5, 9, 14, and 30 dpf. 313 314 At all examined ages, the jaw joint morphology of wild-type and heterozygous mutants 315 appeared indistinguishable. At 5 dpf, homozygous nkx3.2 Δ JRS1/ Δ JRS1 mutants tended to display 316 a reduced RAP on Meckel's cartilage and the cartilage of the palatoquadrate often failed to 317 form a pronounced convex joint process ( Figure 5D). The result of this was a lack of a clearly 318 defined concave-convex articulation in the jaw joint. Additionally, while no mutants 319 displayed a complete fusion between the Meckel's and palatoquadrate cartilages, there was 320 often a small number of chondrocytes (white arrowhead, Figure 5D) connecting the two elements. Mouth gape angle did not appear to be affected, suggesting that the joint was at 322 least partially flexible at this stage.  wild-type and heterozygous siblings ( Figure 5G). Meckel's cartilage and the palatoquadrate 337 appeared to be fully separated, and the palatoquadrate joint process was more convex, 338 resulting in a more concave-convex shape to the joint articulation. The RAP began to ossify 339 as normal in all genotypes. By 30 dpf, the jaw joint of homozygous mutants was 340 indistinguishable from wild-type and heterozygote siblings, with ossification of Meckel's 341 cartilage, the RAP, and the palatoquadrate progressing normally ( Figure 5J). 342

343
The relative expression levels of nkx3.2 were quantified in JRS1 mutants using qPCR to 344 determine the contribution of JRS1 to gene expression at 6 dpf. Larvae were individually 345 genotyped using tail tissue and then pooled into 3 biological replicates of 7 larvae per 346 genotype. Following RNA extraction and cDNA synthesis, qPCR was performed and nkx3.2 347 expression of heterozygous and homozygous JRS1 mutants was compared relative to wild-348 type. No significant differences between the genotypes were detected at this time point 349 ( Figure 5K). 350 Meckel's cartilage shape in 9 dpf wild-type, nkx3. To quantify the subtle differences in posterior Meckel's cartilage shape in JRS1 mutant 362 larvae, the heads 9 dpf nkx3.2 +/+ (n=12), nkx3.2 +/ Δ JRS1 (n=10), and nkx3.2 Δ JRS1/ Δ JRS1 (n=8), were 363 imaged using OPT, and the whole heads (Supplementary Figure 4) and both jaw joints were 364 rendered as maximum projections for a total of 60 jaw joints. The shape of Meckel's 365 cartilage at the joint interface was analysed using 2D geometric morphometrics ( Figure 6A), 366 confirming the previously described observations that there was no significant difference 367 between wild-type and heterozygous mutants (P=0.20) while homozygous mutants differed 368 significantly compared to them both (P=0.0015). Analysis of morphospace also confirmed 369 the tendency for homozygous mutants to display a reduced RAP resulting in a less concave 370 surface interfacing with the palatoquadrate ( Figure 6A). This phenotype is not fully 371 penetrant, as there is some overlap between the range of wild-type and mutant shapes, 372 consistent with the variable severity of mutant joint phenotypes seen at 5 dpf. Averaged 3D 373 renderings of the jaw joint at 9 dpf also displayed the partial fusion between Meckel's 374 cartilage and the palatoquadrate previously described at 5 dpf ( Figure 6D, G). Finally, the JRS1 enhancer was knocked out in zebrafish to reveal a transient jaw joint 388 dysmorphology and partial fusion. 389

390
We combined gene synteny analysis with non-coding sequence conservation analysis with 391 mVISTA and detected a conserved putative nkx3.2 enhancer sequence JRS1 in a number of 392 gnathostome species (Figures 1, 2). Next, we applied Tol2-mediated transgenesis to test 393 As such, the identification of JRS1 as a jaw-joint specific enhancer of Nkx3.2 that is 449 conserved in gnathostomes but absent from the jawless hagfish suggests that JRS1 evolved 450 in the gnathostome stem group after the split with Cyclostomata, and may have played a key 451 role in the early evolution of jaws. It is also noteworthy that JRS1 appears to be functionally 452 conserved in mammals, as the primary jaw joint region, along with its Nkx3.2 expression 453 domain has moved to the middle ear, forming the joint between the malleus and incus 454 The jaw joint phenotype in homozygous JRS1 mutants appears to be rescued by 486 approximately 14 dpf, suggesting that either nkx3.2 expression recovers after an initial 487 presumed decrease caused by the absence of JRS1, or that late nkx3.2 expression is less 488 important for shaping the developing jaw joint, and that other factors take over. At 6 dpf, we 489 could not detect any significant differences in nkx3.2 expression levels, supporting the former hypothesis. If this is the case, it suggests that it takes approximately one week for the 491 phenotype to be rescued once normal nkx3.2 expression is restored. It also indicates that 492 JRS1 is most important to the early jaw joint expression of nkx3.2, prior to 6 dpf, and may 493 contribute relatively less to later expression. Castellanos  This work highlights the importance of screening enhancer deletion mutants for phenotypic 515 effects at multiple developmental stages, as the dynamic nature of gene expression by 516 temporally-specific enhancers can mask or otherwise rescue phenotypes observed earlier in 517 development. We encourage future studies to assess enhancer activity and gene expression 518 to quantify these temporal dynamics, especially in known multi-enhancer systems. 519 520 In summary, we report the discovery and functional characterisation of the first described 521 Nkx3.2 enhancer, JRS1, with activity specific to the developing primary jaw joint. As this 522 enhancer is conserved in sequence and activity in gnathostomes including both bony and 523 cartilaginous fish, we conclude that it arose early in gnathostome evolution and may have 524 been one of the earliest novel cis-regulatory elements to facilitate the evolution of jaws from 525 the ancestral jawless state by driving the localisation of Nkx3.2 into the jaw joint and 526 contributing to the early jaw joint morphology. Further study of JRS1 and other as-yet 527 undiscovered Nkx3.2 enhancers will provide new insights into developmental regulatory 528 network responsible for the evolution of gnathostome jaws. 529 530

Ethics statement 532
All animal experimental procedures were approved by the local ethics committee for animal 533 research in Uppsala, Sweden (permit number C161/4 and 5.8.18-18096/2019). All 534 procedures for the experiments were performed in accordance with the animal welfare 535 guidelines of the Swedish National Board for Laboratory Animals. 536

Conserved synteny analysis 538
The synteny analysis was performed on the genomic regions containing the Nkx3.2 gene 539 (alternative names Nkx3-2 and Bapx1). Synteny data including upstream and downstream 540 genes from Nkx3.2 was extracted from the Ensembl database for several genomes: human, 541 Homo sapiens (GRCh38.p12); mouse, Mus musculus (GRCm38.p6); koala, Phascolarctos were modified for the T7 synthesis needs, and the guide core sequence were performed. 634 Reaction products were used as a template for the in vitro transcription (HiScribe T7 High 635 Yield RNA Synthesis Kit, England Biolabs) and purified. Cas9 mRNA was prepared from the p-636 T3Ts-nCAs9 plasmid (46757 Addgene) including the restriction enzyme digestion (XbaI, New 637 England Biolabs), in vitro transcription (mMESSAGE mMACHINE T3 Transcription Kit, Life 638 Technologies) and product purification. products were run on 1% agarose gels to determine their size. PCR product representing 653 enhancer deletion allele was 352 bp in length, and easily distinguishable from 790 bp wild-654 type allele. Four deletion-positive F0 zebrafish were outcrossed with Tg(fli1a:GFP) fish 655 (Lawson and Weinstein, 2002) to produce four F1 lines, which were raised to adulthood and genotyped by fin clip. One of the four heterozygous lines (deletion allele uu3731, designated 657 ΔJRS1) was selected for further analysis. These F1 nkx3.2 +/ Δ JRS1 fish were incrossed to produce 658 F2 nkx3.2 Δ JRS1/ Δ JRS1 fish for analysis. JRS1 deletion was confirmed by Sanger sequencing 659 (Eurofins Genomics) using the same PCR primers. GTCTCGGTGAGTTTGAGGGA-3'. Amplicon sizes for these target genes were 148bp and 699 187bp, respectively. A dissociation step was performed at the end of the analysis to verify 700 the specificity of the products, and standard curves were generated from pooled cDNA from 701 all samples in triplicate for each target gene to verify the efficiency of the primers (R 2 >0.99). 702 nkx3.2 expression was normalised to rpl13a levels and relative quantification of gene expression was calculated using the Pfaffl method (Pfaffl, 2001), displaying the fold 704 difference in heterozygous and homozygous mutants relative to wild-type, which was set to 705 1.0. FDR-adjusted P-values (Benjamini and Hochberg, 1995)  and bivariate plots of the PC axis were generated using the borealis package (Angelini, 2021). 723 FDR-adjusted P-values (Benjamini and Hochberg, 1995) are reported from pairwise analyses 724 of group means after accounting for allometric differences. 725