Foxc1 establishes enhancer accessibility for craniofacial cartilage differentiation

The specification of cartilage requires Sox9, a transcription factor with broad roles for organogenesis outside the skeletal system. How Sox9 and other factors gain access to cartilage-specific cis-regulatory regions during skeletal development was unknown. By analyzing chromatin accessibility during the differentiation of neural crest cells into chondrocytes of the zebrafish head, we find that cartilage-associated chromatin accessibility is dynamically established. Cartilage-associated regions that become accessible after neural crest migration are co-enriched for Sox9 and Fox transcription factor binding motifs. In zebrafish lacking Foxc1 paralogs, we find a global decrease in chromatin accessibility in chondrocytes, consistent with a later loss of dorsal facial cartilages. Zebrafish transgenesis assays confirm that many of these Foxc1-dependent elements function as enhancers with region- and stage-specific activity in facial cartilages. These results show that Foxc1 promotes chondrogenesis in the face by establishing chromatin accessibility at a number of cartilage-associated gene enhancers.


Introduction
Cartilage is the first skeletal type to be specified in the vertebrate body, providing important templates for later bone development and providing flexibility at joint surfaces and within the nose, ear, ribs, and larynx. The transcription factor Sox9 is essential for chondrogenic differentiation in all vertebrates examined (Bi et al., 1999;Lefebvre et al., 1997;Mori-Akiyama et al., 2003;Yan et al., 2005), yet it also has widespread roles outside the skeletal system, including the reproductive system, kidney, liver, and skin (Jo et al., 2014). How Sox9 is directed to a chondrogenic program in trunk mesoderm and cranial neural crest-derived cells (CNCCs) has remained unclear. Sox9 is known to directly bind to a number of cis-regulatory elements adjacent to chondrogenic genes, including Col2a1, Col10a1, and Acan (Askary et al., 2015;Dy et al., 2012;Lefebvre et al., 1997;Ohba et al., 2015). It is, however, dispensable for chromatin accessibility at these same elements (Liu et al., 2018), suggesting that other unknown factors may first open chromatin at chondrogenic enhancers for later activation by Sox9.
Forkhead-domain (Fox) family transcription factors are excellent candidates for establishing chromatin accessibility at chondrogenic enhancers. In the endoderm lineage, HNF3/FoxA binds closed chromatin at enhancers and makes these more accessible (Cirillo et al., 2002). Foxd3 has been similarly proposed to establish chromatin accessibility in the early neural crest lineage (Lukoseviciute et al., 2018). In mouse, loss of Foxc1 results in widespread cartilage and bone defects (Kume et al., 1998), including impaired tracheal and rib cartilages (Hong et al., 1999), loss of calvarial bone due to premature ossification (Rice et al., 2003;Sun et al., 2013;Vivatbutsiri et al., 2008), syngnathia (Inman et al., 2013), and disruption of endochondral bone maturation (Yoshida et al., 2015). In zebrafish, loss of both Foxc1 paralogs (foxc1a and foxc1b) results in severe reductions of dorsal cartilages of the upper face, which are preceded by reduced expression of several Sox9 targets, including col2a1a, acana, matn1, and matn4 (Xu et al., 2018). Chromatin immunoprecipitation followed by deep sequencing (ChIP-Seq) using a Sox9 antibody in dissected mouse rib and nose cartilage revealed enrichment of Fox binding motifs within Sox9bound cis-regulatory sequences near chondrogenic genes (Ohba et al., 2015). Here, we use profiling of chromatin accessibility in wild-type and mutant zebrafish facial cartilages and find that Foxc1 paralogs are required for accessibility and activity of a number of cartilage enhancers. These findings support a model in which Foxc1 promotes chondrogenesis by establishing selective accessibility of cartilage enhancers in CNCC-derived mesenchyme.

Results and discussion
Chromatin accessibility landscape in chondrocytes of the zebrafish face In order to identify potential cis-regulatory elements important for facial cartilage development, we performed a genome-wide analysis of chromatin accessibility in chondrocytes from 72 hr post-fertilization (hpf) zebrafish. We labeled chondrocytes by co-expression of sox10:Dsred and col2a1a:GFP transgenes ( Figure 1A) and isolated double-positive and control double-negative cells by fluorescence-activated cell sorting (FACS). We then subjected these cells to a modified 'micro' version of the assay for transposase-accessible chromatin followed by next-generation sequencing (mATACseq). In order to focus on potential distal cis-regulatory elements, we excluded accessible regions within 1 kb upstream or 0.5 kb downstream of transcription start sites. This analysis yielded 33,679 distal accessible elements, with 5736 elements over-enriched in chondrocytes and 8955 elements underenriched in chondrocytes ( Figure 1B,C). As confirmation of sequence quality, the cartilage-specific R2 enhancer of the Collagen Type II alpha 1 a gene (col2a1a) is enriched in chondrocyte-specific elements ( Figure 1-figure supplement 1A). De novo motif analysis of the top 2000 chondrocyteenriched regions using HOMER recovered motifs for Sox, Fox, Nfat, Zfx, and Nkx transcription factor families ( Figure 1D; Supplementary file 1A). Sox9 ChIP-Seq of mouse chondrocytes had previously eLife digest Animals with backbones (or vertebrates) have body shape determined, in part, by their skeletons. These emerge in the embryo in the form of cartilage structures that will then get replaced by bone during development. The neural crest is a group of embryonic cells that can become different tissues. In the head, it forms the cartilage scaffold for some of the facial bones and the base of the skull.
During this process, a protein called Sox9 is required for neural crest cells to morph into cartilage. This transcription factor binds to regulatory sequences in the genome to turn cartilage genes on. But Sox9 is also required to form non-cartilage tissues in organs such as the liver, lung, and kidneys. How, then, does Sox9 only turn on the genes required for cartilage formation in the embryonic face? This specificity can be controlled by which regulatory sequences Sox9 can physically access in a cell: controlling which regulatory sequences Sox9 can access determines which genes it can activate, and which type of tissue a cell will become. Xu, Yu et al. wanted to understand exactly how Sox9 switches on the genes that turn neural crest cells into facial cartilage. They studied the genomes of zebrafish embryos, which have a cartilaginous skeleton similar to other vertebrates, and found out which areas were accessible to transcription factors in the neural crest cells that became facial cartilage. Analyzing these regions suggested that sites where Sox9 could bind were often close to binding sites for another protein, called Foxc1. When zebrafish embryos were genetically modified to inactivate Foxc1 proteins, many of the regulatory sequences in cartilage failed to become accessible, and the cartilaginous skeleton did not form properly. These results support a model in which Foxc1 opens up the genomic regions that Sox9 needs to bind for cartilage to form, as opposed to the regions that Sox9 would bind to make different organ cell types.
The findings of Xu, Yu et al. uncover the stepwise process by which cartilage cells are made during development. Further research based on these results could allow scientists to develop new ways of replacing cartilage in degenerative conditions such as arthritis.
revealed Nfat and Fox motifs as the second and third most co-enriched with Sox motifs (Ohba et al., 2015), and the Nkx motif might reflect the role of Nkx3.2 in cartilage differentiation (Provot et al., 2006). Consensus sequences for Sox, Fox, and Nfat motifs were highly similar between zebrafish and mouse (Figure 1-figure supplement 1B), despite our zebrafish analysis focusing on all accessible regions in facial chondrocytes and mouse analysis focusing on only Sox9bound regions in rib chondrocytes. These striking motif similarities indicate strong conservation of the cartilage gene regulatory network between fish and mammals, and between the face and rib, and strongly suggest that many of the identified chondrocyte-specific accessible regions in zebrafish likely function as chondrocyte enhancers. In addition, Gene Ontology (GO) analysis of the nearest genes to the chondrocyte-enriched elements revealed cartilage development as one of the top six associated terms ( Figure 1E). We also recovered terms for neural crest cell migration and dorsal/ ventral pattern formation, likely reflecting retention of enhancer accessibility linked to the neural crest origins and later dorsoventral arch patterning of the precursors of facial cartilage.

Progressive establishment of cartilage-associated chromatin accessibility
We next investigated when cartilage-associated chromatin accessibility is established in relation to CNCC development. CNCCs are first specified at~10.5 hpf at the border of the neural keel, finish their migration into the pharyngeal arches by~20 hpf, and then show the first histological signs of cartilage development in the jaw at~52 hpf (Schilling and Kimmel, 1997). CNCC-derived arch  ectomesenchyme cells can be uniquely identified by co-expression of fli1a:GFP and sox10:Dsred transgenes at 36 and 48 hpf (Askary et al., 2017), stages just prior to cartilage differentiation ( Figure 2A). We therefore performed mATACseq on FACS-purified fli1a:GFP+; sox10:Dsred+ cells at these stages. 4915 of 14,623 elements gaining accessibility from 36 to 48 hpf were linked to GO terms including skeletal system development and cartilage development, and were enriched for Sox and Fox transcription factor binding motifs. 4323 elements with decreasing accessibility were linked to several GO terms related to cellular migration and enriched for motifs of neural crest-associated transcription factors (Nr2f, Lhx, Olig2, Hox, Ets), consistent with decommissioning of enhancers involved in earlier neural crest migration ( Figure 2-figure supplement 1).
To more specifically understand the timing of cartilage-associated chromatin accessibility, we next analyzed chondrocyte-specific accessible elements from our 72 hpf dataset for their accessibilities at 36 and 48 hpf ( Figure 2B). Of the 5,736 elements enriched in chondrocytes at 72 hpf, only 6% (356) had peak accessibility at 36 hpf with no further increase in accessibility by 48 hpf ('Group I'). In contrast, 48% (2741) displayed increased accessibility between 36 and 48 hpf ('Group II'), and 46% (2639) between 48 and 72 hpf ('Group III'). For Group I, de novo motif enrichment revealed predicted binding sites for members of the Nfat, Fox, Lhx, Nr2f, Meis, and Pax families, and significant GO terms included neural crest migration, cell migration in general, and dorsal/ventral pattern formation ( Figure 2C,D). Combined with the known involvement of Foxd3 (Montero-Balaguer et al., 2006;Stewart et al., 2006), Lhx6/8 (Denaxa et al., 2009), Nr2f1/2/5 ( Barske et al., 2018), Meis2 (Machon et al., 2015), and Pax9 (Nakatomi et al., 2010) in CNCC specification, migration, and dorsal-ventral arch patterning, many Group I elements likely represent retention of cis-regulatory elements involved in the earlier specification, migration, and regional patterning of CNCCs. Group II and Group III elements share many common predicted transcription factor binding motifs, including Sox dimer, Fox, Nfat, and Ap1 motifs previously described for mouse cartilage ( Figure 2E,G; Supplementary file 1B; Ohba et al., 2015). An Nkx motif was recovered only for Group II (p=10 À58 , 36% of targets), an Egr motif was enriched for Group II (p=10 À65 , 55% of targets) versus Group III (p=10 À20 , 16% of targets), and a Tead motif was enriched for Group III (p=10 À63 , 36% of targets) versus Group II (p=10 À20 , 3% of targets). GO analysis for linked genes also revealed terms related to skeletal system development (Group II) and cartilage development (Group III), as well as more general terms such as transcription and cell differentiation ( Figure 2F,H). We therefore conclude that the majority of chondrocyte-specific elements gain accessibility after pharyngeal arch formation and that transcription factor binding motifs change during cartilage differentiation. For example, enrichment of Nkx motifs in Group II elements might reflect the role of Nkx3.2 in limiting chondrocyte maturation (Provot et al., 2006) and promoting joint formation (Miller et al., 2003), while the preferential enrichment of the Tead motif, which is linked to growth-associated Hippo signaling (Ota and Sasaki, 2008), in Group III elements might reflect the later proliferative expansion of chondrocytes.

Requirement of Foxc1 for chromatin accessibility at a subset of cartilage elements
We had previously found that Foxc1 genes are essential for cartilage development in the upper face (Xu et al., 2018), and both our mATACseq analysis of zebrafish chondrocytes and published Sox9 ChIP-seq analysis in mouse (Ohba et al., 2015) reveals co-enrichment of Sox and Fox motifs in accessible regions near known cartilage genes. In zebrafish Foxc1 (foxc1aÀ/À; foxc1bÀ/À) mutants, cartilages of the upper/dorsal face fail to develop (Xu et al., 2018). In order to isolate the dorsal arch CNCC precursors affected in Foxc1 mutants, we used a pou3f3b:Gal4; UAS:nlsGFP (pou3f3b>GFP) dorsal CNCC transgenic line (Barske et al., 2020)    Of 22,323 elements at 48 hpf, 72% were similarly accessible, 10% more accessible in dorsal CNCCs, and 18% more accessible in pan-CNCCs (Figure 3-figure supplement 2B). A comparison of the 5,736 regions with specific accessibility in cartilage at 72 hpf revealed high correlation between dorsal-and pan-CNCCs at 36 hpf, with 96% displaying similar accessibility (r = 0.92, Figure 3-figure supplement 3A,C). By 48 hpf, however, we observed notable differences between accessibility of cartilage-specific elements, with 43% displaying greater accessibility in pan-CNCCs and only 0.2% displaying greater accessibility in dorsal CNCCs (r = 0.71, Figure 3-figure supplement 3B,D). The decreased accessibility of cartilage-specific elements in dorsal CNCCs at 48 hpf supports previous studies that dorsal chondrocytes develop later than other chondrocytes in the zebrafish face (Barske et al., 2016;Schilling and Kimmel, 1997).
Analysis of cartilage-associated elements in Foxc1 mutants revealed that 10% (120/1221) of elements that have established peak accessibility by 36 hpf (i.e. Group I) and 41% (636/1556) of elements that increase accessibility between 36 and 48 hpf (i.e. Group II) had reduced accessibility in Foxc1 mutants ( Figure 3E, Figure 3-figure supplement 4). De novo motif analysis of Foxc1-dependent and Foxc1-independent elements in Group I showed enrichment of Sox (p=10 À12 , p=10 À76 ), Tead (p=10 À14 , p=10 À20 ), Mtf (p=10 À12 , p=10 À22 ), and Ets (p=10 À12 , p=10 À61 ) motifs in both ( Figure 3F, Supplementary file 1C). Whereas several types of Fox motifs were uncovered in both, a closer analysis revealed enrichment of Foxa2 and Foxd3 motifs only in Foxc1-dependent elements (p=10 À45 , 38% of targets; p=10 À16 , 18% of targets), and a Foxo1 motif only in Foxc1-independent elements (p=10 À42 , 49% of targets). Similarly in Group II elements, a Foxa1 motif was enriched only in Foxc1-dependent elements (p=10 À59 , 38% of targets), and Foxh1 and Foxp1 motifs only in Foxc1independent elements (p=10 À22 , 22% of targets; p=10 À19 , 7% of targets) ( Figure 3G, Supplementary file 1C). Although Foxc1 motifs were not in the database used for motif predictions, the selective presence of Foxa1/2 motifs in Foxc1-dependent elements is consistent with previous reports that the consensus sequence for Foxc1-bound peaks is nearly identical to Foxa1/2 motifs (Wang et al., 2016). Sox motifs were similarly enriched in both Group I and Group II Foxc1-dependent and -independent elements, and Nfat and Zfx motifs in Group II Foxc1-dependent and -independent elements. An Insm (p=10 À39 , 30% of targets) motif was uncovered only in Foxc1independent Group I elements and an Nkx motif (p=10 À35 , 52% of targets) only in Foxc1-independent Group II elements. Whereas 33% (207/636) of Foxc1-dependent Group II elements had both Sox and Fox predicted motifs, only 16% (151/920) of Foxc1-independent Group II elements had both. These findings suggest that Foxc1-dependent and -independent cartilage elements may be commonly bound by Sox9 but likely differ in co-binding by Foxc1 and additional co-factors. For example, the presence of the Nkx motif only in Foxc1-independent elements suggests that it could be an alternative co-factor for Sox9, in line with the known roles for Nkx3.2 in chondrocyte biology (Provot et al., 2006).

Validation of Foxc1-dependent cartilage enhancers in zebrafish transgenesis assays
To verify whether Foxc1-dependent cartilage elements identified by mATACseq are chondrogenic enhancers, we tested the ability of individual elements in combination with an E1b minimal promoter to drive cartilage expression of green fluorescent protein (GFP) in zebrafish transgenic assays (Figure 4-figure supplement 1A). We tested 22 Foxc1-dependent Group II elements near 15 different genes, which included elements linked to genes with known cartilage function (ucmab, matn4, matn1, lect1, epyc, col9a1a, col9a3, sox10, acana, foxa3, mia) and others with unknown cartilage function (si:dkey33i1l.4, gas1b, lefty2, slc35d1a) (Supplementary file 2). We observed that 59% (13/ 22) of elements drove GFP expression in facial chondrocytes at 6 dpf. These included intronic elements within sox10, lect1, col9a3, col9a1a, and slc35d1a; distal 5 0 elements near ucmab, epyc, mia, acana, and matn4; a distal 3 0 element near gas1b; and a promoter-associated element for lect1. We  . Whereas an element in the first intron of sox10 drove uniform cartilage-specific expression, most elements drove expression in specific sub-regions or particular differentiation stages of cartilage. A distal 5 0 element of ucmab drove expression in chondrocytes of multiple joints in the zebrafish head, an intronic element of lect1 drove chondrocyte expression only in the jaw joint and hyomandibular-otic connection, a promoter-associated element of lect1 drove expression more strongly in the hyoid joint and hyomandibular-symplectic connection (though also more broadly in chondrocytes), a distal 5 0 element of acana drove restricted expression at the hyoid joint, and a distal 3 0 element of gas1b drove restricted expression at the Meckel's-Meckel's joint and connection between the hyomandibular cartilage and opercular bone. Reciprocally, elements associated with col9a3, epyc, mia, matn4, and col9a1a were expressed in chondrocytes but generally excluded from joint regions (particularly apparent at the hyoid joint and hyomandibular-symplectic connection). An intronic element for slc35d1a drove expression in both pharyngeal cartilage and muscle. Furthermore, we found that three enhancer transgenes with diverse expression patterns (broad sox10, joint-restricted ucmab, and joint-excluded epyc enhancers) all displayed reduced activity specifically in the dorsal cartilage regions affected in Foxc1 mutants ( Figure 4B). We also confirmed by in situ hybridization that four genes linked to lect1,col9a3,epyc) showed reduced expression in cartilageforming regions of the dorsal arches of Foxc1 mutants ( Figure 4C), similar to our previous results for col2a1a, matn1, matn4, and acana (Xu et al., 2018).
Testing of four Foxc1-independent Group II elements revealed two that drove cartilage expression (distal 5 0 elements of gas1b and matn4), one that drove cartilage, bone, and ligament expression (distal 5 0 element of sparc), and one with no activity (promoter-associated element of sox10) (Figure 4-figure supplement 1). Thus, the majority of tested Foxc1-dependent and -independent Group II elements are equally capable of driving cartilage expression. We also tested five Foxc1independent Group I elements and found that two showed arch CNCC expression at 36 hpf and minimal cartilage expression at 6 dpf (prrx1a and emx3 elements), one no expression at 36 hpf and ligament expression at 6 dpf (satb2), and two no expression at either stage (prrx1a, cd248a) (Figure 4-figure supplement 2). In contrast, we did not observe 36 hpf arch CNCC expression in any of the 16 Group II elements that drove cartilage expression at 6 dpf. As most of the Group I elements tested had decreased accessibility from 36 to 48 hpf, it seems likely that these elements represent either neural crest elements in the process of decommissioning (e.g. prrx1a and emx3) or elements prefiguring non-cartilage fates (e.g. satb2). Thus, most cartilage enhancers appear to gain accessibility after 36 hpf, with activity of late-opening cartilage enhancers in diverse locations, indicating that global cartilage expression patterns are achieved in part through the summation of enhancers with more restricted activity.

Conclusion
Our findings indicate that Foxc1 function is required for the accessibility of close to half of chondrocyte enhancers in the zebrafish face. Given expression and function of Foxc1 in diverse cartilages (e.g. limb, rib, tracheal) in mouse, it seems likely that Foxc1 has a similar function in chondrocyte enhancer accessibility throughout the body. The co-enrichment of predicted Sox and Fox binding sites in 33% of chondrocyte enhancers suggests a model in which Foxc1 promotes Sox9 binding to     these same enhancers by increasing chromatin accessibility, though this would require direct Sox9 binding studies to validate. However, not all Foxc1-dependent enhancers have predicted Fox binding motifs, which might suggest indirect regulation of chromatin accessibility in at least some cases. In addition, many enhancers do not appear to require Foxc1 activity. It is possible that other members of the Fox family compensate, such as Foxf1/2 in the facial midline (Xu et al., 2018) and Foxa2/ 3 during later hypertrophic maturation (Ionescu et al., 2012). Alternatively, there may be other cofactors that mediate chondrocyte enhancer accessibility and activation, such as Nkx3.2. We did not detect any obvious differences in the types of enhancer-proximal genes or the patterns and stages associated with Foxc1-dependent versus -independent enhancer activity. Further work will also be needed to understand the mechanism by which Foxc1 promotes chondrocyte enhancer accessibility, and the extent to which this reflects direct binding of Foxc1 to chondrogenic enhancers. Foxc1 lacks a chromatin modifying domain (Yoshida et al., 2015) and therefore would need to interact with a co-factor to directly open chromatin. Foxc1 could also act to maintain open chromatin, as shown for Foxa1 in the liver (Reizel et al., 2020). Given that both Foxc1 and Sox9 have expression in many tissues outside the skeletal system, it will also be important to determine whether additional factors help to further restrict their activity to chondrocyte enhancers within skeletogenic mesenchyme.  (Xu et al., 2018); Tg(sox10:Dsred) el110 and Tg (fli1a:EGFP) y1 (Askary et al., 2017); Tg(col2a1aBAC:GFP) el483 (Paul et al., 2016); and Tg(UAS: nlsGFP;a-crystallin:Cerulean) el609 and pou3f3b Gal4ff-el79 (Barske et al., 2020). For enhancer transgenic lines, we synthesized accessible elements with flanking attB4 and attB1 sequences using IDT gBlocks and cloned these into pDONR-P4-P1R using the Gateway Tol2kit (Invitrogen) to create p5E enhancer constructs (Kwan et al., 2007). We then combined p5E constructs with pME-E1b-GFP, p3E-polyA, and pDestTol2AB2 using LR clonase. Final DNA constructs were microinjected with transposase RNA (30 ng/ml each) into one cell stage zebrafish embryos. In most cases, multiple independent stable alleles per construct were analyzed in the F1 generation (Supplementary file 2).

Confocal imaging
In situ hybridization and transgenic embryos were imaged with a Zeiss LSM800 confocal microscope. To facilitate embryo collection at the 36 hpf time point, embryos were moved at 27 hpf to an incubator set at 22˚C to delay their development such that they reached 36 hpf the following morning. Cell dissociation and FACS were performed as previously described (Askary et al., 2017). Around 5,000 cells of each sample were centrifuged at 500 g for 20 m at 4˚C, and the pellet was suspended with 20 mL of lysis buffer (10 mM Tris-HCl [pH 7.4], 5 mM MgCl 2 , 10% DMF, 0.2% N-P40) by pipetting 6-10 times to release the nuclei without purification. The cell lysate was then mixed with 30 mL reaction buffer (10 mM Tris-HCl [pH 7.4], 5 mM MgCl 2 , 10% DMF, and homemade Tn5 Transposase) by vortexing for 5 s. The reaction was incubated at 37˚C for 20 min, followed by DNA purification using a Qiagen MiniElute kit. Purified DNA fragments were used to construct mATACseq libraries as previously described (Buenrostro et al., 2013) and sequenced using the NextSeq 500 platform (Illumina) with a minimum of 50 million paired-reads/sample. Two biological replicates of mATACseq experiments were performed for each condition.

Data analysis and statistics
The Encode analysis pipeline (https://github.com/ENCODE-DCC/chip-seq-pipeline) for ATACseq (Davis et al., 2018) was used with small modifications. The raw reads were trimmed to 37 bp and aligned to the zebrafish GRCz10 genome assembly by STAR aligner (Dobin et al., 2013). PCR duplicates, and the reads that aligned to 'blacklist regions' (Amemiya et al., 2019), were removed, and then peaks were called by model-based analysis of ChIP-Seq (MACS2) (Zhang et al., 2008), with p=10 À7 cutoff and disabled dynamic lambda option (-nolambda) for individual replicates. Only peaks common for two biological replicates (Irreproducibility Discovery Rate < 0.1) were kept for further analysis. For visualization, bigwig files were generated from duplicates-removed bam files with bedtools and bedgraphtobigwig, and normalization was based on total read numbers. Heatmaps were generated with deeptools2 (Ramírez et al., 2016) based on the normalized bigwig signal files. Individual genomic loci were examined by IGV (Broad Institute). For quantitatively comparing the accessibility of distal regulatory elements between two conditions, peaks from two conditions were first merged, and only distal regulatory regions (1 kb upstream or 0.5 kb downstream from the transcription start sites) were kept for comparison. Raw read counts from each replicate were computed using 'bedtools multicov' function (bedtools multicov -bams input.bam -bed peak.bed > count.txt). Raw read count matrices from four replicates (two replicates per condition) were analyzed using DESeq2 package (Love et al., 2014). FDR = 0.1 was used as the cut off to filter the peaks that were differentially accessible between two conditions. Volcano plots were generated using ggplot2 package in R. Fold change and adjusted p value were outputs from DESeq2 package in R, which is the same analysis behind the heatmaps in Figure 3. We used a Àlog10 adjusted p value greater than one as the cutoff for determining significantly changed peaks. HOMER (Heinz et al., 2010) was used to identify de novo motifs and their associated p values. David 6.8 GO analysis was performed on the web interface (https://david.ncifcrf.gov/) based on nearest neighbor genes for all differentially accessible elements. Pearson correlation was calculated using cor() function in R.

Materials and correspondence
Requests for material should be directed to J. Gage Crump (gcrump@usc.edu).

Additional files
Supplementary files . Supplementary file 1. Transcription factor binding motifs in cartilage-accessible elements. Detailed lists of predicted motifs for cartilage-accessible elements in 72 hpf controls (A), and based on timing of gain of accessibility (B), and dependence on Foxc1 function (C).
. Supplementary file 2. Summary of transgenic analysis of cartilage-accessible elements. Description of independent alleles for each enhancer transgenic line, genomic coordinates of elements tested, and activity in cartilage in 6 dpf zebrafish.
. Supplementary file 3. In situ probes. Details on primer sequences used to amplify probe regions for cloning, enyzmes used to linearize the probe plasmids, and the types of RNA polymerase used to synthesize RNA probes for in situ hybridization.
. Transparent reporting form

Data availability
Chromatin accessibility data have been deposited in GEO under accession number GSE157575.
The following dataset was generated: Author (