Skull development is a complex process that involves ongoing interaction between the bones of the skull and cranial soft tissues. The cranial vault is comprised of intramembranous bones joined by sutures of dense fibrous tissue that accommodate the growing brain. Bone is added at these sutures during growth, and the skull eventually ossifies completely. Craniosynostosis (CS), the premature closure of one or more of the cranial vault sutures is a common congenital anomaly. Approximately 80% of CS occurs as an isolated anomaly, non-syndromic craniosynostosis (NSC), without major associated malformations (1
). Rare mutations in the FGFR2, TWIST1
, FREM1, LRIT3, EFNA4,
duplications have been reported in a minor fraction of NSC cases (2
). The remainder (approximately 20%) of CS cases are syndromic, occurring with one or more additional major malformations, due to single-gene mutations in one of at least eight genes (FGFR1, FGFR2, FGFR3, TWIST1, EFNB1, POR, MSX2
) primarily involving the coronal sutures (10
The most common type of NSC is synostosis of the sagittal suture (sNSC), accounting for 40 to 58% of all cases (13
). The etiology of sNSC is not well understood, although the published literature suggests that it is a multifactorial condition, where both genetic and environmental factors play a role, as indicated by a higher rate of concordance in monozygotic as compared to dizygotic twins (30% vs. 0% in sNSC) (15
), increased male to female ratio (3:1) (1
), and a high recurrence risk in affected families (17
). Several environmental factors have also been associated with sNSC, including parity, prematurity, intrauterine constraint, and maternal tobacco or nitrosatable drug use (18
In an attempt to identify susceptibility loci for sNSC, we recruited 201 case-parent trios and 13 nuclear families with at least one individual affected with sNSC (68% non-Hispanic white [NHW], 32% mixed ethnicity families). These participants were enrolled and evaluated through the collaborative effort of the International Craniosynostosis Consortium (genetics.ucdmc.ucdavis.edu/icc.cfm
). Samples were genotyped on the Illumina 1M Human Omni1-Quad array. Final criteria for the current analyses were based on strict sNSC phenotype (no additional synostoses, congenital anomalies or developmental delay present in the proband or affected sibling) and NHW ethnicity as determined by self-reported ancestry and confirmed by principal-component analysis (PCA) using EIGENSOFT 3.0 (22
). Application of these criteria reduced our discovery population to 130 case-parent trios.
After stringent quality control procedures, we retained genotypes for 915,307 SNPs (914,402 autosomal and 905 pseudo-autosomal) which were analyzed using the transmission disequilibrium test (TDT) as implemented in PLINK v1.07 (23
). This resulted in the identification of 21 SNPs , 18 on 20p12.3 and 3 on 7p14.3, which reached the genome-wide significance threshold of P < 5 × 10−8
(, Supplementary Fig. 1
). The 18 chromosome 20 SNPs were in two LD blocks (92 kb and 19 kb, on chromosome 20p12.3; , Supplementary Fig. 2
, ), with the most significant SNP, rs1884302 (P
= 1.13 × 10−14
), located approximately 345 kb downstream of BMP2.
Conditional analysis using extended TDT (ETDT) in UNPHASED (24
) confirmed that rs1884302 is leading these signals (data not shown). Except for a predicted mRNA, BC038533
, no other known genes, CNVs, transcripts, miRNA or regulatory elements were identified in this region. However, at least four highly conserved genome regions with multiple predicted sites for transcription factor binding were identified by bioinformatic analysis in the vicinity of rs1884302 (Supplementary Fig. 3
). Exonic resequencing of BMP2
in our cohort did not identify any variants with obvious deleterious effects. The three SNPs on chromosome 7p14.3 span a region of 167 kb within BBS9
introns 4 and 15 (, Supplementary Fig. 4
). RT-PCR sequencing of cDNA from eight sNSC calvarial osteoblast cell lines discovered six BBS9
variants, but three were determined to be previously reported SNPs and all were benign based on bioinformatic protein prediction tools and conservation analysis (see online Methods). We also searched for two-way gene-gene interactions using the 100 most significant TDT GWAS SNPs using a conditional logistic regression model (25
). No results were significant after correcting for multiple testing (1247 tests, Bonferroni P
≤ 4 × 10−5
, Supplementary Table 1
Manhattan plot of P-values obtained from the genome-wide TDT of 130 trios (N=914,402)
Regional association plots for associations with sNSC at genome-wide significance
GWAS, replication study and meta-analysis results for most significant SNP on chromosomes 20 and 7
Considering the male predominance of sNSC, we performed a separate TDT analysis for 104 NHW trios with a male proband. The 21 SNPs on chromosomes 20 and 7 remained genome-wide significant, with an additional pseudo-autosomal marker reaching a low, but not genome–wide statistically significant P-value (rs2522623, P = 5.75 × 10−5) and thus did not meet our criteria for replication. This marker is located in protocadherin 11 X-linked (PCDH11X), which to date has not been reported to play a functional role in sNSC. Analysis of a larger cohort may demonstrate a stronger association to this locus.
From a total of 46 autosomal SNPs that showed significant and suggestive associations with P
< 1 × 10−5
, 25 SNPs were prioritized for the replication study (Supplementary Table 2
). In order to confirm the observed signals, we genotyped a NHW replication population of 172 unrelated cases with sNSC and 548 unaffected controls. We considered our data as a replication when the direction of effect of the alleles surpassing nominal significance was consistent between the discovery and replication datasets. All seven genotyped markers on chromosome 20 and four out of five markers on 7p14.3, with one SNP failing genotyping, were successfully replicated with meta-analysis P-
(Supplementary Table 2
, Supplementary Table 3
). There was little or no evidence of heterogeneity of the odds ratio in the two different sample sets for the eleven replicated SNPs. Additionally, genotypic TDT (gTDT) analysis indicated the effect of SNPs was consistent with an additive model (Supplementary Table 4
). We detected that, in most cases, transmitted minor alleles conferred a 2.36-4.38 fold risk of sNSC at the chromosome 20 markers, whereas the chromosome 7 transmitted minor alleles were associated with a decreased risk of sNSC (Supplementary Table 3
). We were unable to replicate other signals from the GWAS TDT with P
< 1 × 10−5
. These associations could be false positives or were missed due to differences between the study design and statistical methods used for the discovery and replication stages. We also relied on self-reported maternal ethnicity instead of ancestry informative markers for the replication cohort, which limited our ability to precisely assign ethnicity to the case-control study. The effect of gene-environment interaction or the relatively small sample of the replication cohort may also have affected the findings.
Imputation on borderline significant regions on chromosomes 3, 7, 8, 15, 22 was performed using BEAGLE v3.3.2 (26
) to identify additional genome-wide significant SNPs. Fifty-two intronic SNPs in DLG1
reached genome-wide significance (Supplementary Table 5
), but we previously failed to replicate this association in our case/control sample (Supplementary Table 2
). The lack of replication of a common variant may be due to our small sample size, differences in study design between our GWAS and replication datasets, or rare variants that contribute to the risk. Imputation was carried out using the 1000Genomes data for the reference panel, which does include some sequence variants with minor allele frequencies < 1%. Given that sNSC has a prevalence of 1 in 5000, we suspect that rare functional variants not identified by 1000Genomes may be involved and sequencing of the associated region will be necessary. No statistically significant signals were observed after the imputation at the other loci.
Using these modestly-sized populations we discovered and replicated common genetic variants in at least two loci that contribute to sNSC. The chromosome 20 locus is located in the vicinity of the BMP2
gene that plays a role in skeletal development. BMP2 is a member of the TGFβ superfamily that activates SMAD1, 5, and 8, and together with other osteogenic transcription factors such as RUNX2, ultimately regulating osteoblast development (27
An interesting observation in humans and animal models is that mutations with opposite effects in the same gene result in craniosynostosis (accelerated ossification) or parietal foramina/large fontanelles (deficient ossification) phenotypes. This holds true for mutations in at least 4 genes: MSX2
) and Nell1
) Consistent with this, gain-of-function mutations in ALX4
, a gene associated with parietal foramina, have been identified in sNSC probands by our group (36
). Similarly, a variant in BMP2
, p.S37A, has been associated with osteoporosis (37
): thus, if p.S37A results in loss-of-function, one can speculate that gain-of-function variants causing overexpression of BMP2
would result in a phenotype of accelerated ossification and potentially CS. Additionally of interest, a limb specific regulatory element located 110 kb downstream of BMP2
is mutated in brachydactyly type A2 (38
). Therefore, we predict that a skull-specific BMP2
regulatory element is located within or near the association region. Also, non-coding SNPs can affect gene expression variation of one or many genes. Even though none of our top hits were found in the eQTL databases, this may be due to the fact that no human eQTL studies have been conducted in the skull suture tissues.
The GWAS and the replication analysis also showed evidence for association on chromosome 7p14 within an intronic region of BBS9.
BBS9 is a component of the multiprotein complex BBSome which plays critical roles in intraflagellar transport, or moving cargo molecules in and out of cilia. Loss-of-function mutations in BBS9
have been found in patients with Bardet-Biedl syndrome (BBS) (39
). BBS arises from defects of primary cilia that function as a platform for many important signaling events involving at least PDGFRα, hedgehog, epidermal growth factor, and 5-HT6
serotonin receptor (40
). Although BBS does not present with suture phenotypes, it presents with a skeletal phenotype (polydactyly) and orofacial anomalies (39
). It is noteworthy that aberrant signaling of PDGFRα and epidermal growth factor or hedgehog has been implicated in craniofacial development (45
). In addition, it has been shown that FGF-signaling regulates cilia length and function during development of zebrafish and Xenopus
). Thus, it is possible that CS-causing FGFR gain-of-function mutations also work via aberrant assembly and/or signaling of cilia. It remains to be seen whether gain-of-function mutations or overexpression of BBS9
would influence suture fusion. Three of the genes in the vicinity of the association on chromosome 7 (RP9P
) are pseudogenes and the remaining two, RP9
, have no known role in skeletal development and cause retinitis pigmentosa and hemolytic anemia, respectively (51
). This region is located approximately 14 Mb telomeric of TWIST1
and thus is unlikely to impact its function.
Our findings implicate variants within or near genes with relevant function in skeletal development, suggesting a major genetic effect for sNSC. A larger sample size will be needed to detect additional susceptibility genes with moderate genetic effects. The lack of Mendelian segregation suggests a more complex pattern of inheritance for sNSC. The causative variants may be of low frequency and in linkage disequilibrium with the common alleles associated with sNSC in this study. Another possibility is that one or more environmental factors interacting with genetic variants contribute to sNSC. Less likely, but still possible are oligogenic inheritance, gene-gene interactions, imprinting, epigenetic changes, or incomplete penetrance of the variants.
Biologically plausible candidate genes, BMP2 and BBS9, both involved in skeletal development, were identified in our study. Functional studies of these genes, such as RNA sequencing and possibly tiling arrays of the regions to look for non-coding or microRNAs, are warranted. The potential contribution of these loci to other types of NSC is worthy of investigation. Our results represent the first major step toward deciphering the genetic etiology of the sNSC.