Effects of SOX9 dosage changes on RE chromatin accessibility

To assess the effect of SOX9 dosage changes on chromatin accessibility, we carried out the assay for transposase-accessible chromatin with sequencing (ATAC-seq) on SOX9-tagged CNCCs with six different SOX9 dosages achieved by varied dTAG^V-1 concentrations (Fig. ​(Fig.1d),1d), as well as on WT CNCCs treated with either dimethylsulfoxide or the highest dTAG^V-1 concentration (500nM). Principal component (PC) analysis on ATAC counts per million (CPM) values of the 151,457 reproducible peak regions (which are candidate REs and are herein referred to as REs) revealed a batch-independent dTAG^V-1 (and thus SOX9 dosage titration) effect in PC space (Extended Data Fig. ​Fig.2a2a).

Open in a separate window

Extended Data Fig. 2

Effects of SOX9 dosage changes on chromatin accessibility.

(a) Principal component analysis of ATAC-seq counts per million (CPM) of all 151,457 REs across all CNCC samples. Shapes indicate the dTAGV-1 concentration treated for 48h. Colors indicate the combination of hESC line from which CNCCs were derived and differentiation batch (S9c1/2=SOX9-tagged clone1/2). Arrow indicates the SOX9 dosage effect. (b) Volcano plot of 500nM dTAGV-1 treatment on two SOX9-tagged (left) or two WT (right) CNCC differentiation replicates for all 151,457 REs. -log10(p-value) (y-axis) represents the unadjusted two-sided p-value from DESeq2, point color represents Benjamini-Hochberg adjusted p-value. (c) Distributions of fold-changes versus full (100) SOX9 dosage for all REs for which SOX9 dosage explains a significant (5% FDR; red, 16,538; blue, 27,334) or nonsignificant (grey, n=107,585) amount of variance (likelihood ratio test, LRT), stratified by the direction of change in full SOX9 depletion. Boxplot center represents median, box bounds represent 25th and 75th percentiles, and whiskers represent 5th and 95th percnetiles. (d) Same fold-change values as in (c) for a random subset (n=10,000) of significant REs, plotted as a heatmap and clustered by row based on Kendall distance. (e) Example ATAC-seq browser tracks of individual RE accessibility at different SOX9 dosages (y-axis, normalized coverage in 10bp bins), averaged across all replicates at each dosage.

WT CNCCs treated with 500nM dTAG^V-1 clustered with untreated CNCCs in PC space and had no significantly (5% false discovery rate (FDR)) changed REs, as compared to 6,169 changed REs from two SOX9-tagged replicates treated with 500nM dTAG^V-1 (Extended Data Fig. ​Fig.2b),2b), indicating minimal off-target effects. Plotting each SOX9-tagged sample’s loading in this PC direction versus SOX9 dosage revealed a nonlinear relationship, indicated by a lower Aikake information criterion (AIC) for a nonlinear Hill equation than a linear function (Fig. ​(Fig.2a).2a). These results also suggest a largely monotonic effect of SOX9 dosage on individual RE accessibility, which we confirmed by pairwise comparisons between all reduced SOX9 dosages and full dosage (Extended Data Fig. 2c,d).

Open in a separate window

Fig. 2

Most SOX9-dependent REs are buffered in their response to SOX9 dosage changes, with a sensitive subset.

a, Loadings from PC analysis of ATAC-seq CPM of all 151,457 REs across all CNCC samples (see Extended Data Fig. ​Fig.2a),2a), corrected for differentiation batch and plotted as a function of estimated relative SOX9 dosage (shown as percentage relative to no dTAG^V-1). Black line represents Hill equation fit. Points are biological (differentiation or clone) replicates. b, Example ATAC-seq browser tracks of individual RE accessibility at different SOX9 dosages (yaxis, normalized coverage in 10-bp bins, identical range in all browser tracks), averaged across all replicates at each dosage. c, Schematic of approach to model nonlinearity of SOX9-dependent REs. a and b in linear model refer to slope and intercept, respectively. a and b in Hill equation refer to ED₅₀ and Hill exponent, respectively. d, Illustration of different ED₅₀ and Hill exponent values on theoretical Hill equation curves. e, Individual REs from b with replicates, fitted by Hill equation (black line). f, Histogram of ΔAIC of all 35,713 SOX9-dependent REs. Red line indicates ΔAIC=2. g, Scatterplot of ED₅₀ and Hill exponent of 23,414 SOX9-dependent REs with good fit (P<0.05 for both parameters).

Individual REs showed distinct responses to SOX9 dosage, with some showing constant decreases in accessibility correlated with SOX9 dosage, and others showing buffered responses (that is, minimal accessibility changes until SOX9 dosage is greatly reduced; Fig. 2b,e), with a similar variety for upregulated REs (Extended Data Fig. ​Fig.2e).2e). Previous studies have observed similar variation in dose–response curves, either among select targets of the immune TF NF-κB⁴¹ or in cytokine-induced signaling⁴². To model these responses, we first defined all SOX9-dependent REs as those responding significantly (5% FDR) to full depletion of SOX9 over 48h, using all seven replicates in either condition. This yielded 35,713 REs, of which 20,346 decreased and 15,367 increased in accessibility following SOX9 depletion. Most SOX9-dependent REs were better fitted by the Hill equation than a linear model (median linear minus Hill ΔAIC of 5.2, 73.9% with ΔAIC>2; Fig. ​Fig.2f).2f). To allow direct comparisons among REs, we fitted the Hill equation to all SOX9-dependent REs for all subsequent analyses.

The Hill equation yields two key parameters: the empirical dose 50 (ED₅₀) representing the SOX9 dosage at which the RE reaches half of its maximal levels, and the Hill exponent, which indicates how switch-like the RE response is (Fig. ​(Fig.2d).2d). In this study, we define sensitivity (or its inverse, buffering) based on the RE response to decreasing SOX9 dosage from 100%. Higher ED₅₀ values (at constant Hill exponent) indicate increased sensitivity, while higher Hill exponents (at constant ED₅₀) indicate decreased sensitivity. Both values varied between REs (Fig. ​(Fig.2g),2g), but the ED₅₀ was substantially more correlated with an alternative measure of sensitivity than the Hill exponent (Extended Data Fig. ​Fig.3;3; Spearman ρ of −0.96 and −0.45 for ED₅₀ and Hill, respectively), indicating that it is the main determinant of sensitivity/buffering. Of all SOX9-dependent REs, 26,026 (73%) have ED₅₀<30 (buffered), 5,276 (14%) of REs have ED₅₀ between 30 and 40 (moderately sensitive), and 4,411 (13%) of REs have ED₅₀>40 (highly sensitive; Supplementary Table 1). The proportion of downregulated or upregulated REs in each of these groups is consistent (about 68% downregulated versus about 32% upregulated). Together, these results indicate a range in RE responses to SOX9 dosage, with most SOX9-dependent REs buffered against changes in SOX9 dosage but some showing more sensitive responses.

Open in a separate window

Extended Data Fig. 3

Fitted ED₅₀ values are well-correlated with orthogonally calculated measures of sensitivity.

For all SOX9-dependent REs with good Hill equation fits (p<0.05 for both ED₅₀ and Hill exponent), correlation between either ED₅₀ (left, Spearman ρ −0.961) or Hill exponent (right, Spearman ρ −0.457) and buffering index calculated at 50% SOX9 dosage. See Methods for details on calculation of buffering index – 0 means no buffering (effect of 100 to 50% SOX9 dosage on RE accessibility is 50% of effect of 100 to 0% SOX9 dosage), 100 means full buffering (no effect of 100 to 50% SOX9 dosage on RE accessibility, but substantial effect of 100 to 0% SOX9 dosage).

Features affecting RE sensitivity to SOX9 dosage

We next sought to identify regulatory features associated with RE sensitivity to SOX9 dosage. For the remainder of this paper, we use a bootstrapping approach when comparing ED₅₀ values between groups of REs/genes to incorporate fitting uncertainty (n=200 bootstraps; Methods). We reasoned that the SOX9-dependent REs comprised a mix of direct effects of SOX9 regulation and indirect effects acting through other TFs. Direct SOX9 effects should arise rapidly after full SOX9 depletion whereas indirect effects should be delayed. We therefore carried out ATAC-seq 3h after 500nM dTAG^V-1 treatment of SOX9-tagged CNCCs (yielding full SOX9 depletion within 1h; Extended Data Fig. ​Fig.4a).4a). Of the 35,713 48-h SOX9-dependent REs, 9,279 showed significant (5% FDR) accessibility changes at 3h, of which almost all (96.3%) were decreases (Fig. ​(Fig.3a).3a). Relative to delayed and/or upregulated REs, rapidly downregulated REs were substantially more likely to harbor the SOX9 palindrome sequence motif (Fig. ​(Fig.3b),3b), as well as SOX9 binding as assessed by V5 chromatin immunoprecipitation and sequencing (ChIP–seq; Extended Data Fig. ​Fig.4b).4b). These results are consistent with SOX9 acting as a direct activator of REs at rapidly downregulated sites (n=9,279) and indirectly regulating delayed and/or upregulated sites (n=26,434). Compared to delayed and/or upregulated REs, rapidly downregulated sites were substantially more sensitive to SOX9 dosage (that is, had higher ED₅₀ values; Fig. 3c,d).

Open in a separate window

Fig. 3

Features affecting sensitivity of the RE response to SOX9 dosage.

a, Scatterplots of full SOX9 depletion effect on chromatin accessibility at 48h (xaxis) versus 3h (yaxis) for all 48-h SOX9-dependent REs. b, Barplot indicating fraction of REs containing a SOX9 palindrome (3–5-bp spacing), stratified by time delay of response. c, ED₅₀ of REs stratified by time delay of response. d, Median dosage curves based on the median fitted ED₅₀ and Hill exponents for each group. a.u., arbitrary units. N in c,d indicated in b. e,f, ED₅₀ of rapid-down REs (likely direct SOX9 targets) stratified by SOX9 motif type, with motif position weight matrices on the left (e), or overlap with ChIP–seq peaks for TWIST1 (x axis) TFAP2A (color) and NR2F1 (shape) (f). n of groups from left to right in e: 2,263, 1,315, 360, 5,221. n of groups from left to right in f: 1,874, 736, 985, 1,805, 774, 602, 296, 2,087. g, Median dosage curves as in d for the indicated combinations of SOX9 motif and binding of other TFs. Points and error bars in (c,e,f) represent median and 95% confidence intervals as computed by 200 bootstraps (see Methods).

Open in a separate window

Extended Data Fig. 4

Rapidly responding REs are likely direct targets of SOX9.

(a) mNeonGreen fluorescence intensity (at least 5,000 cells per histogram) in SOX9-tagged or WT CNCCs with the times of treatment by dTAGV-1. (b) V5 ChIP-seq signal from CNCCs with V5-tagged SOX9 present (‘SOX9-tagged V5’) or absent (‘WT’) plotted over sets of SOX9-dependent REs as defined in Fig. ​Fig.3a.3a. (c) Hill exponent of rapid down REs stratified by SOX9 motif type, with motif position weight matrices as in Fig. ​Fig.3e.3e. N of groups from left to right: 2,263, 1,315, 360, 5,221. (d) For the SOX9 palindrome motif at with a 0–10bp spacing between the inverted repeats (y-axis), rapid down REs were stratified on the basis of that motif match. N of blue color groups from bottom to top: 448, 1,249, 1,010, 1,740, 2,628, 1,593, 1,157, 931, 839, 1,249, 1,233. N of red color groups from bottom to top: 7,532, 7,933, 8,200, 7,410, 6,458, 7,568, 8,034, 8,289, 8,371,. Points and error bars represent median and 95% confidence intervals as computed by bootstrap (see Methods).

While direct versus indirect regulation is one contributor to RE sensitivity to SOX9 dosage, there is additional variation in sensitivity among the 9,279 direct SOX9 target sites: 4,266 (46%) have ED₅₀<30 (buffered), 2,760 (30%) have ED₅₀ between 30 and 40 (moderately sensitive), and 2,253 (24%) have ED₅₀>40 (highly sensitive). We therefore sought to identify additional features associated with variation in sensitivity among all direct SOX9 target sites. REs containing the full SOX9 palindrome with spacing of 3–5base pairs (bp) were more sensitive than sites containing either one or several partial palindromes, with REs containing no detected motif being the least sensitive (Fig. ​(Fig.3e).3e). The 3–5-bp SOX9 palindrome was also associated with a modest increase in the Hill exponent, consistent with the palindrome’s reported requirement for cooperative SOX9 binding⁴³ (Extended Data Fig. 4c,d). Thus, motif type and resultant SOX9 binding mode modulate RE sensitivity to SOX9 dosage among its direct targets. Among direct SOX9 targets, those with larger effects of SOX9 depletion were most sensitive to SOX9 dosage (Extended Data Fig. ​Fig.5a5a).

Open in a separate window

Extended Data Fig. 5

Additional features affecting sensitivity of the RE response to SOX9 dosage changes among direct SOX9 targets.

ED₅₀ of rapid down SOX9-dependent REs, stratified by (a) magnitude of change in response to full SOX9 depletion, (b) presence of Coordinator (TWIST1), NR2F1, or TFAP2A sequence motif matches, baseline levels of (c) H3K27ac or (d) chromatin accessibility, or (e) the combination of SOX9 motif type and TWIST1/TFAP2A binding by ChIP-seq. For (a), (c), and (d), higher deciles mean higher values, and N of each decile is 928 (except for deciles 1 and 6 for which N=927). For (b), N of groups from left to right: 1,539, 879, 1,735, 1,200, 1,693, 718, 952, 563. For (e), N of red circle groups from left to right: 578, 353, 87, 1,603; N of blue circle groups from left to right: 660, 389, 110, 1,712; N of red triangle groups from left to right: 454, 208, 54, 650; N of blue triangle groups from left to right: 691, 375, 107, 1,248. Points and error bars represent median and 95% confidence intervals as computed by bootstrap (see Methods).

We assessed whether additional factors beyond SOX9 motif type could modulate RE sensitivity to SOX9 dosage among its direct targets. We focused on binding by other TFs, specifically TWIST1, NR2F1 and TFAP2A, as they have well-known roles in CNCCs and their binding in hESC-derived CNCCs has previously been characterized^29,38. Binding of other TFs at SOX9 direct target sites substantially decreased RE sensitivity to SOX9 dosage; the strongest effects were seen for TWIST1 and TFAP2A, with minor effects of NR2F1 at TWIST1- and TFAP2A-bound REs (Fig. ​(Fig.3f).3f). We replicated this result using TF sequence motifs (Extended Data Fig. ​Fig.5b).5b). Baseline levels of both the active histone mark H3K27ac and chromatin accessibility were negatively correlated with sensitivity (Extended Data Fig. 5c,d). REs containing the SOX9 palindrome motif also unbound by other TFs were most sensitive to SOX9 dosage (Fig. ​(Fig.3g3g and Extended Data Fig. ​Fig.5e).5e). Together, these results indicate that at least three features independently contribute to variation in RE sensitivity to SOX9 dosage: direct versus indirect regulation by SOX9; among directly regulated SOX9 targets, the type of SOX9 motif and resulting binding mode; and binding of other key CNCC TFs. Thus, REs for which SOX9 is likely the primary TF directly driving accessibility are most sensitive to SOX9 dosage.

To understand mechanisms underlying buffering against changes in SOX9 dosage among its direct targets, we carried out ChIP–seq of SOX9 (using the V5 tag) and TWIST1 in SOX9-tagged CNCCs with four different SOX9 concentrations achieved using dTAG^V-1. We grouped direct SOX9 targets into three bins based on ATAC-seq sensitivity to SOX9 dosage and plotted their ChIP–seq fold changes at each SOX9 concentration versus 100%. REs most sensitive to SOX9 dosage in accessibility are enriched for the SOX9 palindrome motif and have a lower fraction and level of TWIST1 binding (Extended Data Fig. 6a,c). These REs were also most sensitive in SOX9 binding, whereas the accessibility-buffered REs were also buffered in SOX9 binding. At about 50% SOX9 dosage, buffered sites retained nearly unperturbed levels of SOX9 binding (Extended Data Fig. ​Fig.6b).6b). TWIST1 binding showed similar responses to SOX9 dosage (Extended Data Fig. ​Fig.6d),6d), with no increase at partial SOX9 dosage, as might have been expected by compensatory buffering. Notably, stronger SOX9 perturbation (<20%) resulted in diminished TWIST1 binding at both buffered and sensitive sites (Extended Data Fig. ​Fig.6d),6d), consistent with a model in which SOX9 and TWIST1 binding is co-dependent at a subset of direct SOX9 target sites, but their synergistic function buffers co-regulated REs against small TF dosage changes (Extended Data Fig. ​Fig.6e6e).

Open in a separate window

Extended Data Fig. 6

Binding of SOX9 and TWIST1 in response to SOX9 dosage changes.

(a) Fractions of rapidly downregulated REs with SOX9 palindrome match (top) or TWIST1 binding (bottom), stratified by ATAC-seq sensitivity (ED₅₀) to SOX9 dosage (colors). N of groups by color: red, 823; green, 2,473; blue, 5983. (b) Distributions of ATAC-seq (left) or SOX9 (V5) ChIP-seq (right) fold-changes vs full SOX9 dosage at each concentration for rapidly downregulated REs with SOX9-dependent ChIP-seq signal (DESeq2 log2FoldChange<0 in 0 vs 100 SOX9 dosage), stratified based on ATAC-seq sensitivity to SOX9 dosage (colors). N of groups by color: red, 106; green, 391; blue, 666. P-values from two-sided Kruskal-Wallis tests comparing each instance of three groups. (c) Unperturbed TWIST1 ChIP-seq signal at TWIST1-bound, rapidly downregulated REs stratified based on ATAC-seq sensitivity to SOX9 dosage (colors). N of groups by color: red, 513; green, 1318; blue, 1956. Boxplot center represents median, box bounds represent 25th and 75th percentiles, and whiskers represent 5th and 95th percentiles. (d) TWIST1 ChIP-seq fold-changes vs full SOX9 dosage at same groups of REs as in (c). In (b) and (d), points and error bars represent median and 25th and 75th percentiles of distribution. (e) Models for RE buffering by synergistic TF functions. P-values from two-sided Kruskal-Wallis tests comparing each instance of three groups.

Effects of SOX9 dosage on gene expression

We next assessed the gene expression response to SOX9 dosage by RNA sequencing (RNA-seq) analysis of the same SOX9 dosage series. This response was largely monotonic (Extended Data Fig. ​Fig.7a)7a) and varied in its shape between individual genes (Fig. ​(Fig.4a4a and Extended Data Fig. ​Fig.7b).7b). SOX9 dosage changes had overall nonlinear effects on gene expression in PC space (Extended Data Fig. ​Fig.7c),7c), and most (70.3%) of the 1,232 SOX9-dependent genes (of which 688 decreased and 544 increased following full depletion, 5% FDR) were better fitted by a Hill than a linear equation (Fig. ​(Fig.4b4b and Supplementary Table 2), with variability in the ED₅₀ and Hill exponents (Fig. ​(Fig.4c).4c). Of all SOX9-dependent genes, 76% have ED₅₀<30 (buffered), 12% have ED₅₀ between 30 and 40 (moderately sensitive), and 12% have ED₅₀>40 (highly sensitive; Supplementary Table 2). Thus, like REs, most genes are buffered against SOX9 dosage, and a minor subset is sensitive.

Open in a separate window

Fig. 4

RE dose–response curves partially predict the shape of gene dose–response curves.

a, Examples of genes with sensitive (left) or buffered (right) responses to SOX9 dosage changes, as assessed by RNA-seq. Black lines represent Hill equation fits. b, Histogram of ΔAIC of all 1,232 SOX9-dependent genes, calculated as in Fig. ​Fig.2c.2c. Red line indicates ΔAIC=2. c, ED₅₀ and Hill exponent of 832 SOX9-dependent genes with good fit (P<0.05 for both parameters). d, Schematic of approach to predict RNA level changes based on ABC contribution scores and fold changes of REs at each SOX9 dosage. e, Distributions of observed (left) or predicted (right) fold changes versus full SOX9 dosage at each concentration, normalized to full-depletion fold change and stratified by ED₅₀ (colors). Only genes transcriptionally downregulated by 24-h full SOX9 depletion (assessed by metabolic labeling, SLAM-seq) are analyzed. n of groups by color: red, 66; green, 77; blue, 57. Points and error bars represent median and 25th and 75th percentiles of distribution. Indicated P values are from two-sided Kruskal–Wallis test comparing distributions of the three colored groups.

Open in a separate window

Extended Data Fig. 7

Effects of SOX9 dosage changes on gene expression.

(a) Heatmap of fold-changes versus full (100) SOX9 dosage for all all genes for which SOX9 dosage explains a significant (5% FDR) amount of variance (likelihood ratio test), clustered by row based on Kendall distance (b) Examples of genes upregulated in response to SOX9 depletion with buffered (left) or sensitive (right) responses. Black and blue lines represent Hill and linear fits, respectively (c) Principal component analysis of RNA-seq counts per million (CPM) across all SOX9-tagged and WT CNCC samples. Shapes indicate the dTAGV-1 concentration treated for 48h. Colors indicate the combination of hESC line from which CNCCs were derived and differentiation batch (S9c1/2=SOX9-tagged clone1/2). Arrow indicates the SOX9 dosage effect.

We investigated whether RE responses can predict the responses of their cognate target genes. We focused on the subset of SOX9-dependent genes showing changes in nascent transcription in response to 3-h or 24-h SOX9 depletion, assessed by thiol (SH)-linked alkylation for the metabolic sequencing of RNA (SLAM-seq)⁴⁴. Of the 1,232 48-h SOX9-dependent genes, 122 (62 down, 60 up) responded significantly (10% FDR) at 3h, and 395 (206 down, 189 up) responded at 24h (Supplementary Table 2). Effect sizes at 24h were correlated with, albeit larger than, those at 3h (Extended Data Fig. ​Fig.8a).8a). Known direct SOX9 targets such as COL2A1 (ref.⁴⁵) responded at 24h but not 3h (Extended Data Fig. ​Fig.8b),8b), suggesting a time lag between chromatin and transcriptional effects. Accordingly, we sought to predict the gene expression responses for genes responding transcriptionally to 24h of SOX9 depletion.

Open in a separate window

Extended Data Fig. 8

Measuring and predicting response of transcriptionally regulated SOX9 target genes.

(a) Scatterplot of effects of full SOX9 depletion for 3h (x-axis) versus 24h (y-axis) on nascent transcription, as assayed by SLAM-seq, for all SOX9-dependent genes (that is responding to full SOX9 depletion for 48h in RNA-seq). Y=x line in red. (b) Effects of full SOX9 depletion on transcription of COL2A1, a known direct target of SOX9. Points represent SLAM-seq counts from biological replicates, adjusted p-value calculated by DESeq2. (c, d) For all REs responding to full depletion of SOX9 at 48h, the effects of 3h (c) or 24h (d) full SOX9 depletion on chromatin accessibility (ATAC-seq, x-axis) or H3K27ac levels (ChIP-seq, y-axis) is plotted. (e) Distributions of observed (left) or predicted (right) fold-changes vs full SOX9 dosage at each concentration, stratified based on direction of transcriptional response to full SOX9 depletion (colors). N of groups by color: red, 184; grey, 10,339; blue, 197. Points and error bars represent median and 25th and 75th percentiles of distribution. ** p=2.1e-12, *** p<2.2e-16, two-sided Kruskal-Wallis test comparing the three groups. (f) Examples of predictions for a buffered (left) or sensitive (right) gene. (g) Median absolute deviation between observed and predicted dosage response curves for transcriptionally downregulated genes, stratified by number of SOX9-downregulated REs within 100kb of TSS. Spearman rho for correlation and associated two-sided p-value are shown. N of groups from left to right: 227, 238, 173, 162, 105, 77, 55, 35, 32, 26, 49. Boxplot center represents median, box bounds represent 25th and 75th percentiles, and whiskers represent 5th and 95th percentiles. (h) Median ED₅₀ of REs within 100kb of the TSS of transcriptionally downregulated genes stratified by sensitivity to SOX9 dosage. N of groups from left to right: 27, 27, 27, 28. Points and error bars represent median and 95% confidence intervals as computed by bootstrap (see Methods). All reported correlation coefficients in this figure are Spearman rho values.

The activity-by-contact (ABC) model predicts RE–gene connections by computing the contribution of each RE to transcription (ABC score) as its ‘activity’ (combination of accessibility and H3K27ac) divided by its contact (estimated by chromatin conformation capture or a genomic distance–power law function), normalized to the contributions of other REs⁴⁶. We used ABC to predict the effect of multiple REs changing in ‘activity’ at each SOX9 dosage. Although ‘activity’ includes H3K27ac levels, effects of full SOX9 depletion on accessibility and H3K27ac were highly correlated (Extended Data Fig. 8c,d). Thus, the fold change in the expression of a gene is predicted as the average fold change in accessibility at all nearby REs, weighted by the contribution of each RE (ABC score; Fig. ​Fig.4d4d).

We first assessed prediction of directionality of expression changes, comparing observed and predicted responses for transcriptionally upregulated or downregulated genes as well as SOX9-independent genes. Predicted responses significantly stratified these genes in the same manner as observed responses, most accurately for downregulated genes (Extended Data Fig. ​Fig.8e).8e). We thus focused on predicting differences in SOX9 dosage sensitivity among downregulated genes. Predicted responses separated SOX9-downregulated genes binned by their observed sensitivity to SOX9 (Fig. ​(Fig.4e),4e), although to a lesser extent than the observed responses. Inspection of top genes indicated accurate prediction of both buffered (TENT5B) or sensitive (SOX5) responses (Extended Data Fig. ​Fig.8f).8f). Genes with less accurate predictions have fewer SOX9-dependent REs in their vicinity (Extended Data Fig. ​Fig.8g),8g), suggesting that the model does not fully capture complex cis-regulatory landscapes with multiple inputs. Nevertheless, genes with more sensitive nearby REs were more sensitive to SOX9 dosage (Extended Data Fig. ​Fig.8h).8h). Together, these results indicate that broadly, the RE response to SOX9 dosage (sensitive or buffered) translates into the expression response of cognate genes based on the contribution of the RE to the transcription of that gene.

The pro-chondrogenic program is sensitized to SOX9 dosage

We next sought to assess the impact of SOX9-sensitive genes on cellular phenotypes, focusing on chondrogenic differentiation potential as SOX9 functions in both entry into and continuation of chondrogenesis⁴⁷. Genes with both cartilage development functions and increased expression during in vitro chondrogenesis (‘pro-chondrogenic genes’) showed substantially higher ED₅₀ values than other gene groups (Fig. ​(Fig.5a).5a). Examples include the collagen-encoding genes COL11A1 (highly sensitive) and COL2A1 (moderately sensitive; Fig. ​Fig.5b),5b), as well as genes encoding other transcriptional regulators such as SOX5 (Extended Data Fig. ​Fig.8f).8f). SOX9-upregulated genes did not yield a similar pattern of sensitivity (Extended Data Fig. ​Fig.9a),9a), suggesting that pro-chondrogenic functions of SOX9 may be especially dosage sensitive. Gene Ontology analysis of all moderately and highly sensitive (ED₅₀>30) genes revealed enrichment for cartilage condensation function as well as additional pathways with important roles in CNCCs, such as transforming growth factor beta and bone morphogenetic protein, but also neuronal/glial-related pathways unlikely to have important functions in mesenchymal CNCCs (Supplementary Table 3).

Open in a separate window

Fig. 5

The pro-chondrogenic program is sensitive to SOX9 dosage.

a, ED₅₀ of SOX9-downregulated genes stratified by presence in the ‘Cartilage development’ Gene Ontology (GO) category (xaxis), and expression change in chondrocytes compared to CNCCs (color, data from ref.²⁹). n of groups from left to right: 157, 269, 241, 6, 4, 11. Points and error bars represent median and 95% confidence intervals as computed by 200 bootstraps (see Methods). b, Examples of known pro-chondrogenic genes that are moderately or highly sensitive to SOX9 dosage. c, Schematic of approach to titrate SOX9 dosage during 21-day chondrogenic differentiation. d, Sulfated glycosaminoglycan (sGAG, representative of mature cartilage) content at day 21 of chondrogenesis as a function of SOX9 dosage as estimated in Extended Data Fig. ​Fig.9b.9b. Black curve represents Hill equation fit. e, Median dosage curves based on fitted ED₅₀ and Hill exponents for all REs and genes, pro-chondrogenic genes (purple, labeled group in a) and sGAG content (from d).

Open in a separate window

Extended Data Fig. 9

Effects of SOX9 dosage on chondrogenesis.

(a) ED₅₀ of SOX9-upregulated genes stratified by presence in the ‘Cartilage development’ Gene Ontology (GO) category (x-axis), and expression change in chondrocytes compared to CNCCs (color, data from Long et al.²⁹). N of groups from left to right: 94, 217, 204, 3, 9, 17. Points and error bars in (a,e) represent median and 95% confidence intervals as computed by bootstrap (see Methods). (b) Fluorescence intensity at day 10 (red) or 21 (blue) of chondrogenesis in SOX9-tagged chondrocytes as a function of dTAGV-1 concentration. gMFI, geometric mean. (c) Sulfated glycosaminoglycan (sGAG, representative of mature cartilage) at day 21 of chondrogenesis in WT CNCCs treated with DMSO or 500nM dTAGV-1 (N=5 for each group). Bars represent mean, p-value from two-sided T-test. (d) Fitted ED₅₀ values for all SOX9-downreguated REs (N=20,346) and other genes (N=688), pro-chondrogenic genes (N=11), and sGAG content (experiment depicted in Fig. ​Fig.5d5d with N=6 replicates per SOX9 concentration). Points (median) and error bars (95% confidence) for genes and REs computed by bootstrap, ED₅₀ point estimate and 95% confidence interval for sGAG content estimated from Hill equation model fit. (e) ED₅₀ by craniofacial disorder association for genes upregulated upon SOX9 depletion. Gene-craniofacial disorder associations determined as in Fig. ​Fig.6a.6a. N of groups from left to right: 508, 20, 9, 5. Points and error bars represent median and 95% confidence intervals as computed by bootstrap (see Methods).

To test whether increased sensitivity of pro-chondrogenic genes results in increased sensitivity of chondrogenesis, we titrated SOX9 dosage to five distinct levels both before and during 21-day differentiation of CNCCs to chondrocytes²⁹ (Fig. ​(Fig.5c5c and Extended Data Fig. ​Fig.9b).9b). To quantify functional chondrogenesis, we measured total levels of sulfated glycosaminoglycans (sGAGs), linear polysaccharides that mark extracellular matrix of mature cartilage, using a colorimetric assay. This revealed a nonlinear relationship between SOX9 dosage and functional chondrogenesis (Fig. ​(Fig.5d),5d), with no effect of the highest dTAG^V-1 concentration (500nM) on WT CNCC chondrogenesis, indicating minimal off-target effects (Extended Data Fig. ​Fig.9c).9c). The SOX9 dosage–sGAG curve more closely matched the curve for pro-chondrogenic genes than for other genes or REs (Fig. ​(Fig.5e5e and Extended Data Fig. ​Fig.9d).9d). Thus, in vitro chondrogenesis is sensitized to SOX9 dosage, more so than most genes or REs, at least partly owing to the heightened sensitivity of important pro-chondrogenic genes.

Statistics and reproducibility

No statistical method was used to predetermine sample size. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment. Facial scans were excluded if participants were laughing, crying or otherwise emoting or if the non-rigid registration failed. Facial scans from participants with missing covariate information (for example, age or sex) were additionally removed. Two RNA-seq samples (WT_R8_5e-7M and S9CC47_R6_5e-7M) were identified as extreme outliers in initial PC analysis. This was confirmed to be due to a library quality issue, and so these samples were excluded before any further downstream analyses. For modeling of individual REs/genes as a function of SOX9 dosage, outlier samples, defined as z-score greater than 3, were removed from the analysis of that RE/gene.

Cell culture

Female H9 (WA09; RRID: CVCL_9773) hESCs were obtained from WiCell and cultured in either mTeSR1 (Stem Cell Technologies 85850) for at least one passage before differentiation into CNCCs or mTeSR Plus (Stem Cell Technologies 100–0276) for gene editing, single-cell cloning, expansion and maintenance. hESCs were grown on Matrigel growth factor reduced basement membrane matrix (Corning 354230) at 37°C. hESCs were fed every day for mTeSR1 or every 2days for mTeSR Plus and passaged every 5–6 days using ReLeSR (Stem Cell Technologies 05872).

HEK293FT cells were obtained from Invitrogen (R70007) and cultured in complete medium (DMEM-HG (GE Healthcare Life Science SH30243.01), 10% FBS, 1× Non-essential amino acids (Gibco 1114-0050), 1× GlutaMAX (Gibco 4109-0036), 1× antibiotic–antimycotic (Gibco 1524-0062)). Cells were fed every other day and passaged every 2–3 days using trypsin–EDTA (Gibco 25200072).

AAV production and titration for CRISPR–Cas9 genome editing

Left and right homology arms (about 1kilobase (kb)) surrounding the SOX9 stop codon, flanking the linker–FKBP12-FV36–linker–mNeonGreen–linker–V5–stop tag, were cloned into an AAV backbone (pAAV-GFP, Addgene plasmid no. 32395). This vector plasmid, along with the AAV6 packaging plasmid pDGM6 (Addgene, plasmid no. 110660), was transfected into 70–80% confluent, early-passage HEK293FT cells seeded 24h before transfection at about 8–9 million cells per 15-cm plate, and changed with fresh medium 2–6h before transfection. For each 15-cm plate (2 per individual AAV6 preparation), the transfection mix was: 22μg pDGM6, 6μg vector plasmid, 120g polyethyenimine (Sigma-Aldrich 408719), and Opti-MEM (Gibco 3198-5070) to 1ml. At 24-h post-transfection, cells were changed into slow-growth medium (same as complete medium but with 1% FBS). Cells were collected 48h after changing to slow-growth medium with the AAVpro Purification Kit Midi (Takara, 6675) as per the manufacturer’s instructions.

Titration of purified AAV6 was carried out by quantitative PCR. Briefly, a previously flash-frozen and thawed 10-μl aliquot of virus was treated with TURBO DNase (Invitrogen, AM2238) as per the manufacturer’s instructions to digest unpackaged DNA. DNase was inactivated by 0.001M EDTA (final concentration) and incubation at 75°C for 10min. Virus DNA was released by proteinase K treatment (1M NaCl, 1% w/v N-lauroylsarcosine, 100μgml⁻¹ proteinase K (Invitrogen, 25530049)) at 50°C for 2h to overnight. Samples were then boiled for 10min, and diluted twice in H₂O such that the final dilution was 1:200,000. DNA standards comprising 10¹⁰–10³ molecules were prepared using AAV6 backbone plasmids containing inverted terminal repeats. Quantitative PCR was carried out on standards and test samples using the Lightcycler 480 Probes Master kit (Roche, 04707494001) with inverted terminal repeat probe and primer sequences indicated in Supplementary Table 5.

Generation of CRISPR–Cas9 and AAV genome-edited cell lines

hESCs were pre-treated with 10μM RHO/ROCK pathway inhibitor Y-27632 (Stem Cell Technologies, 72304) for 2–24h, collected and brought to single cells with Accutase and vigorous pipetting, and about 800,000 were nucleofected with a Cas9–sgRNA RNP complex using the Lonza 4D Nucleofection system. RNP consisted of 17μg Sp-Cas9 HiFi (IDT) and 300pmol sgRNA duplex (sequence in Supplementary Table 5). Cells were plated on Matrigel-coated plates with mTeSR Plus with 10μM Y-27632 and AAV at desired multiplicity of infection (typically about 25,000). Cells were changed into mTeSR Plus with 10μM Y-27632 but no AAV 4–24h after initial plating, and an additional equal volume of mTeSR Plus with no Y-27632 was added 2days later. Subsequent feedings were carried out with no Y-27632 until cells approached confluency, at which point cells were again collected and dissociated to single cells with Accutase (after 10μM Y-27632 pre-treatment) and 500–1000 cells were plated per well of a 6-well plate. Cells were then expanded until colonies were of sufficient size to pick, before which cells were again pre-treated with 10μM Y-27632 for 2–24h. Colonies were picked into 24- or 48-well plates without Y-27632 and allowed to expand for about 5days and passaged 1:2 using ReLeSR, with one half plated on another 24- or 48-well plate and the other half used for lysis with QuickExtract (Lucigen, QE9050). Genotyping PCR was carried out with one primer outside the homology arms and one primer inside the opposite homology arm (see Supplementary Table 5 for sequence). Clones containing the desired knock-in were expanded and used for genomic DNA extraction with the Quick-DNA miniprep kit (Zymo D3024), followed by the same genotyping PCR and Sanger sequencing to confirm knock-in.

Differentiation of hESCs to CNCCs and chondrocytes

hESCs were grown for 5–6 days until large colonies formed, and then they were disaggregated using collagenase IV and gentle pipetting. Clumps of about 200 hESCs were washed in PBS and transferred to a 10cm Petri dish in neural crest differentiation medium (1:1 ratio of DMEM-F12 and Neurobasal, 0.5× Gem21 NeuroPlex supplement with vitamin A (Gemini, 400-160), 0.5× N2 NeuroPlex supplement (Gemini, 400-163), 1× antibiotic–antimycotic, 0.5× Glutamax, 20ngml⁻¹ bFGF (PeproTech, 100-18B), 20ngml⁻¹ EGF (PeproTech, AF-100-15) and 5μgml⁻¹ bovine insulin (Gemini Bio-Products, 700-112P)). After 7–8 days, neural crest emerged from neural spheres attached to the Petri dish, and after 11 days, neural crest cells were passaged onto fibronectin-coated 6-well plates (about 1million cells per well) using Accutase (Sigma-Aldrich A6964) and fed with neural crest maintenance medium (1:1 ratio of DMEM-F12 and neurobasal, 0.5× Gem21 NeuroPlex supplement with vitamin A (Gemini, 400-160), 0.5× N2 NeuroPlex supplement (Gemini, 400-163), 1× antibiotic–antimycotic, 0.5× Glutamax, 20ngml⁻¹ bFGF, 20ngml⁻¹ bFGF EGF and 1mgml⁻¹ BSA (Gemini)). After 2–3 days, neural crest cells were plated at about 1million cells per well of a 6-well plate, and the following day cells were fed with neural crest long-term medium (neural crest maintenance medium+50pgml⁻¹ BMP2 (PeproTech, 120-02)+3μM CHIR-99021 (Selleck Chemicals, S2924; BCh medium)). After transition to BCh medium, CNCCs at subsequent passages were plated at about 800,000 cells per well of a 6-well plate. CNCCs were then passaged twice to passage 4, at which depletion experiments were carried out, or cells were further differentiated to chondrocytes. For depletion experiments, dTAG^V-1 (Tocris, 6914) at a range of concentrations was added to BCh medium, with an equivalent amount of dimethylsulfoxide (DMSO) as vehicle control.

To differentiate CNCCs to chondrocytes, passage 3 CNCCs were passaged to passage 4, seeded at about 250,000 CNCCs per well of a 12-well plate, and grown for 3days in BCh medium. Then, CNCCs were transitioned to chondrocyte medium without TGFb3 (ChM: DMEM-HG, 5% FBS, 1× ITS premix, 1mM sodium pyruvate, 50μgml⁻¹ ascorbic acid, 0.1μM dexamethasone and 1× antibiotic–antimycotic), with or without dTAG^V-1. The following day, cells were fed with chondrocyte medium with TGFb3 (ChMT: ChM+10ngml⁻¹ TGFb3), with or without dTAG^V-1. Cells were fed every subsequent 3days with ChMT. Cells were collected at day 10 and/or day 21 of the differentiation.

sGAG quantification

Total sGAG levels per well of chondrocytes independently differentiated from CNCCs for 21 days, representing mature cartilage formation, were quantified using the Blyscan glycosaminoglycan assay (Biocolor). Briefly, collagen in the extracellular matrix was digested by washing cells with PBS and then adding 1ml of Papain digestion buffer per well of a 12-well plate. Cells were incubated at 65°C for 3h with gentle agitation every 30min, then 0.5ml additional digestion buffer was added and lysate was moved to Eppendorf tubes and incubated at 65°C overnight. Quantification of sGAG content from about 10μl of the lysates was carried out as per the manufacturer’s instructions, and the volume of each lysate was measured separately and used to infer the total sGAG content of the entire well.

Flow cytometry

CNCCs were collected for flow cytometry using Accutase and quenching with FACS buffer (5% FBS in PBS). Chondrocytes were incubated in digestion medium (DMEM-KO, 1mgml⁻¹ Pronase (Roche, 11459643001), 1mgml⁻¹ collagenase B (Roche, 11088815001), 4Uml⁻¹ hyalauronidase (Sigma, H3506-500mg)) for about 1h with gentle agitation every 15min. Digested cells were then washed twice in PBS. Flow cytometry was used to measure mNeonGreen fluorescence after excluding doublets and debris based on forward and side scatter (Beckman Coulter Cytoflex). Fluorescence values were summarized per biological replicate using geometric means. The relative SOX9 level as percentage of the SOX9-tagged, unperturbed (treated with DMSO) sample was calculated by first subtracting the geometric mean fluorescence of the untagged (WT) sample from both the unperturbed and dTAG^V-1-treated sample, and then dividing the dTAG^V-1-treated sample fluorescence by the unperturbed sample fluorescence.

Protein collection and western blotting

Cells were washed with PBS and scraped into RIPA buffer (50mM Tris, 150mM NaCl, 1% NP-40, 0.1% Nadeoxycholate, 0.1% SDS in H₂O with 1× protease inhibitor cocktail (Sigma-Aldrich 4693132001)), incubated on ice for 10min, and sonicated to disrupt pelleted DNA using Bioruptor Plus (Diagenode). Sonicated lysates were incubated on ice for 10min, and centrifuged at 16,000g for 10min at 4°C to pellet debris. Supernatants were normalized to the same protein content using the Pierce BCA Protein Assay kit (ThermoFisher, 23225), mixed with 4× SDS sample loading buffer (Invitrogen NP0007) and 0.1M dithiothreitol (DTT), and boiled for 7min. Samples were separated on Tris-glycine polyacrylamide gel electrophoresis (PAGE) gels in 1× Tris-glycine buffer with 0.1% SDS, transferred in 1× Tris-glycine buffer with 20% methanol, blocked in 5% milk+1% BSA in PBST, immunoblotted with either SOX9 antibody (1:1,000, Sigma-Aldrich AB5535) or β-actin antibody (1:20,000, Abcam ab49900) overnight at 4°C, probed with the appropriate secondary, developed using Pierce ECL Western Blotting Substrate (ThermoFisher, 32106), and imaged using an Amersham ImageQuant 800 system (Cytiva).

RNA isolation and preparation of RNA-seq libraries

Total RNA was extracted from CNCCs using Trizol reagent (Invitrogen) followed by Quick-RNA Miniprep kit (Zymo) with on-column DNase I digestion. Unstranded mRNA libraries were prepared with the NEBNext Ultra II RNA Library Prep Kit for Illumina (NEB no. E7770S/L).

Metabolic RNA labeling and preparation of SLAM-seq libraries

4-Thiouridine was incorporated into nascent transcripts by incubating CNCCs with BCh medium containing 100μM 4-thiouridine, as well as DMSO or 500nM dTAG^V-1 depending on experimental condition, for 2h. Plates were covered in foil and handling was carried out in a hood with no light where possible. For 3- and 24-h depletion experiments, labeling was started at 1 and 22h after dTAG^V-1 addition, respectively.

Total RNA was extracted using Trizol reagent, phenol–chloroform extraction was carried out, and the aqueous phase was used as input to the Quick-RNA Miniprep kit. During RNA extraction with Quick-RNA Miniprep kit, 0.1mM DTT was added to the RNA wash and RNA pre-wash buffers, but the on-column DNase I step was skipped. RNA was eluted in H₂O with 1mM DTT, quantified with Qubit RNA Broad Range assay (ThermoFisher, Q10211), and >2μg total RNA was used as input to the alkylation reaction. Alkylation was carried out in dark tubes after which light exposure was allowed, and after quenching RNA was purified and subjected to on-column DNase I digestion using the RNA Clean & Concentrator-5 kit (Zymo, R1013).

A 500ng quantity of alkylated RNA was used as input to QuantSeq 3′ mRNA-Seq Library Prep Kit FWD with unique dual index add-on (Lexogen, 113.96), with 15 cycles of PCR amplification. Library size distributions were confirmed by separation on a PAGE gel and staining with SYBRGold and pooled on the basis of quantifications from Qubit dsDNA High Sensitivity Kit (ThermoFisher Q32854). Pooled libraries were sequenced using Novaseq 6000 platform (2× 150bp).

ATAC-seq collection and library preparation

CNCCs were incubated with BCh medium containing 200Uml DNase I (Worthington, LS002007) for 30min and collected using Accutase. Viable cells were counted using a Countess Automated Cell Counter (Invitrogen), and 50,000 viable cells were pelleted at 500RCF for 5min at 4°C and resuspended in ATAC-resuspension buffer (10mM Tris-HCl pH 7.4, 10mM NaCl, 3mM MgCl₂ in sterile water) containing 0.1% NP-40, 0.1% Tween20 and 0.01% digitonin and incubated on ice for 3min. Following wash-out with cold ATAC-resuspension buffer containing 0.1% Tween20, cells were pelleted and resuspended in 50μl transposition mix (25μl 2× TD buffer, 2.5μl transposase (100nM final), 16.5μl PBS, 0.5μl 1% digitonin, 0.5μl 10% Tween20, 5μl H₂O) and incubated for 30min at 37°C with shaking. The reaction was purified using the Zymo DNA Clean & Concentrator kit and PCR-amplified with NEBNext High-Fidelity 2× PCR Master Mix (NEB, M0541L) and primers as defined in ref.⁶². Libraries were purified by two rounds of double-sided size selection with AMPure XP beads (Beckman Coulter, A63881), with the initial round of 0.5× sample volume of beads followed by a second round with 1.3× initial volume of beads. Library size distributions were confirmed by separation on a PAGE gel and staining with SYBRGold and pooled on the basis of quantifications from Qubit dsDNA High Sensitivity Kit. Pooled libraries were sequenced using the Novaseq 6000 platform (2×150bp).

ChIP and library preparation

One fully confluent 10-cm plate of cells was crosslinked per ChIP experiment in 10ml PBS with 1% methanol-free formaldehyde for 10min and quenched with a final concentration of 0.125M glycine for 10min with nutation. Crosslinked cells were scraped into tubes with 0.001% Triton X in PBS, washed with PBS without Triton, pelleted by centrifugation, flash-frozen in liquid nitrogen and stored at −80°C. Samples were defrosted on ice and resuspended in 5ml LB1 (50mM HEPES-KOH pH7.5, 140mM NaCl, 1mM EDTA, 10% glycerol, 0.5% NP-40, 0.25% Triton X-100, with 1× cOmplete Protease Inhibitor Cocktail and phenylmethylsulfonyl fluoride) and rotated vertically for 10min at 4°C. Samples were centrifuged for 5min at 1,350g at 4°C, and resuspended in 5ml LB2 (10mM Tris, 200mM NaCl, 1mM EDTA, 0.5mM EGTA, with 1× cOmplete Protease Inhibitor Cocktail and optionally 1mM phenylmethylsulfonyl fluoride) and rotated vertically for 10min at 4°C. Samples were centrifuged for 5min at 1,350g at 4°C, resuspended in 300μl LB3 per sonicated sample, and incubated for 10min on ice. Samples were sonicated in 1.5ml Bioruptor Plus TPX microtubes (Diagenode, c30010010-50) on Bioruptor Plus for 10 cycles of 30s on–30s off. Every 5 cycles, samples were lightly vortex and briefly centrifuged. Samples were diluted in additional LB3 to 1ml, pelleted at 16,000 RCF for 10min, and the supernatant was removed. Triton X-100 was added to 1%.

To check DNA size distribution and quantity, a 10-μl aliquot of sonicated chromatin from each sample was diluted to 100μl in elution buffer (50mM Tris, 10mM EDTA, 1% SDS) with 0.0125M NaCl and 0.2mgml⁻¹ RNase A and incubated at 65°C for 1h, followed by addition of proteinase K to 0.2mgml⁻¹ and an additional 1h of 65°C incubation. DNA was purified using Zymo DNA Clean & Concentrator Kit with ChIP DNA Binding Buffer (Zymo, D5201-1-50) and size distribution and quantity was assessed by separation on a 1% agarose gel and Qubit HS DNA kit, respectively. Qubit measurements were used to normalize samples to the same DNA concentration.

Following normalization, the chromatin was divided for input (2%) and ChIP samples. A minimum of 25μg DNA was used for histone ChIP analyses, and 50μg for V5 ChIP analyses. A 5μg quantity of anti-H3K27ac (Active Motif, 39133) antibody (1:200 dilution) or a 10μg quantity of anti-V5 (Abcam, ab9116 or ab15828) or TWIST1 (Abcam, ab50887) antibody (1:100 dilution) was added per ChIP sample, and incubated overnight at 4°C. Protein G Dynabeads (ThermoFisher) were first blocked with Block solution (0.5% BSA (w/v) in 1× PBS) and then added to cleared chromatin to bind antibody-bound chromatin for a 4–6h incubation. Chromatin-bound Dynabeads were washed at least 6 times with chilled RIPA wash buffer (50mM HEPES-KOH pH7.5, 500mM LiCl, 1mM EDTA, 1% NP-40, 0.7% Nadeoxycholate), followed by a wash with chilled TE+50mM NaCl. Chromatin was eluted for 30min in elution buffer (50mM Tris, 10mM EDTA, 1% SDS) at 65°C with frequent vortexing. The ChIP and input samples were then incubated at 65°C overnight to reverse crosslinks (12–16h). Samples were diluted and sequentially digested with RNase A (0.2mgml⁻¹) for 2h at 37°C followed by proteinase K (0.2mgml⁻¹) for 2h at 55°C for 2–4h to digest protein. ChIP and input samples were purified by Zymo DNA Clean & Concentrator Kit with ChIP DNA binding buffer.

For library preparation, samples were quantified by Qubit dsDNA HS assay kit, and 10–50ng of ChIP DNA was used for library preparation with end repair, A-tailing and adaptor ligation (NEB). Following USER enzyme treatment, libraries were cleaned up with one round of single-side AMPure XP bead clean-up, and then amplified to add indices using NEBNext Ultra II Q5 Master Mix and NEBNext Multiplex Oligos for Illumina kit (NEB, E7335S) with 4–10 cycles (as determined by input amounts from NEB protocol). ChIP libraries were purified by two rounds of double-sided AMPure XP bead clean-up (0.5× then 0.4× initial sample volume of beads added) to remove large fragments and deplete adaptors. Library concentration and quality within ChIP or input groups was assessed by Qubit dsDNA HS assay kit and separation on a PAGE gel, and used to pool within ChIP or input groups. KAPA quantitative PCR was used to pool across ChIP or input groups. Pooled libraries were sequenced using the Novaseq 6000 platform (2×150bp).

Sequencing data pre-processing

Quantification and statistical analysis

PRS endophenotype definition and GWAS

Phenotyping

GWAS

The scores on the 30 PRS-driven univariate traits, for which a significant difference was observed, were combined into a single phenotype matrix ([N×M] with N=8,246 controls and M=30 facial segments). This matrix was tested for genotype–phenotype associations in a multivariate meta-analysis framework using canonical correlation analysis (canoncorr in Matlab 2017b). However, instead of carrying out a separate GWAS per facial segment, information across multiple segments is now combined into a single multivariate GWAS. As canonical correlation analysis does not accommodate adjustments for covariates, we removed the effect of relevant covariates (sex, age, age-squared, height, weight, facial size, four genomic ancestry axes, camera system), on both the independent (SNP) and the dependent (facial shape) variables using partial least-squares regression (plsregress from Matlab 2017b) before GWAS.

The US and UK subsamples served both as identification and replication sets in a two-stage design, after which the P values were meta-analyzed using Stouffer’s method^67,68. Per SNP, the lowest P value was selected (meta_US versus meta_UK) and compared against the genome-wide Bonferroni threshold (5×10⁻⁸). We observed 1,767 SNPs at the level of genome-wide significance, which were clumped into 22 independent loci as follows. Starting from the lead SNP (lowest P value), SNPs within 10kb or within 1Mb but with r²>0.01 were clumped into the same locus represented by the lead SNP. Next, considering only the lead SNPs, signals within 10Mb and with an r²>0.01 were merged. Third, any locus with a singleton lead SNP was removed.

Post-GWAS analyses

To define facial shape, GWAS SNPs that affect facial shape in either a PRS-like or non-PRS-like manner, we obtained the combined list of SNPs affecting normal-range variation (and not orofacial clefting) in facial shape from ref.²⁴. We tested each of these SNPs directly for association with the PRS endophenotype from the above-described GWAS. SNPs with a Bonferroni-corrected Pvalue<0.01 (corresponding to an uncorrected P-value cutoff of 7×10⁻⁵) in either the UK or UK cohort from the PRS GWAS were considered significant. We then considered all SNPs in linkage disequilibrium (r²>0.5, identified with SNiPA⁶⁹ ‘Proxy Search’ tools using 1000 Genomes Phase 3 v5 European reference panel) with either set of facial GWAS SNPs (no PRS endophenotype association or significantly associated). SOX9-dependent REs containing these SNPs were assigned as ‘PRS-like’ or affecting other aspects of facial shape according to the type of linked SNP they contained.

We thank members of the laboratories of J.W. and J.K.P., as well as James Ferrell, for helpful comments. This research was supported by National Institutes of Health (NIH) grants R01HG008140 and R01HG011432 (J.K.P.), and R35GM131757-1 (J.W.), and financial support from the Howard Hughes Medical Institute and the Nomis Foundation. S.N. was supported by a Helen Hay Whitney Fellowship. S.K. was supported by a Damon Runyon Fellowship (DRG-2420-21). J.W. was supported by a Lorry Lokey endowed professorship and a Stinehart Reed award. The KU Leuven research team (P.C., H.H., H.S.M.) was supported by the NIH (1-R01-DE027023, 2-R01-DE027023), The Research Fund KU Leuven (BOF-C1, C14/15/081 and C14/20/081) and The Research Program of the Research Foundation - Flanders (Belgium) (FWO, G078518N). R.A.S., B.H. and O.D.K. were supported by NIH U01DE024440. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government – department EWI. Some of the computing for this project was carried out on the Sherlock cluster. We thank Stanford University and the Stanford Research Computing Center for providing computational resources and support that contributed to these research results. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.