Affinity-optimizing enhancer variants disrupt development

Abstract

Enhancers control the location and timing of gene expression and contain the majority of variants associated with disease^1-3. The ZRS is arguably the most well-studied vertebrate enhancer and mediates the expression of Shh in the developing limb⁴. Thirty-one human single-nucleotide variants (SNVs) within the ZRS are associated with polydactyly^4-6. However, how this enhancer encodes tissue-specific activity, and the mechanisms by which SNVs alter the number of digits, are poorly understood. Here we show that the ETS sites within the ZRS are low affinity, and identify a functional ETS site, ETS-A, with extremely low affinity. Two human SNVs and a synthetic variant optimize the binding affinity of ETS-A subtly from 15% to around 25% relative to the strongest ETS binding sequence, and cause polydactyly with the same penetrance and severity. A greater increase in affinity results in phenotypes that are more penetrant and more severe. Affinity-optimizing SNVs in other ETS sites in the ZRS, as well as in ETS, interferon regulatory factor (IRF), HOX and activator protein 1 (AP-1) sites within a wide variety of enhancers, cause gain-of-function gene expression. The prevalence of binding sites with suboptimal affinity in enhancers creates a vulnerability in genomes whereby SNVs that optimize affinity, even slightly, can be pathogenic. Searching for affinity-optimizing SNVs in genomes could provide a mechanistic approach to identify causal variants that underlie enhanceropathies.

Fig. 1. An ETS-A site in the ZRS enhancer contains two human variants that are associated with polydactyly, both of which subtly increase ETS binding affinity.

a, The human ZRS contains five known and functionally validated ETS sites, ETS1–ETS5, all of which have suboptimal affinity. We identify a new site, ETS-A, which has a relative affinity of 0.15. Six HOX sites (yellow) and one HAND2 site (pink) have also been previously identified. Thirty-one SNVs associated with polydactyly are found in humans (black bars), and SNVs are also found in other species such as cats, mice and chicks (green bars). b, Two human SNVs associated with polydactyly, denoted French 2 and Indian 2, occur within the ETS-A site. Both SNVs lead to a subtle increase in relative affinity of ETS-A to 0.24 and 0.26 respectively. c–f, The ETS-A sequence in a reference (Ref) mouse (c) drives the expression of Shh (d) and Ptch1 (e) restricted to the posterior domain of the developing limb bud in E11.75 and E12.0 embryos, respectively, as shown by in situ hybridization. f. Skeletal staining shows a WT mouse hindlimb with normal digit morphology. g–j, The French 2 SNV (g) drives the ectopic expression of Shh (h) and Ptch1 (i) in the anterior limb bud of homozygous embryos (arrow) in addition to the normal domain of posterior expression. j, A mouse hindlimb homozygous for French 2 has an extra triphalangeal thumb. k–n, The Indian 2 SNV (k) drives the ectopic expression of Shh (l) and Ptch1 (m) in homozygous embryos (arrow). n, A hindlimb from an Indian 2 homozygous mouse has an extra triphalangeal thumb. We did not calculate n for Shh because the expression is highly dynamic and thus hard to accurately capture; instead, we calculate the n of Ptch1 as a readout of Shh.

Fig. 2. Synthetic changes to the ETS-A site that create a 0.25-affinity site and a LOF site cause predicted phenotypes.

a–d, The Syn 0.25 ETS-A site (a), which has an affinity of 0.25, drives the ectopic expression of Shh (b) and Ptch1 (c) in the anterior limb bud, in addition to the normal domain of expression in the posterior limb bud. d, Skeletal staining of a homozygous mouse hindlimb shows an extra triphalangeal thumb. e–h, The LOF ETS-A site (e) drives the normal expression of Shh (f) and Ptch1 (g) in the posterior limb bud. h, Homozygous mice have normal digit morphology. We did not calculate n for Shh because the expression is highly dynamic and thus hard to accurately capture; instead, we calculate the n of Ptch1 as a readout of Shh.

Fig. 3. All mice with the approximately 0.25-affinity ETS-A are indistinguishable in terms of the penetrance, laterality and severity of polydactyly.

a, Penetrance and laterality of phenotypes seen in French 2, Indian 2 and Syn 0.25 heterozygous (Het) and homozygous (Homo) mice. There is no significant difference in penetrance and laterality between heterozygotes and between homozygotes of French 2, Indian 2 and Syn 0.25 mice (Fisher’s exact test). Polydactyly occurs more frequently on the right hindlimb in both heterozygotes and homozygotes. In Syn 0.25/LOF mice, polydactyly is less penetrant than in heterozygous Syn 0.25 mice (Syn 0.25/WT). NS, not significant. b, Polydactyly phenotypes seen in French 2, Indian 2 and Syn 0.25 heterozygotes (dashed lines) and homozygotes (solid lines), as well as in Syn 0.25/LOF mice. There is no significant difference in the phenotypes seen between the French 2, Indian 2 and Syn 0.25 lines (Fisher’s exact test). TT, triphalangeal toe. Biorender.com was used to generate images in this figure. Source data

Fig. 4. A greater increase in affinity at the ETS-A site causes more severe and penetrant polydactyly and long-bone defects.

a, The ETS-A reference (Ref) sequence. b, In situ hybridization of Shh in WT hindlimb at E11.75 shows a domain of expression in the posterior of the limb bud. c, Skeletal staining of a WT mouse hindlimb showing five digits. d, WT forelimb showing five digits e, WT hindlimb focused on tibia and fibula. f, The Syn 0.52 ETS-A sequence. g, In situ hybridization of Shh in a Syn 0.52 homozygous mouse hindlimb bud at E11.75. Ectopic expression in the anterior of the limb bud is indicated by an arrow. h, Skeletal staining of a homozygous Syn 0.52 ETS-A mouse hindlimb showing six digits. The toe and the extra digit are both triphalangeal. i, Homozygous Syn 0.52 mouse forelimb showing six digits. Both the thumb and the extra digit are triphalangeal. j, Hindlimb focused on the tibia and fibula in a Syn 0.52 ETS-A homozygous mouse. The tibia is severely shortened. The hindlimbs in e,j were imaged with the same magnification. See Extended Data Fig. 6 for more details.

Fig. 5. Affinity-optimizing SNVs drive GOF expression in the ZRS and IFNβ enhanceosome.

a, Schematic of the human ZRS showing 19 putative ETS sites and 17 SNVs that increase the affinity of ETS sites by 1.6-fold or more. ETS sites identified from previous studies (ETS0–ETS5) are labelled^,. We annotated changes in expression caused by these SNVs on the basis of results from a saturation mutagenesis MPRA and our study. b,c, Box plots showing all SNVs tested within MPRA mutagenesis experiments, their significance and their effects on expression. The bounds of the box plots define the 25th, 50th and 75th percentiles, and whiskers are 1.5× the interquartile range. The y axis is the log of the adjusted P value (P_adj) with the direction of expression change (positive values indicate an increase in expression and negative values indicate a decrease in expression). Dashed horizontal lines indicate significance thresholds at P = .05. Each dot represents a tested SNV. SNVs that have a significant increase in expression are shown in green, those with no significant change in expression are grey and those with a significant decrease in expression are red. For each plot, we show SNVs that increase the affinity of the transcription factor, SNVs within the TFBS that do not alter affinity and all other SNVs. b, Analysis of the ETS TFBS in a saturation mutagenesis MPRA on a 485-bp region of the ZRS enhancer. c, Analysis of the IRF TFBS in a saturation mutagenesis MPRA on the IFNβ enhanceosome. We used one-tailed Mann–Whitney U tests to determine any significant enrichment for GOF enhancer activity. Source data

Fig. 6. Affinity-optimizing SNVs drive GOF expression in a wide variety of disease-associated enhancers.

a–d, Analysis of MPRAs for a variety of enhancers in different cell types. Box plots showing all SNVs tested within MPRA mutagenesis experiments and their significance and effects on expression. The bounds of the box plots define the 25th, 50th and 75th percentiles, and whiskers are 1.5× the interquartile range. One-tailed Mann–Whitney U test. a,b, Analysis of saturation mutagenesis MPRA assays of 11 disease-associated enhancers comparing all SNVs, SNVs within TFBSs that do not alter affinity and SNVs that increase the affinity of TFBSs for ETS (a) and AP-1 (b). c,d, Analysis of MPRA comparing the effect of SNVs within lymphoblastoid regulatory elements for the ETS TFBS (c) and the AP-1 TFBS (d). e, Filtering for affinity-optimizing (aff-opt) SNVs significantly increases our ability to predict causal GOF enhancer variants. The bar graph shows the percentage of all SNVs that lead to GOF expression relative to the percentage of affinity-optimizing SNVs that lead to GOF expression. Green bars indicate SNVs that cause GOF expression within analysed MPRA datasets; yellow bar indicates SNVs that cause GOF expression in our the current study—namely, French 2 and Indian 2. Fisher’s exact test was used to determine any significant enrichment for GOF expression in the all and aff-opt categories: **P < 0.001, *P < 0.01. Affinity-optimizing SNVs are those that lead to a fold change of at least 1.6 for ETS because this is the fold change for French 2, Indian 2 and Syn 0.25. Affinity-optimizing SNVs for AP-1 and IRF are those that cause a fold change of at least 1.5. Source data

Extended Data Fig. 1. Class I ETS family members have conserved DBDs.

DBDs of ETS-1 transcription factors. Numbers show the location of contacts between the DNA-binding site and the protein in human and mouse^,,. The DBDs are highly conserved across human, mouse and flies across other class I ETS family members^,,.

Extended Data Fig. 2. PBM binding affinities correlate with the in vivo ETS-1 ChIP signal in various cell types.

a, ETS-1 ChIP–seq signals from mouse primary B cells (y axis) and PBM mouse ETS-1 affinities (x axis) show a strong Spearman’s rank correlation of 0.88 (ref. ). To show that the correlation is not heavily dependent on values from affinity bins 0.1 and 1.0, we show the Spearman’s rank correlation between cumulative affinity bins below the graph. For example, the Spearman’s correlation using bins 0.1, 0.2, 0.3, 0.4 is 0.703. b, The ETS-1 ChIP signal is low for all random non-ETS k-mers c, Correlation of ChIP signal and PBM affinity for other ETS-1 ChIP assays in mouse primary B cells, human natural killer cells and human THP-6 cell lines^–.

Extended Data Fig. 3. Conservation of ZRS ETS sites between humans and mice.

a, The human and mouse ZRS sequence is highly conserved. The human ZRS has 19 putative ETS-binding sites, all with affinities equal to or lower than 0.52. The mouse ZRS has 6 functionally validated and 12 putative ETS sites; 15 of these are conserved in location and affinity with the human ZRS. b, The ETS-A site and surrounding sequence show perfect conservation between mouse and human, as indicated by asterisks. Blue box highlights the ETS-A sequence within the human and mouse ZRS.

Extended Data Fig. 4. EMSA shows the binding of human and mouse ETS-1 to the ETS-A site and ETS-A variants.

a, WT ETS-A probe sequence can bind to the ETS-1 DBD to generate a band, the intensity of which decreases as the biotin-labelled probe is outcompeted by a non-biotin competitor probe. A single-bp change in the ETS-A site (LOF) leads to a loss of binding signal. b, EMSAs for WT, French 2, Syn 0.25, Indian 2 and Syn 0.52 sequences. The total amount of bound probe relative to the unbound probe is not statistically different between French 2, Indian 2 and Syn 0.25 sequences, suggesting that all three sequences have the same affinity. P = 0.18, one-way ANOVA. The Syn 0.52 sequence binds more strongly than the 0.25 or 0.15 sequences do. EMSAs were performed independently twice and both replicates show similar results. For gel source data, see Supplementary Fig. 1.

Extended Data Fig. 5. Ptch1 in situ hybridization in the hindlimb bud and forelimb bud of transgenic mice.

Embryos were collected at around the 53-somite stage for Ptch1 in situ hybridization. a, Ptch1 expression is restricted to the posterior domain in WT hindlimb buds. Ectopic Ptch1 expression can be seen in the anterior domain of homozygous French 2 (b), Indian 2 (c), Syn 0.25 (d) hindlimb buds. LOF homozygotes (e) and embryos with Syn0.25/LOF alleles (f) do not have ectopic Ptch1 expression as predicted. g, Syn 0.52 hindlimb buds have larger domain of ectopic Ptch1 expression than embryos with approximately 0.25 affinity ETS-A sites. h, In WT forelimb buds, Shh and Ptch1 expression is restricted to the posterior domain. i, Syn 0.52 homozygotes display ectopic Shh expression and Ptch1 expression in the forelimb buds. The Shh in situ hybridizations shown were performed on embryos at around the 48-somite stage. All limb bud images were acquired and cropped using the same settings.

Extended Data Fig. 6. Syn 0.52 mice show highly penetrant polydactyly in both forelimb and hindlimb, as well as tibial hemimelia.

a, In the hindlimb, heterozygotes and homozygotes have 100% penetrance in polydactyly. In the forelimb, homozygotes have 100% penetrance whereas heterozygotes have 93.6% penetrance. b, Tibial hemimelia is observed in 95.3% of homozygotes but no heterozygotes. c, Digit phenotypes on the hindlimbs of all mice studied. d, Digit phenotypes on the forelimbs of Syn 0.52 mice. WT, LOF, French 2, Indian 2 and Syn 0.25 have no forelimb phenotypes. TT denotes triphalangeal toe or thumb. Source data

Extended Data Fig. 7. Affinity-optimizing SNVs are significantly associated with GOF enhancer activity.

Enhancers containing SNVs that increase ETS affinity are significantly enriched in GOF expression relative to all other SNVs within a 2% sequence mutagenesis assay of the ZRS enhancer. The bounds of the box plots define the 25th, 50th and 75th percentiles, and whiskers are 1.5 × the interquartile range. One-tailed Mann–Whitney U test. Source data

Extended Data Fig. 8. EMSA shows stronger binding of the human HOXA13 and HOXD13 DBDs to the Dutch 2 variant relative to the WT sequence.

a, WT probe sequence can bind to the human HOXA13 DBD to generate a band, the intensity of which decreases in the presence of the non-biotin competitor probe. The Dutch 2 variant binds to HOXA13 more strongly than does the WT allele. A single-bp change that ablates HOXA13 binding shows no binding. b, EMSA with the same WT and Dutch 2 probe sequences as in a, performed with the human HOXD13 DBD. The Dutch 2 variant binds to HOXD13 more strongly than does the WT allele. EMSAs were performed independently twice and both replicates show similar results. For gel source data, see Supplementary Fig. 1.

Extended Data Fig. 9. Experimental details for MPRA performed with 11 disease-associated enhancers.

Eleven disease-associated enhancers tested in saturation mutagenesis MPRAs; table modified from a previous study. Two different MYC enhancers were assayed. UC88 is an ultraconserved enhancer. The MPRA for each enhancer is performed within cell lines relevant to the disease studied, as detailed.

Extended Data Fig. 10. Affinity-optimizing eQTL variants are enriched in GOF target gene expression.

a, Seven ETS affinity-optimizing variants within lymphoblastoid regulatory element MPRA drive significant GOF reporter expression. Five of these are associated with significant eQTL differential expression and all of these SNVs drive increased target gene expression. Dotted lines indicate thresholds for significance at P < 0.05 (significant GOF in green, and significant LOF in red). b, eQTL analysis for all MPRA variants tested within the lymphoblastoid regulatory elements (regardless of variant effect within MPRA assay)^,. ETS affinity-optimizing SNVs are enriched in GOF target gene expression. The bounds of the box plots define the 25th, 50th and 75th percentiles, and whiskers are 1.5 × the interquartile range. One-tailed Mann–Whitney U test. c–e, Genome-wide analyses for all significant eQTLs (P < 0.01) within the lymphoblastoid cell line show that affinity-optimizing variants are enriched for GOF expression. Line indicates mean values. One-tailed Mann–Whitney U Test. c, ETS SNVs with ≥1.59 affinity fold change. d, ETS SNVs with ≥3.46 affinity fold change; this is the fold change of 0.52 mice relative to WT 0.15 ETS-A. e, AP-1 SNVs ≥1.5 affinity fold change. Source data