Introduction

Craniosynostosis is the result of premature fusion of one or more cranial sutures. The prevalence is 1 per 2000 to 2500 live births. All cranial sutures could be affected but with different frequencies. Scaphocephaly defined by premature closure of the sagittal suture is the most common craniosynostosis, with a frequency of 41% to 54% [1,2,3]. The frequency of trigonocephaly is approximatively 15% but has increased during the last 2 decades, both in Europe and the United States [4, 5].

In coronal craniosynostosis, an underlying genetic cause can be found in 60% of non-syndromic bicoronal cases and 30% of unicoronal cases; in most of the isolated sagittal craniosynostosis cases, the molecular etiology is unknown [1]. Some rare variants were identified in SMAD6 [6, 7], FGFR2 [8, 9], AXIN2 [10] and TWIST1 [11]. Isolated trigonocephaly is caused by loss-of-function variants in SMAD6 [6, 7] and FREM1 [12].

We aimed to determine the molecular etiology in children affected by isolated sagittal and metopic craniosynostosis based on the known genes in these defects [4, 13]. Our approach was based on target panel sequencing in a cohort of infants who underwent craniosynostosis surgery. Children with genetic craniosynostosis were followed after surgery to evaluate their evolution, including neuropsychological development. Our results highlight the use of molecular analysis in the context of isolated sagittal and/or metopic craniosynostosis to enhance an understanding of the pathophysiology of midline craniosynostosis.

Patients and methods

Patients

Our cohort included all infants with isolated midline synostosis (sagittal or metopic) who were followed in the craniosynostosis national reference center at Lyon University Hospital from 2018 to 2020. All infants presented a clinically and radiologically confirmed premature fusion of the sagittal or the metopic suture. No teratogenic exposure was reported in all cases. infants with clinical and radiological signs of syndromic craniosynostosis were excluded. All infants were examined by a multidisciplinary team including neurosurgeons, clinical geneticists and neuropsychologist. The molecular analysis was approved by the French ethics committee (no. IRB00011687; College de neurochirurgie IRB #1: 2022/05) and was performed after parents gave signed informed consent for genetic testing.

Genetic analysis

Next-generation sequencing (NGS)

DNA was extracted from whole blood samples by using the Qiasymphony instrument (Qiagen, Courtaboeuf, Paris) according to the manufacturer’s protocol. All DNA was screened by targeted-NGS with a panel of 31 genes associated with craniosynostosis (ALPL, ALX3, ALX4, ASXL1, AXIN2, BCL11B, CDC45, EFNA4, EFNB1, ERF, FGF10, FGFR1, FGFR2, FGFR3, FREM1, GLI3, IFT122, IL11RA, MSX2, PHEX, RECQL4, RUNX2, SMAD6, SPECC1L, SIX1, STAT3, TCF12, TWIST1, WDR19, WDR35, ZIC1). The targeted regions were exons and splicing junctions at +/−50 pb designed from hg19/GRCh37 genome. The NGS technology was based on the SureSelectQXT Reagent kits (Agilent, Les Ulis, France) which provide DNA library preparation solutions with hybrid capture-based target enrichment. They include all reagents for library preparation, hybridization, and capture. The post-capture pooling of enriched libraries was sequenced on a Miseq sequencer (Illumina, Paris). The MiseqReporter (Illumina) integrated in the Miseq sequencer was used to perform the variant calling and annotation. Sequence results were obtained after aligning the fastqs, mapping, variant filtering and calling by using SeqNext software (JSI Medical Systems, Ettenheim, Germany) based on Smith-Waterman and Burrows Wheeler Aligner algorithms. The copy number variation (CNV) was also calculated with the CNV module of this software. In detail, each run contains control sequences that allow the internal normalization of the coverage of the analyzed sequences within the same run. The number of reads covering each region of interest (ROI) of patient DNA is normalized by the median of the reads number of the control-ROI present in the same run, resulting in the relative coverage for the target ROI of patient sample (ROI P). The same calculation is performed for the control sample (without CNV), resulting to the relative coverage for the target ROI of control sample (ROI C). The ROI P is normalized by the ROI C (ROI P/ROI C), resulting in the ratio relative coverage (ratio ROI). The ratio ROI is then stated as a percentage (%). The coverage of the targeted sequencing was >99.9% for depth 30X. The annotated variants were filtered according to minor allelic frequency (MAF), and variants with maximal MAF < 1% were retained by using gnomAD (https://gnomad.broadinstitute.org/).

The others were considered polymorphisms and were removed. The unknown variants with maximal MAF < 1% were classified by using the Alamut interface (SOPHiA Genetics, Lausanne, Switzerland) with SIFT (http://sift.jcvi.org) and Polyphen 2 (http://genetics.bwh.harvard.edu/pph2). The variant pathogenicity was estimated by using CADD (https://cadd.gs.washington.edu). The analysis of splicing variants involved using MaxEntScan, NNSPLICE and GeneSplicer.

If necessary, the variants were searched in the databases ClinVar (http://www.ncbi.nlm.nih.gov/clinvar) and professional Human Gene Mutation Database (http://www.hgmd.cf.ac.uk/ac/index.php). Variants were interpreted in accordance with the American College of Medical Genetics and Genomics guidelines (Richards et al., 2015): all heterozygous identified variants were classified in 3 classes defined as variants of unknown significant (VUS) or likely pathogenic or pathogenic variants.

Sanger sequencing was used to confirm the identified variants of interest by NGS and to define the RT-PCR sequence by using Life Technologies products and software on an ABI3130 sequencer (ThermoFischer, Les Ulis, France).

Transcriptional study of splicing variants

The splicing variant study for SMAD6 was performed from paxgene blood samples. RNA was extracted with the RNeasy kit (Qiagen, les Ulis, France), followed by the high-capacity cDNA Reverse Transcription Kit (Thermofisher). cDNA was amplified with primers designed in exons 1 and 3 of SMAD6 (ex1F:CTCAAAACCGTCACGTACTCG and ex4R:CTGCCCTGAGGTAGGTCGTA). The PCR products were extracted from gel and sequenced by the Sanger technique.

Array-comparative genomic hybridization (CGH) analysis

The array-CGH involved using the SurePrint G3 human CGH Microarray Kit 180 K according to the manufacturer’s instructions (Agilent Technologies, France). The analysis for copy number changes involved Agilent CytoGenomics v.5 (Agilent Technologies) with genome version hg19. CNV was classified according to the Database of Genomic Variants (http://projects.tcag.ca/variation) to exclude common polymorphic CNVs with frequency >1% in healthy controls.

Neuropsychological testing

Neuropsychological data in children with a genetic variant were compared with those of a control group of 70 children with a single sagittal suture or metopic suture without identified genetic anomalies or co-morbidities, who benefitted from a systematic neuropsychological evaluation with the same tests during the same period in our craniofacial unit. All children with a positive genetic finding underwent neurodevelopmental evaluation. The type of test was adapted to the age of the child at the time of the neuropsychological assessment.

Neuropsychological development was assessed in all children before age 6 years with the French adaptation [14] of the Child Development Inventory (CDI) of Ireton. Briefly, the CDI allows for obtaining information from the parents about their child’s development. The CDI is designed to identify developmental risks in young children and help identify domains for which prevention actions could provide appropriate help, as parents are encouraged to become more involved in observing and understanding their child’s development. Several domains are isolated and provide several specific scores: social (SO), self-help (SH), gross motor (GM), fine motor (FM), expressive language (LEX), language comprehension (LCO) and a general development (GD).

Moreover, children over 2.5 years of age, were also tested with Wechsler’s scales, WPPSI-IV or WISC-V according to their age.

GD and IQ were pooled for the comparative analysis due to their high correlation of 0.85 (Duyme, 2010a) and the ratio of actual development to theoretical development calculated. We considered a ratio <85% as indicator of a high risk of delayed development; and a ratio <70% of a very high risk of delayed development (Duyme, 2010a, Ireton).

Results

Characteristics of the cohort

We included 101 infants with isolated trigonocephaly or scaphocephaly. No infant had coexisting metopic and sagittal craniosynostosis. All included craniosynostosis cases were isolated without any other significant phenotypic features. The M/F sex ratio was 1.89. We included 58 infants with scaphocephaly (38 boys) and 43 with trigonocephaly (28 boys). Craniosynostosis, whether metopic or sagittal synostosis, was more frequent in boys than girls. All patients underwent cranial remodelling before age 12 months, then were followed yearly. The mean follow-up was 24.8 months.

Molecular analysis of the cohort

Among our 101 infants, 13 carried a total of 13 variants; that is, 12.9% of the infants carried a variant in genes known to be involved in craniosynostosis.

The M/F sex ratio was 1.6, with 5 scaphocephaly cases (2 boys) and 8 trigonocephaly cases (6 boys).

Seven infants carried variants in SMAD6, 2 in FGFR2, 1 in TWIST1, 1 in FREM1, 1 in ALX4 and 1 in TCF12. All variants were detected at the heterozygous level in genes associated with an autosomal dominant craniosynostosis. No homozygous or compound heterozygous genotypes were identified. All retained variants had maximal MAF < 0.025% and all variants classified as likely pathogenic presented a gnomAD MAF below <2 × 10−5. (Tables 1, 2).

Table 1 characteristics of the patients with SMAD6 variants.
Full size table
Table 2 characteristics of the patients with variants other than SMAD6.
Full size table

To complete the analyses of the variants identified in affected infants, we used the observed/expected score (oe) suggested in gnomAD. This score is important to interpretate the pathogenic effect of loss-of-function variants because the lower the oe, the less tolerance (threshold <0.35).

SMAD6 association with midline craniosynostosis

In male infants affected by a trigonocephaly, we identified 2 frameshift variants in SMAD6 that resulted in premature stop codons (NM_005585.4: c.465_471del; NM_005585.4: c.1296dupC). Both variants were classified as “pathogenic” even if the oe was high (1.97). This gene is clearly associated with metopic premature synostosis. Indeed, this variant has been described as pathogenic in craniosynostosis [6] and was identified in radioulnar synostosis [15] and bicuspid aortic valve [16]. Moreover, the nonsense variant NM_005585.4: c.43 C > T, p.(Arg15*) was identified in a girl with trigonocephaly; this variant was reported in premature fusion of metopic suture [7], so the variant was considered pathogenic.

Three missense variants (NM_005585.4: c.793 C > T, p,(His265Tyr); c.817 G > A, p.(Glu273Lys) and c.857 A > G, p.(Asp286Gly)) were found in trigonocephaly cases. The 3 variants had high CADD Phred scores of 25 and 27, beyond the threshold of 20, which is considered damaging. The two variants (p,(His265Tyr) and p.(Glu273Lys)) characterized by a very low MAF (<0.000045) were predicted as deleterious by multiple computational algorithms including predictions from CADD, Polyphen-2 and SIFT, MetalLR, Fathmm and REVEL. Both of these variations affect highly conserved amino acids located in the MH1 domain which may responsible for an instability of the protein as shown by Calpena et al. for other closed missense variants [6]. In addition, the p.(His265Tyr) residue is involved in zinc metal binding required for transcription factor activity in UniProt (https://www.uniprot.org/) and at the nucleotide level, the c.817 G > A is found in the exon-intron junction that could affect splicing. Therefore, we considered these variants as probably pathogenic according to ACMG guidelines. The NM_005585.4: c.857 A > G, p.(Asp286Gly) variant presents a higher MAF than the previously mentionned SMAD6 variants and therefore, exceeds the provided AFmax threshold (0.000045) recommend in the publication of Calpena et al.5. In addition, the p.Asp286 residue is not conserved and is located in a low conserved region. The prediction softwares (CADD, Polyphen-2 and SIFT, MetalLR, Fathmm and REVEL) classified it as “damaging”. Therefore, this variant was considered as a VUS (Table 1).

In one case of scaphocephaly, we detected an intronic variant that affects one of the 4 highly conserved nucleotides and predicts leading to aberrant splicing at the donor site (NM_005585.4: c.874 + 4 A > G) according to in silico prediction software (Maxent, Genesplicer, NNsplice). A functional study confirmed the impact of the splicing variant from blood samples. The result showed an exon 2 skipping at the heterozygous level as demonstrated in Fig. 1. The p,(His265Tyr) and p.(Glu273Lys) variants were inherited by the mother who was clinically asymptomatic, whereas the analysis segregation for the third patient (NM_005585.4: c.857 A > G, p.(Asp286Gly)) was incomplete because the mother was negative and the father was unavailable Moreover, segregation analysis was not informative because the parents were clinically asymptomatic. Also, this variant was considered probably pathogenic according to results of functional study and the possibility of incomplete penetrance. Finally, among the 7 identified variants in SMAD6, except for NM_005585.4: c.465_471dup and NM_005585.4: c.43 C > T, we found 5 variants that were novel (Table 1).

Fig. 1: Transcription analysis.
figure 1

Effect of NM_005585.4: c.874 + 4 A > G on splicing (A) Maxent prediction software: decreased use of the expected donor site and increased use of a potential downstream donor site, (B) electrophoresis of RT-PCR product showing one band for control and two bands for the patient 5, (C) Sanger sequencing showed that the variant led to exon 2 skipping of SMAD6 transcript.

Full size image

FGFR2 pathway associated with scaphocephaly

In the FGFR pathway, we found p.(Val274Ile) and p.(Pro23His) variants located in one of 3 Ig-like C2 in FGFR2 which have been already described in syndromic and isolated craniosynostosis [9, 17]. These variants were considered deleterious according to all in silico software analyses (CADD, polyphen-2, SIFT, REVEL and MetaLR) and their MAF was very low, especially for p.(Val274Ile) variant. In addition, the p.Val274 residue has been previously implicated in Crouzon syndromic with the following variants: NM_000141.5:c.820_824delGTAGAinsTT p.(Val274_Glu275delinsLeu) [17] and NM_000141.5: p.(Val274_Glu275insGlyGlyAspLeu) [18] suggesting a probable pathogenicity. Furthermore, these variants are found in one of the two IgIII domains most involved in syndromic craniosynostosis. However, the parents of the proband showed no clinical signs, So, p.(Val274Ile) and p.(Pro23His), associated with scaphocephaly were classified as VUS in craniosynostosis.

A variant was found in TWIST1 that codes for a transcription factor with an inhibitor effect on the FGFR2 pathway [19]. This p.(Gly57Ala) variant in TWIST1 was never listed in any database and was considered probably pathogenic according to some in silico software analyses (Table 2). Exceptional missense variants closed to the identified variant and outside the bHLH domain were already repertoried [11, 20]. In addition, the mother’s proband was negative for this variant and his father died before the molecular diagnosis Therefore, according to ACMG guidelines, we considered this variant to be classified as VUS.

Moreover, we identified a heterozygous NM_207036.1: c.971-8 A > G variant in TCF12 in a patient with scaphocephaly. This gene codes for a transcription factor that can heterodimerize to TWIST1. The intronic variant in TCF12 is located at an acceptor site and has a pathogenic effect according to in silico prediction software analyses (Maxent, Genesplicer, NNsplice). Indeed, this variant abolished the use of the expected splicing site in c.971 with maximal score −100% and created a new splicing acceptor site in c.971-7 (Table 2). The family analysis was not contributive.

Another heterozygous missense p.(Arg216Trp) variant was identified in ALX4 in an isolated case of scaphocephaly. It is located at the same position as the NM_021926.3: c.646 C > T, p.(Arg216Gly) variant, which is associated with a moderate form of frontal dysplasia type 1 with minute parietal foramina [21]. The presence of enlarged parietal foramina is explained by a loss-of-function variant [22], and the craniosynostosis is due to the gain-of-function variant [23, 24]. Although this p.(Arg216Trp) variant in ALX4 is considered deleterious according to all in silico software analyses (CADD, polyphen-2, SIFT), the heterozygous p.(Arg216Trp) variant inherited by the unaffected mother of the infant was classified as VUS in this context. The index case and his father carried the variant NM_006492.2:c.510 G > C, p.(Lys170Asn) in ALX3. This variant presented very low frequency in gnomAD (all populations: 0.0016%; African: 0.018%; European: 0.00090%) and was considered probably pathogenic according to prediction software analyses (SIFT, CADD, Polyphen-2). Both paralogous ALX3 and ALX4 are involved in craniofacial development [25] and both could be targets of TWIST1 [26, 27]. So, variants identified in scaphocephaly could belong to genes involved in the FGFR2 pathway.

FREM1 associated with trigonocephaly

We identified a large heterozygous deletion, NM_144966.5: c. (? _ 1); (1393 + 1_1394-1)del, in FREM1 in an infant with trigonocephaly. The deletion that affected exons 1 to 9 of FREM1 was confirmed by array-CGH: Arr(GRCh37) 9p22.3(14848725-14986398)x1 with size 120 kb. This deletion including only FREM1 was inherited by the infant’s asymptomatic father. Also, this deletion was previously described as a risk factor in isolated trigonocephaly [12].

Neuropsychological testing

Neurodevelopmental assessment involved a neuropsychological test after surgery in 13 children aged 18 months to 6.8 years. Mean age at neuropsychological testing was 3.2 years (SD 1.9). Patients who completed the CDI were 18 months to 5 years, 7 months old (mean 2.0 years [SD 0.9]). Patients who completed the Wechsler Scales were ages between 3.2 years and 6.8 years old (mean 4.2 years, sd 1.6 years).

The control group included 70 children. In total, 51 children had scaphocephaly (40 boys; age 2 to 16.5 years; mean 6.7 years [SD 3.0]) (Table 3). Control children with scaphocephaly who completed the CDI were 1.5 to 4.5 years (mean 3.5 years [SD 1.2]) and those with scaphocephaly who completed the Wechsler scales were 3.8 to 16.5 years (mean 7.9 years [SD 2.5]). A total of 19 control children had trigonocephaly (11 boys; age 2 to 16.5 years; mean 6.4 years [SD 2.6]). Children with trigonocephaly who completed the CDI were 1.7 to 5.7 years old (mean 3.9 years [SD 3.0]) and those with trigonocephaly who completed the Wechsler scales were 5.7 to 11.7 years (mean 7.4 years [SD 1.7]).

Table 3 Distribution of children by percentiles of development.
Full size table

We analyzed the distribution of children as a function of their intellectual level, assessed as well by GD or IQ, between 4 percentiles as represented in Table 3.

Children with variants were more frequently under the 50th pecentile of development than those without variants (p < 0.05) (Table 4).

Table 4 Distribution of children with and without variants by percentiles of development.
Full size table

The rate of neuropsychological development under the 50th percentile was higher for children with than without variants (Fig. 2).

Fig. 2: Distribution of children with or without a variant according to their developmental level, below or under the 50th percentile.
figure 2

Unlike the control group without mutation, the variant group is mostly under the 50th percentile.

Full size image

We used the 15th percentile of development as a cutoff to observe different profiles. Children with scaphocephaly and development under the 15th percentile were close to 15% for the global development as in the normal population (Table 5), but more frequently under the 15th percentile for fine motor. For children with trigonocephaly, developmental was frequently under the 15th percentile for language.

Table 5 For the 2 groups (scaphocephaly with variants and trigonocephaly with variants), number and ratio of children under the 15th percentile for Child Development Inventory (CDI) scores.
Full size table

For our variants group, early development is more sensitive to difficulties for expression acquisitions, self help and also fine motor (Table 6). However, children with SMAD6 and trigonocephlay show more difficulties specifically for langage and self help.

Table 6 For the group with SMAD6, number and ratio of children under the 15th percentile for Child Development Inventory (CDI) scores.
Full size table

Discussion

We evaluated the presence of variants in genes described in isolated sagittal and metopic synostosis in infants who had undergone surgery at the craniosynostosis national reference center of Lyon University Hospital from 2018 to 2020. A total of 101 infants with midline craniosynostosis were studied. The male preponderance found in midline craniosynostosis agrees with the Norwegian study [28]. The overall variant detection frequency was 12.9%. Also, we showed a higher variant frequency in boys affected by isolated sagittal or metopic synostosis than girls.

The SMAD6 gene emerged as the major gene in isolated trigonocephaly, with a positive result in 7/43 (16.2%) children. It was the third independent cohort showing a substantial number of metopic cases associated with SMAD6 variants [6, 7]. Six variants can be classified as probably pathogenic and participate in the phenotype. For the seventh, only further screening studies could provide information to confirm or deny its pathogenic role on the median suture.

We observed an incomplete penetrance in trigonocephaly associated with SMAD6. These observations agree with the published results of The Oxford team in which the penetrance reached only 16% [6]. The segregation analyses were performed for all recruited infants affected by an isolated form and who carried a variant in a gene involved in craniosynostosis; however, we observed extreme phenotypic variability or incomplete penetrance in most families. This well-known observation in isolated craniosynostosis underlines the difficulties in interpreting the segregation analysis in the presence of variants and highlights the importance to perform tests in unaffected parents and relatives for genetic counseling. Also, a digenic inheritance [29], genetic modifiers [7] or environmental and epigenetic factors [30] might be suggested to explain the incomplete penetrance. Studies in this research domain need to be led to understand this common characteristic of non-syndromic craniosynostosis.

For children with a SMAD6 variant, the neurodevelopmental evolution showed frequently a delayed acquisition of language, both in expression and comprehension, even if self-help or global development seemed unaffected. A neuropsychological follow-up should be recommended in children with SMAD6 variant and some help could be provided in the early stages of language acquisition before the lecture learning.

Our findings, in accordance with the observations made by Wu et al. [31], confirm that SMAD6 variants has an impact on the neurodevelopment. However, our results differ in the type of affected profiles probably due to the young ages of our patients. We can thus hypothesize that the impact of SMAD6 variant may vary according to the age and to the stage of neurodevelopment of the child.

The deletion in the FREM1 gene in trigonocephaly was described in only one publication with a mice mouse model that confirmed the impact of FREM1 in metopic fusion [12]. However, a recent publication contradicted its pathologic role in trigonocephaly and considered that FREM1 has no effect on metopic suture [32]. Actually, the deletion should be considered a genetic susceptibility factor.

Concerning FGFR2, we identified 2 VUS associated with scaphocephaly; this phenotype is exceptional and was previously described [8, 9]. Interestingly, the p.(Pro23His) variant is located in the same region as the p(Cys62Arg), both associated with a scaphocephaly [9]. It seems therfore that a missense variation in this specific domain could lead to a mild phenotype with an extreme clinical variability or incomplete penetrance as already observed by McGillivray et al. for a variant located in tyrosine kinase domain [8]. Differently, with respect to the p.(Val274Ile) variant, its position, frequency and prediction of pathogenicity support that this variant is probably responsible for the phenotype. Thus, the FGFR2 pathway seems not involved in fusion of metopic suture as compared with the SMADs pathway.

In addition, heterozygous missense variants in TWIST1 located outside the bHLH domain and closed to the variant we identified were already associated with Saethre Chotzen syndrome [11, 20]. Also, heterozygous missense variants in TWIST1 could be responsible for non-syndromic craniosynostosis with premature fusion of sagittal suture [33]. In our case, the incomplete penetrance and low pathogenicity of the missense variant not listed in the databases suggest that this variant could be a susceptibility factor in craniosynostosis. Further functional studies of TWIST1 missense variants outside the bHLH domain could be interesting to evaluate their pathogenicity as previously done by Funato et al. [20].

The ALX4 variants associated with the craniosynostosis type scaphocephaly are described as gain-of-function variants [23, 24] that emphasize the activity of the ALX4 transcription factor in osteoblast differentiation [34]. Currently, nine probably pathogenic or VUS associated with sagittal non-syndromic craniosynostosis have been described in literature (Supplementary Table) by different research teams [23, 24, 35,36,37]. As with the other non- syndromic craniosynostosis, an incomplete penetrance characterized ALX4 gene. In our study, the p.(Arg216Trp) variant located in DNA binding of ALX4 appears to be associated with craniosynostosis. At this residue, the p.(Arg216Gly) variant in ALX4 was found associated with a moderate form of frontonasal dysplasia type 2 with minute parietal foramina [25]. However, this opposite phenotype was not noted in our patient 13. It is important to consider that for a given amino acid such as Arginine, located in a DNA binding domain, the effect differs depending on the amino acid modification. For instance, the substitution Arg for Gly is unable to bind to DNA compared to the substitution Arg for Trp which retains a strong DNA binding capacity [38]. Thus, this observation could be a clue to understand the impact of the variation we identified in patient 13. Especially since ALX4 may be a potential target of TWIST1 or even TCF12, both of which are required for calvarial development and suture closure [26, 39].

In conclusion, our results highlight the importance of genetic testing in children presenting isolated midline craniosynostosis for better understanding the pathophysiology of craniosynostosis and confirm the major role of SMAD6 in metopic non-syndromic craniosynostosis.