Regulatory elements in SEM1-DLX5-DLX6 (7q21.3) locus contribute to genetic control of coronal nonsyndromic craniosynostosis and bone density-related traits

Abstract

Purpose: The etiopathogenesis of coronal nonsyndromic craniosynostosis (cNCS), a congenital condition defined by premature fusion of 1 or both coronal sutures, remains largely unknown.

Methods: We conducted the largest genome-wide association study of cNCS followed by replication, fine mapping, and functional validation of the most significant region using zebrafish animal model.

Results: Genome-wide association study identified 6 independent genome-wide-significant risk alleles, 4 on chromosome 7q21.3 SEM1-DLX5-DLX6 locus, and their combination conferred over 7-fold increased risk of cNCS. The top variants were replicated in an independent cohort and showed pleiotropic effects on brain and facial morphology and bone mineral density. Fine mapping of 7q21.3 identified a craniofacial transcriptional enhancer (eDlx36) within the linkage region of the top variant (rs4727341; odds ratio [95% confidence interval], 0.48[0.39-0.59]; P = 1.2E-12) that was located in SEM1 intron and enriched in 4 rare risk variants. In zebrafish, the activity of the transfected human eDlx36 enhancer was observed in the frontonasal prominence and calvaria during skull development and was reduced when the 4 rare risk variants were introduced into the sequence.

Conclusion: Our findings support a polygenic nature of cNCS risk and functional role of craniofacial enhancers in cNCS susceptibility with potential broader implications for bone health.

Keywords: Coronal Nonsyndromic; Craniosynostosis; DLX6 DLX5; GWAS; Regulatory elements; SEM1.

Graphical abstract

Figure 1

Summary of the genome-wide analysis of coronal nonsyndromic craniosynostosis. A. Manhattan plot of the meta-analysis of discovery and replication genome-wide association analyses using common variants (minor allele frequency > 1%). The y-axis shows the −log₁₀ transformed P value of each variant association found using a standard-error-weighted approach and controlling for population stratification, and the x-axis shows the chromosomal position. Variants crossing the genome-wide significance threshold of P < 5E−08 are color coded in red, and those with P < 5E−06 are in green. The top signals are annotated with the closest genes. Inset: quantile-quantile plot showing distribution of expected P values under the null model (red-dotted line) vs observed P values (black dots). B. Regional plot of the 4 top independent genomic association signals from the European meta-analyses. The y-axis shows −log₁₀P values for individual variants annotated with the genes in the selected genomic interval. The top variants are marked as purple diamonds and other variants in pairwise linkage disequilibrium (r²) with the top variant, based on the 1000 Genomes Project Phase 3 European reference samples, are color coded as per the scale in legend.

Figure 2

Multilocus analysis of the top susceptibility loci. A. Polygenic risk score calculated using genome-wide association study summary statistics from our discovery cohort was used to predict the risk for coronal nonsyndromic craniosynostosis in the replication cohort. Best fit model was achieved by PRsice with 106 variants with association P < 5E−05. Inset: difference in means of Polygenic risk score values in the craniosynostosis cases (blue) and the control group (yellow) are shown. B. Knowledge graph connecting the identified genes with shared enriched functional terms from Enrichr. In the network, identified genes are represented as orange ovals, whereas shared enriched annotations from Enrichr are shown as blue rectangles. Known physical interactions between the protein products of the identified genes are depicted by red lines, connections to functional terms are depicted by gray lines, and related terms are connected by blue lines.

Figure 3

Phenome-wide association analysis of rs4727341, the top risk variant. Phenotype-wide association analysis of complex traits associated with rs4727341. Summary statistics from the UK Biobank, FinnGen, and genome-wide association study catalog repositories were downloaded from Open Target (https://genetics.opentargets.org/). Only traits with P value < .005 are shown in the diagram. x axis shows traits and y axis shows the variant’s P value of association to each trait. The circles are color coded by the trait category (see legend) as reported in Open Target website. The red dashed line shows the significance threshold corrected for the number of traits shown. In the figure, heel bone mineral density and other traits appear multiple times since the association was reported in many independent studies/publications as follows: heel bone mineral density (Heel BMD): ¹GCST006979, ²GCST006288, ³NEALE2_3148_raw, ⁴NEALE2_78_raw (t score automated), ⁵NEALE2_4125_raw (t score automated right), ⁶NEALE2_4124_raw (right); cortical surface area: ¹GCST010282_20: pars triangularis, ²GCST010701: MOSTest, ³GCST010697: min P and ⁴GCST90091060; other heel measurements are also shown—heel broadband ultrasound attenuation (heel bua): NEALE2_3144_raw: direct entry and NEALE2_4120_raw: right); and heel quantitative ultrasound index (heel qui): NEALE2_3147_raw: direct entry and NEALE2_4123_raw: right.

Figure 4

Enhancer analysis of the SEM1 locus. A. A regional association plot surrounding the top risk variant, rs4727341 (shown as a purple diamond). The x axis represents a 0.2 Mb region, 100 kb upstream and downstream of the lead variant; the y axis shows −log₁₀P values for individual variant associations from the European meta-analysis annotated with the genes in the selected genomic interval. Pairwise linkage disequilibrium (LD) (r²) with the lead variant color coded based on the 1000 Genomes Project Phase 3 European reference samples. B. Zoom-in of the top signal region to highlight the rs4727341-LD region (r² > 0.9) (dashed black box). The y axis shows −log₁₀P values for individual variant associations from the European meta-analysis as above. Pairwise LD (r²) with the lead variant color coded based on the 1000 Genomes Project Phase 3 European reference samples. C. Zoom-in genomic region annotated with −log₁₀P values for individual 6 kb sliding windows from rare variant TDT aggregate analysis of rare variants in family-based study (purple dots). The dot is represented at the start of each window. D. Predicted craniofacial regulatory elements near rs4727341. Histone modification marks are associated with active craniofacial predicted enhancers (CNCC1 through F2). Highlighted are 3 predicted enhancer candidates (Enhan.), eDlx34, eDlx35, and eDlx36, which were tested in this study (orange bar). Phylop conservation track (Cons.) from University of California Santa Cruz (UCSC) Genome Browser is shown in black color. E. Zebrafish enhancer assays at 3 days after fertilization (dpf). eDlx34 drove green fluorescent protein (GFP) expression in the heart and somitic muscles. eDlx35 drove specific GFP expression in the mandibular arch and branchial arches (basibranchials, hypobranchials, and ceratobranchial 1-5) and notochord. eDlx36 drove specific GFP expression in the premaxillary, maxillary, FNP, and apical region of the skull. F. The effect of 4 rare variants on the in vivo eDlx36 enhancer activity. GFP-positive cells drove by eDlx36 enhancer in the FNP and apical region of the skull at 7 and 14 dpf, whereas the mutated eDlx36 embryos showed low GFP expression and fewer GFP-positive embryos in the craniofacial tissues at 7,14, 24, and 30 dpf.