Results from Genetic Studies in Patients Affected with Craniosynostosis: Clinical and Molecular Aspects

Background: Craniosynostosis (CS) represents a highly heterogeneous genetic condition whose genetic background has not been yet revealed. The abnormality occurs either in isolated form or syndromic, as an element of hundreds of different inborn syndromes. Consequently, CS may often represent a challenging diagnostic issue. Methods: We investigated a three-tiered approach (karyotyping, Sanger sequencing, followed by custom gene panel/chromosomal microarray analysis, and exome sequencing), coupled with prioritization of variants based on dysmorphological assessment and description in terms of human phenotype ontology. In addition, we have also performed a statistical analysis of the obtained clinical data using the nonparametric test χ2. Results: We achieved a 43% diagnostic success rate and have demonstrated the complexity of mutations’ type harbored by the patients, which were either chromosomal aberrations, copy number variations, or point mutations. The majority of pathogenic variants were found in the well-known CS genes, however, variants found in genes associated with chromatinopathies or RASopathies are of particular interest. Conclusion: We have critically summarized and then optimised a cost-effective diagnostic algorithm, which may be helpful in a daily diagnostic routine and future clinical research of various CS types. Moreover, we have pinpointed the possible underestimated co-occurrence of CS and intellectual disability, suggesting it may be overlooked when intellectual disability constitutes a primary clinical complaint. On the other hand, in any case of already detected syndromic CS and intellectual disability, the possible occurrence of clinical features suggestive for chromatinopathies or RASopathies should also be considered.


INTRODUCTION
Premature fusion of calvarial sutures, i.e., craniosynostosis (CS), represents a highly heterogeneous neurocranium malformation (Morriss-Kay and Wilkie, 2005). The disease can be classified following at least three criteria: the origin, the number of affected sutures, and the occurrence of additional clinical features. Regarding the etiology, CS is distinguished as a primary condition, i.e., genetic; or secondary, if it arises from mechanical, metabolic, or hormonal defects. The assessment of affected sutures' numbers allows recognising single or compound CS, and finally, a syndromic CS is described when additional features other than secondary to premature suture fusion occur. Otherwise, CS should be classified as an isolated condition (Wilkie et al., 2007;Johnson and Wilkie, 2011;Lattanzi et al., 2017).
The neonatal skull comprises six calvarial sutures-one metopic, one sagittal, two coronal, and two lambdoid, closing physiologically from 3 months to late 50 years of postnatal life. Calvarial sutures are fibrous junctions, which allow the skull to grow and develop during the expansion of the brain and permit skull compression during delivery (Baer, 1954;Opperman, 2000;Rice, 2008). Thus, the lack of physiological sutural patency impedes the allometric cranial growth, resulting in cranial deformities and, frequently, increased intracranial pressure, i.e., craniostenosis. Consequently, affected patients present with facial dysmorphism, cortex lesion, seizures, intellectual disability, visual and hearing impairments, or breathing difficulties that are all secondary to CS (Renier et al., 1982;Thompson et al., 1995;Tubbs et al., 2001;Gupta et al., 2003;Mathijssen and Arnaud, 2007;Chieffo et al., 2010).
CS affects approximately one in 2,500 births and burdens public health due to the requirement of extensive surgical treatment in the first year of life and multi-level specialist medical care in the subsequent postnatal periods (Wilkie et al., 2010;Lattanzi et al., 2017). Despite recent advancements in genetic diagnostics, the pathogenesis of CS remains still unknown or partially understood. The large cohort screenings reveal genetic etiology in barely 21-62% of all recruited cases, depending on the size of the study, ethnicity of the population, and range of the molecular analysis. Conversely, about 40-80% of CS cases remain molecularly unresolved (Roscioli et al., 2013;Paumard-Hernandez et al., 2015;Timberlake et al., 2016;Lee et al., 2018;Topa et al., 2020). In this paper, we have presented the study results encompassing 166 individuals in whom we had applied a three-step diagnostic algorithm to identify different mutation types. In addition, we have also performed a statistical analysis of the clinical data we had obtained.

PATIENTS AND METHODS
All procedures involving human participants were performed under the Helsinki Declaration. The Institutional Review Board of Poznan University of Medical Sciences granted ethics approval (no. 742/17). All patients agreed to participate in this study, and written informed consent for participation and publishing the information and images in an online open-access publication was obtained from all participants and the parents of minors before genetic testing.

Cohort Description
The patients were recruited provided they were born from pregnancies without exposure to environmental factors potentially causative for CS. Our cohort consisted of 166 individuals (33 patients belonging to 18 families and 133 sporadic patients), of whom 85 were males and 81 females. All patients underwent dysmorphological assessment, which allowed us to section off the two subgroups-isolated (if CS was accompanied only by the secondary defects directly resulting from CS) and syndromic (if CS was accompanied by additional defects, not resulting from CS).

Statistical Analysis
STATISTICA (version 13.3) TIBCO software was used for data analysis. The statistical significance of the phenotypic diversity of the cohort and the differences between the frequency of identified genetic modification in the different phenotypic groups was tested using the nonparametric test χ 2 .

Genetic Studies
We extracted genomic DNA from peripheral blood leukocytes drawn into EDTA-coated tubes using either the manual salting-out method or automated extraction using the MagCore ® HF16 Automated Nucleic Acid Extractor (RBC Bioscience Corp.). Whole blood for lymphocyte culture was drawn into heparin-coated tubes. The study has been divided into three tiers. Tier 1 included karyotyping and screening of the most frequent mutations located in exon no. 7 of FGFR1 (NM_023110.3), exons no. 7 and 8 of FGFR2 (NM_000141.5), and exon no. 7 of FGFR3 (NM_000142.5), and the entire coding sequence of TWIST1 (NM_000474.4). In tier 2, chromosomal microarray analysis and targeted nextgeneration sequencing (NGS) of custom gene panel were performed. Exome sequencing (ES) was applied in Tier 3.

PCR and Sanger Sequencing
We performed molecular screening for all recruited patients (n = 166) utilizing targeted PCR followed by Sanger sequencing. We tested the occurrence of the most frequent, recurrent mutations located in exon 7 of FGFR1 (NM_023110.3), exons 7 and 8 of FGFR2 (NM_000141.5), and exon 7 of FGFR3 (NM_000142.5), and the entire coding sequence of TWIST1 (NM_000474.4). Specific primers for amplification were designed using the online available Primer3 tool v. 0.4.0. For the detailed list of primers, see Supplementary Table S1. PCR products were sequenced using dye-terminator chemistry (kit v.3, ABI 3130XL) and run on automated sequencer Applied Biosystems Prism 3700 DNA Analyzer.

Karyotyping
Whole blood lymphocyte culture was performed following the standard protocol. Next, we used the Giemsa-banding (GTG) technique at 550 band resolution per haploid genome.

Chromosomal Microarray Analysis
Patients suspected of harboring pathogenic copy number variations (CNVs) were tested using (CMA). Depending on the research stages', different CMA formats were used. Firstly, the assay was performed employing a high resolution 1.4 NimbleGen oligonucleotide array comparative genomic hybridization (aCGH; Roche NimbleGen) according to standard protocols provided by the manufacturer. Results were analysed with Deva software (Roche Nimblegen) using the ADM2 segmentation algorithm (Sowińska-Seidler et al., 2018). The chromosomal profile was visualized using SignalMap software (NimbleGen Systems Inc.). Next, we applied SurePrint G3 Human CGH Microarray 8 × 60 k, 4 × 180 k, 1 × 1M arrays (Agilent Technologies). The hybridization signals were detected with SureScan Dx Microarray Scanner (Agilent Technologies) and visualized with the use of Agilent CytoGenomics software (Agilent Technologies) (Bukowska-Olech et al., 2020b). The pathogenicity of CNVs was evaluated using the following tools and available online databases-Cytoscape 3.7.1, ClinGen, DECIPHER, database of Genomic Variants (DGV), Mouse Genome Informatics (MGI), or UCSC Genome Browser applying tracks such as Conservation, VistaEnhancers, ENCODE Regulation or HiC.

NGS of Custom Gene Panel
A cohort with negative results was subjected to targeted nextgeneration sequencing of a custom 225.709 kb in size gene panel (Agilent Technologies). Captured and indexed libraries were sequenced on the previously described Ion Torrent S5 sequencing system. Variants identified by TorrentSuite, as first described by Bukowska-Olech et al., were further analyzed using an extended custom pipeline (Bukowska-Olech et al., 2020c). For prioritizing candidate variants in probably causative genes, local Phen2Gene installation was used. For final SNV/indel prioritization, the updated pipeline combined Exomiser 12.1.0 with ANNOVAR (all noncommercially available databases, downloaded on 16th December 2020), Ensembl/VEP 102.0, and CADD 1.6 Exomiser default phenotype scoring was supplemented with alternate scoring where the original formula used Phen2Gene scores instead (Zhao et al., 2020). Two prioritizers (OMIM and HIPHIVE) were used with Exomiser. For both Phen2Gene and Exomiser, each sample was labeled with Human Phenotype Ontology terms assigned according to clinical notes (manual curation) (Supplementary Table  S2). Common (AF>0.001 in population frequency datasets) and benign (as per strong ClinVar support) alleles were dropped out for reporting. The final pathogenicity of detected variants was analysed in line with the American College of Medical Genetics (ACMG) classification (Richards et al., 2015). Confirmation and segregation studies were performed applying PCR followed by Sanger sequencing as described in section PCR and Sanger sequencing. A list of primers was summarized in Supplementary Table S1.

Exome Sequencing
The coding region and flanking intronic regions were enriched using a custom-designed in-solution exome enrichment (TWIST bioscience, San Francisco, United States) and were sequenced using the Illumina NovaSeq system (Illumina, San Diego, United States). Sequencing reads were demultiplexed using Illumina bcl2fastq2. The removal of the adapter was performed with Skewer. The trimmed reads were mapped to the human reference genome (hg19) using the Burrows-Wheeler Aligner, and variants were called using in-house software. First, Only SNVs and small indels in the coding regions and the flanking intronic regions (±8 bp) with a minor allele frequency (MAF) < 1.5% were evaluated. Second, the known diseasecausing variants (according to Human Gene Mutation database; HGMD) were also evaluated in up to ±30 bp of flanking regions and up to 5% MAF. Downstream analysis was carried out using pipeline described for Tier 2 above. As before, the variant evaluation was based on the ACMG guidelines for interpreting sequence variants. Confirmation and segregation studies were performed applying PCR followed by Sanger sequencing as described in section PCR and Sanger sequencing. A list of primers was summarized in Supplementary Table S1.

Tier 1
PCR followed by Sanger sequencing of the most frequent mutations located within the FGFR1, FGFR2, FGFR3 genes and screening of the entire TWIST1 coding sequence allowed us to diagnose 43 patients from 32 different families. Out of them, three patients carried the same alteration in the 7 th exon of FGFR1 gene, 15 individuals presented with one from 10 mutations in the FGFR2 gene, 10 patients harbored one recurrent variant in the FGFR3, whereas 15 patients harbored nine distinct variants in the TWIST1 gene ( Table 2). Karyotyping revealed three heterozygous deletions: one in locus 7q32.3-q35, second in locus 9p, and third in locus 18q21.32-q23 (Supplementary Table S5). The first two were additionally resized using 4 × 180 k Agilent CMA, while the third was by ES.

DISCUSSION
Craniosynostosis represents a highly heterogeneous medical condition whose etiology has not been yet fully elucidated. The results obtained from cohorts screened worldwide showed that the molecular background could be indicated in merely 21%, to propitiously 62% of patients (Roscioli et al., 2013;Paumard-Hernandez et al., 2015;Timberlake et al., 2016;Lee et al., 2018;Topa et al., 2020). Positive genetic testing is mainly achieved among subgroups with syndromic CS. The reported germline mutations are usually classified as point mutations, however, chromosomal aberrations, copy number variations, minor exonic deletions/duplication, or biallelic inheritance were also reported in the medical literature (Timberlake et al., 2016;Goos and Mathijssen, 2019;Yilmaz et al., 2019).
In this study, we have reported 166 individuals affected with different forms of CS in whom we had applied a three-step diagnostic algorithm (Tier 1-3) (Figure 1). To our best knowledge, this is the first large CS patients' screening in which multi-leveled methods, including chromosomal aberrations and CNVs detection, were applied. The proposed approach allowed us to identify an exact genetic cause in around 43% of all CS patients. Since our multi-leveled molecular diagnostic strategy of CS patients is unique and previously unreported, we were unable to directly compare all the results obtained here with previous similar research studies. This is because all reports that have been submitted so far, were aimed at identifying only point mutations via targeted Sanger sequencing, targeted gene panel NGS or ES. Similar to other researchers analyzing these mutations type occurrence, we have shown that causative variants in the FGFR1, FGFR2, FGFR3, TWIST1, and TCF12 genes account for most common cause of CS (Roscioli et al., 2013;Paumard-Hernandez et al., 2015;Lee et al., 2018;Topa et al., 2020). Notably, FGFR1, FGFR2, FGFR3 variants mainly occur in hotspot positions and, along with TWIST1 gene mutations, were analyzed in Tier 1. Undeniably, this step was crucial in the CS testing algorithm and cost-effective compared to other genetic methods (73% of all diagnoses;-χ 2 (1;71) = 14.43; p < 0.001). On the other hand, four novel variants in the TCF12 gene and additional FGFR2, FGFR3 alterations (other than hot-spots) were found using custom gene panel NGS (Table 2). No other gene included in the applied custom gene panel housed more than one pathogenic variant. Interestingly, based on the medical literature, the EFNB1 gene is usually reported as the sixth most commonly affected gene in CS (Paumard-Hernandez et al., 2015;Miller et al., 2017;Lee et al., 2018). However, we postulate that a very characteristic disease, i.e., craniofrontonasal dysplasia (CFND) resulting from pathogenic variants of the EFNB1, should be considered a standalone genetic disorder in which features other than CS guide the proper diagnosis (Bukowska-Olech et al., 2021). Moreover, CFND represents a classic viscerocranium defect, whereas CS is a neurocranium abnormality. Hence, we have excluded all individuals suggestive of CFND and consequently have not reported EFNB1 mutations in this study.
Next, we have revealed that six patients with syndromic CS carried pathogenic variants in genes involved in epigenetic regulations such as ARID1A (P134), KMT2A (P137), KMT2D (P114), NSD1 (P58, P99, P166) ( Table 2), resulting in Coffin-Siris syndrome type 2, Wiedemann-Steiner syndrome, Kabuki syndrome type 1, and Sotos syndrome type 1, respectively. Such Mendelian disorders, i.e., those resulting from disruptions of epigenetic processes, were termed chromatinopathies. They are all characterized by intellectual disability, immune deficiencies, or skeletal anomalies. However, FIGURE 1 | The scheme of the diagnostic algorithm applied in our study regarding 166 patients affected with craniosynostosis. The algorithm was divided into three tiers.
Occasionally, CS has been reported only in Kabuki syndrome and Sotos syndrome thus far (Tatton-Brown et al., 2005;Martínez-Lage et al., 2010;Topa et al., 2017). The second group of Mendelian diseases in which CS has been noted are RASopathies resulting from the dysregulation of RAS/MAPK pathway. Retrospectively, we have described here one patient carrying CNVs in which BRAF gene deletion occurred (P141) (Bukowska-Olech et al., 2020b). Similar to our finding, other researchers have also highlighted the co-occurrence of both CS and RASopathy, resulting from mutations in other components of RAS/MAPK pathway (Kratz et al., 2009;Takenouchi et al., 2014;Ueda et al., 2017). Custom genes panel allowed us to detect novel pathogenic variants-p.Arg132* in the ERF gene (P129), and p.Ser391* in the ZIC1 gene (P131), resulting in Craniosynostosis 4 and Craniosynostosis 6, respectively (Twigg et al., 2013. Both ERF and ZIC1 are newly recognized CS-related genes, however, only a few cases carrying variants in those two have been reported Miller et al., 2017;Glass et al., 2019). In addition, we have evaluated one variant in the ZIC1 gene as VUS p.Ser404Pro (P130) since it was present in the patient's healthy father. However, this alteration was absent from the gnomAD v3.1.2 database (accession date: 3 December 2021). Finally, ES revealed two additional alterations in individuals with syndromic CS-p.Arg1295* in the MN1 (P63), which was recently discovered, and subsequently linked to CS, and p.Arg305Gly in FAM111A (P136), resulting in Kenny-Caffey syndrome type 2, in which CS represents an unseen clinical feature ( Table 2) (Mak et al., 2020). Importantly, we have noted an apparent gap regarding screening for chromosomal aberrations or CNVs among CS patients (Lattanzi et al., 2012;Poot, 2019). To our knowledge, no major CS groups were analyzed via karyotyping or CMA, therefore most research data describing microscopic chromosomal changes or submicroscopic CNVs causative for CS were published as single case reports (Villa et al., 2007;Marques et al., 2015). Besides, only a few chromosomal aberrations and CNVs known thus far represent recurrent changes underlying CS (e.g., deletions in 7p21, 9p22-p24, and 11q23-q24 or duplication in 5q33.3), however none of them was present in our cohort (Reardon et al., 1993;Shiihara et al., 2004;Poot, 2019). All genomic losses or gains reported in this research were not commonly associated with CS (Budisteanu et al., 2010;Dilzell et al., 2015). Hence, we could not recommend additional loci to be screened regarding the cohort of syndromic CS, especially those associated with intellectual disability. Here, karyotyping followed by GTG banding allowed us to detect two intrachromosomal deletions 7q32.3-q35, and 18q21.32-q23, both resized to chr7:131837067-144607071, and chr18:620405559-80247644, respectively. Next, using CMA, we have detected four CNVs including three duplications-1q22-q23.1 (chr1:155961428-157217426), 2p21 (chr2:44990857-45008348), 17p13.3 (chr17:847955-1641601), and one deletion 5q35.3 (chr5:177277901-177283748) from which the largest encompassed 1.3 Mb, whereas the smallest 5.8 kb (Sowińska-Seidler et al., 2018;Bukowska-Olech et al., 2020a). The detailed mutations description following the International System for Human Cytogenomic Nomenclature (ISCN) was listed in Table 3; for a list of genomic mutations' content, see Supplementary Table  S5. Notably, in most chromosomal aberrations or CNVs identified here, CS occurred as an additional phenotype. The Frontiers in Molecular Biosciences | www.frontiersin.org April 2022 | Volume 9 | Article 865494 8 above findings suggest that karyotyping and CMA cannot be replaced by targeted chromosomal testing when applied in a cohort of syndromic CS associated with an intellectual disability or developmental delay.
Regarding the results presented in this study, we would like to point to the possible underestimated co-occurrence of CS and intellectual disability. Our clinical experience suggests that CS may be overlooked when intellectual disability constitutes a primary clinical complaint. Hence, we recommend calvarial sutures' evaluation in patients with intellectual disability. On the other hand, in any case of already detected syndromic CS and intellectual disability, the possible occurrence of clinical features suggestive for either chromatinopathies or RASopathies should also be considered (Kratz et al., 2009;Cao et al., 2017;Zollino et al., 2017;Davis et al., 2019).
Undeniably, the molecular diagnosis of CS should distinguish its isolated or syndromic form, which presence determines the subsequent diagnostic steps. In addition, syndromic CS should be classified as a disorder associated with intellectual disability or a disorder without intellectual disability ( Figure 2). Targeted PCR and Sanger sequencing of FGFR1, FGFR2, FGFR3, TWIST1, and TCF12 genes resulted in the highest diagnostic rate in our cohort of craniosynostosis patients strongly recommend analysing those genes first (isolated CS and syndromic CS without intellectual disability). Because of many advantages of NGS-based methods, including mosaicism detection, screening those genes using targeted genes panel via NGS would be optimal. However, based on our results, we suggest applying a custom genes panel limited to the fewer genes, such as FGFR1, FGFR2, FGFR3, TWIST1, TCF12, ERF, ZIC1, RECQL4, and NSD1 (Siitonen et al., 2009;Sharma et al., 2013;Twigg et al., 2013Twigg and Wilkie, 2015;Miller et al., 2017;Kutkowska-Kaźmierczak et al., 2018;Lee et al., 2018;Bukowska-Olech et al., 2020c;Topa et al., 2020). In the case of syndromic CS and intellectual disability, the genetic investigation should start from chromosomal aberrations or CNVs detection. Lastly, when ES data bioinformatic analysis is performed, genes associated with RASopathies and chromatinopathies should be considered.
The heterogeneity of CS is enormous, resulting from either anatomical variability of the disorder, in which different types and number of sutures can be affected, or epidemiological aspects, including isolated presentation of CS or occurrence of various accompanying symptoms (Kutkowska-Kazmierczak et al., 2018). It has also been shown that CS may be detected in at least 180 different syndromes, which often are rare and in which CS is not a pathognomic feature. Besides, some researchers have also documented or postulated the association between CS and two-locus inheritance. Consequently, molecular causes of the disease seem to be complex in some CS individuals and, as presented in this study, the pathogenic variant or affected gene may be restricted to only one individual Flaherty et al., 2016;Timberlake et al., 2016, Timberlake et al., 2018Wilkie et al., 2007;Timberlake and Persing, 2018). However, in many CS patients, including our cohort, genetic causes of observed phenotypes remain unrevealed. One may suspect the presence of deep-intronic and regulatory variants, polygenic inheritance, or even epigenetic influences. Another explanation may be the technical limitations of currently available diagnostic methods. It has been shown, for example, that the application of long-read sequencing in NGS-based methods may clarify the genetic background in many unresolved cases (Fujimoto et al., 2021;Hiatt et al., 2021;Miller et al., 2021;Rastegar and Yasui, 2021). Considering the above, the next steps that we should consider to implement for diagnosis of our unsolved cases include whole-genome sequencing, RNA-seq, whole-genome bisulfite sequencing, or long-read ES.
To conclude, our research may constitute a significant source of epidemiological information as we have presented precise phenotypic and genetic data derived from 166 consecutive CS patients of Caucasian origin. We yielded a 43% diagnostic success rate using the presented approach, highlighting the high occurrence of pathogenic variants within "classic" CS genes, i.e., FGFR1, FGFR2, FGFR3, TWIST1, and TCF12. Moreover, we have critically summarized the applied diagnostic methods (Figure 1) and proposed the optimized, cost-effective diagnostic algorithm, which may be helpful in a daily diagnostic routine of various CS' types ( Figure 2).

DATA AVAILABILITY STATEMENT
The datasets generated during and/or analyzed during the current study are available from the corresponding authors on reasonable request.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Institutional Review Board of Poznan University of Medical Sciences (no. 742/17). Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.