Whole Exome Sequencing Analysis in Fetal Skeletal Dysplasia Detected by Ultrasonography: An Analysis of 38 Cases

Background: Skeletal dysplasias (SDs) are a heterogeneous group of genetic disorders that primarily affect bone and cartilage. This study aims to identify the genetic causes for fetal SDs, and evaluates the diagnostic yield of prenatal whole-exome sequencing (WES) for this disorder. Methods: WES was performed on 38 fetuses with sonographically identified SDs and normal results of karyotype and single nucleotide polymorphism (SNP) analysis. Candidate variants were selected by bioinformatics analysis, and verified by Sanger sequencing. Results: WES revealed pathogenic or likely pathogenic variants associated with SDs in 65.79% (25/38) of fetuses, variants of uncertain significance (VUS) in SDs-related genes in 10.53% (4/38) cases, and incidental findings in 31.58% (12/38) fetuses. The SDs-associated variants identified in the present study affected 10 genes, and 35.71% (10/28) of the variants were novel. Conclusion: WES has a high diagnostic rate for prenatal SDs, which improves pregnancy management, prenatal counseling and recurrence risk assessment for future pregnancies. The newly identified variants expanded mutation spectrum of this disorder.


INTRODUCTION
Skeletal dysplasias (SDs) are a heterogeneous group of disorders that affect bone and cartilage and are characterized by abnormal shape, growth, and integrity of the skeleton. Its prevalence is estimated to be about 1 of 4,000-5,000 live births (Orioli et al., 1986;Stoll et al., 1989), and account for nearly 10% of fetal structural malformations detected by ultrasonography (Tang et al., 2020). Currently, SDs comprise 461 types of diseases that are categorized into 42 groups, and pathogenic variants in 437 different genes have been linked to 92% (425/461) of these diseases (Mortier et al., 2019). In addition, pathogenic variants in an individual gene can cause different phenotypes, and similar phenotypes can be resulted from variants in different genes. The high genetic and phenotypic diversity make prenatal diagnosis of SDs by ultrasound or magnetic resonance imaging challenging (Schramm et al., 2009;Geister and Camper, 2015), while molecular diagnosis offers the possibility of identifying the specific entity underlying SDs.
The conventional molecular diagnosis strategies include karyotyping and chromosomal microarray analysis (CMA). Karyotyping has a diagnostic yield of around 30% in fetuses with sonographically identified structural birth defects (Zhang et al., 2017), and CMA provides an additional diagnostic yield of 4-6% in fetuses with an ultrasound anomaly and normal karyotype (Wapner et al., 2012). Thus, using these molecular diagnosis methods, less than half of fetuses with structural anomalies can be diagnosed. In addition, SDs are mainly monogenic disorders (Geister and Camper, 2015;Mortier et al., 2019), and therefore molecular diagnosis methods with high resolution is need for their diagnosis. Next-generation sequencing (NGS) provides the deep, high-throughput, and in-parallel DNA sequencing. Wholeexome sequencing (WES) is a NGS method that sequences exons that contain > 85% of the genetic variants associated with human disease phenotypes (Rabbani et al., 2014). Two recent large-scale prospective studies in fetus with structural anomalies revealed that WES provided an additional diagnostic yield of 8-10% in fetuses with an ultrasound anomaly and normal karyotype and CMA results, and the detection rate is strongly correlated with the number of fetal anomalies (Lord et al., 2019;Petrovski et al., 2019). Lord et al. demonstrated that WES had a diagnosed yield of 15·4%(10/65) in fetuses with skeletal anomalies and normal karyotype and CMA results (Lord et al., 2019). In a cohort of 28 chromosomally normal fetuses with SD, Liu et al. (2019) revealed that 75% (21/28) of cases were detected with mutations in genes related to skeletal diseases by WES.
In the present study, WES was performed to identify genetic causes for 38 fetuses with sonographically identified SDs and normal results of karyotype and single nucleotide polymorphism (SNP) array analysis, and to evaluate the diagnostic rate of WES for prenatal SDs.

Subjects and Sample Collection
The study was approved by the ethics committee of Hunan Provincial Maternal and Child Health Care Hospital. All parents agreed to participate in the study and provided signed informed consent.
Thirty-eight fetuses diagnosed as SDs by prenatal ultrasonography at our center from 2016 to 2019 were enrolled in this study. The sonographic criteria incudes: (1) the presence of short limb deformities in which fetal femur length (FL) < −2 SD or FL below the 5th centile of our reference ranges at midtrimester ultrasonography), and/or (2) the presence of other skeletal anomalies including various deformities, finger/toe deformities, missing fingers/toes, and/or absence of upper/lower limbs. The overall study workflow was shown in Figure 1. The clinical data and the results of prenatal ultrasound examination are shown in Table 1. The fetal amniotic fluid, umbilical cord blood or skin tissue of aborted fetuses, and peripheral blood of family members were collected (Table 1).

Library Construction and Sequencing
Genomic DNA of the 38 fetuses and 7 pairs of parents were extracted using a Qiagen DNA Blood Midi/Mini kit (Qiagen GmbH, Hilden, Germany) as previously described (Peng et al., 2020). DNA was sheared into fragments of around 150 bp and blunt-ended followed by addition with deoxyadenosine at the 3 ends. The DNA library was generated by ligating adaptors to the ends of double DNA strands that enable sequencing. The library was amplified by PCR and then hybridized to a pool of biotinylated oligo probes specific for exons. Streptavidin magnetic beads were used to capture DNA-probes hybrids. The hybrid products were eluted and collected, and were subsequently amplificated by PCR. The libraries were tested for enrichment by qPCR, and size distribution and concentration of the libraries were determined using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, United States). For sequencing, the genomic DNAs were sequenced using NovaSeq 6000 platform (Illumina, San Diego, CA, United States) with 150 bp pair-end reads. Raw image files were processed using CASAVA (v1.82) for base calling and generating raw data.

Data Analysis
Adaptor sequences of the raw data were trimmed at the tails of reads using Cutadapt (v1.15) (Martin M. Cut adapt removes adapter sequences from high-throughput sequencing reads. EMBnet 17:10-12[J]), and then aligned to the human reference genome (hg19/GRCh37) with using Burrows-Wheeler FIGURE 1 | The overall experimental workflow in the present study. ACMG, American College of Medical Genetics and Genomics; Verita trekker R and Enliven R represent patented software of Berry Genomics; CNV, copy number variation; HGMD, Human Gene Mutation Database; InDel, insertion/deletion; Mt, mitochondrial variation; QC, quality control; ROH, runs of homozygosity; SNP, single nucleotide polymorphism; WES, whole-exome sequencing. The formula: Z-score = (XGA -MGA)/SDGA, where XGA is data from other populations at a known gestational age, MGA is the mean value for our population calculated from the reference equations at this gestational age, and SDGA is the SD associated with the mean value at the same gestational age from our population. The fetal biometric measurements of other populations at each gestation were expressed as Z-scores calculated with our reference equations.   (Wang et al., 2010) and the Enliven R Variants Annotation Interpretation System authorized by Berry Genomics. Annotation databases mainly included gnomAD 3 , the 1000 Genome Project 4 , Berry big data population database, dbSNP 5 , OMIM 6 , ClinVar 7 , HGMD 8 , and HPO 9 . The mean depth of coverage for whole exomes by WES was 100× and a mean of 97.5% of bases was covered by at least 20 reads.

Data Filtering and Diagnostic Interpretation of the Filtered Variants
Variants that have any of the following conditions were excluded, including: (1) sequencing depth is ≤ 10; (2) maximum population frequency was larger than 0.01 when referring to dbSNP build 142, the 1,000 Genomes Project and the ExAC Browser; (3) not in coding or exon-intron junction region; and (4) predicted as benign by SIFT, FATHMM, MutationAssessor, CADD, and PolyPhen-2. The remaining variants were identified as potential SNVs. The variants identified were then classified into to five categories according to the American College of Medical Genetics and Genomics (ACMG) recommendations of 2015 (Richards et al., 2015), including pathogenic, likely pathogenic, benign, likely benign, and variants of uncertain significances (VUSs). The ACMG guidelines define 28 criteria that are grouped into seven categories of evidence including pathogenic very strong (PVS), pathogenic strong (PS), pathogenic moderate (PM), pathogenic supporting (PP), benign stand-alone (BA), benign strong (BS), benign supporting (BA). Each of the seven categories includes one or multiple criteria (for example, PS has four criteria PS1-PS4), and each criterion is defined based on a variety of information including: (1) the type of variant and population frequency, (2) the patient's overall clinical presentation and its concordance with the phenotype described in the literature, (3) the family history, the predicted mode of transmission, and the familial segregation, (4) the prevalence of the disease compared with the frequency of the variant in the general population, (5) the predicted effect of the variant by in silico predictive algorithms, and (6) the functional data available in the literature.
By combining these criteria, a variant is classified into one of the five categories.

Verification of Gene Mutations by Sanger Sequencing
The pathogenic or likely pathogenic mutations detected by WES were confirmed by Sanger sequencing.

Genetic Variations
The karyotyping results were normal for the 31 fetuses (data not shown). WES was successfully performed on 31 fetus samples and 7 father-mother-proband samples (fetus 1, 12, 20, 25, 31, 35, 36 and their parents). The results were shown in Figure 2.
A total of 774 potential SNVs were initially identified, including 626 from the fetus only group and 148 from the fetus-mother-father trio group. Among them, 34 variants from the fetal proband-only group and 10 variants from the fetusmother-father trio group were evaluated to be of potential clinical significance. Among the 44 variants, 22 pathogenic variants and 7 likely pathogenic variants from 25 fetuses were evaluated to be directly associated with SDs (Table 1 and Figure 2), thereby yielding a positive detection rate of 65.79% (25/38). 4 VUS in SDs-related genes were detected in 4 fetuses ( Table 1 and Figure 2), yielding an VUS detection rate of 10.53% (4/38). Among the 29 SDs-related pathological/likely pathogenic variants, 20 were de novo variants that followed an autosomal dominant inheritance pattern, including 18 missense mutations, 1 frame-shift mutation, and 1 splice-site mutation, and 9 variants were functional homozygous or compound heterozygous mutations in accordance with their recessive patterns of inheritance, including 4 missense mutations, 1 nonsense mutation, 3 were frameshift mutations, and 1 splice site mutation (Table 1).
When categorizing the fetal samples into short-limb subgroup and non-short limb group, a definitive diagnosis was made in 70.0% (21/30) cases with short-limb and in 62.5% (5/8) cases without short-limb.

DISCUSSION
Skeletal dysplasias are a heterogeneous group of genetic disorders that account for nearly 10% of the fetal structural malformations (Tang et al., 2020). In this study, WES was performed on 38 fetuses with sonographically identified SDs and normal results of karyotype and SNP array analysis. Pathogenic/likely pathogenic variants were identified in 65.79% (25/38) of the fetuses, which indicates a high diagnostic yield of WES for prenatal molecular diagnoses of SDs. 10 novel variants identified in this study expand mutation spectrum of SDs and contribute to the genetic diagnosis and counseling of this disorder. This study also identified pathogenic variants during the WES data analysis that are beyond the scope of the SD for which the fetuses were prescribed WES tests, yielding an incidental detection rate of 31.58% (12/38). Current ACMG guidelines recommend the reporting of incidental findings in clinical exome and genome sequencing, specifically pathogenic variants in 59 medically actionable genes (Kalia et al., 2017), and therefore this high rate of incidental findings will present greater ethical challenges during genetic counseling.
Pathogenic variants from 437 different genes have been associated with SDs (Mortier et al., 2019). In the present study, 28 variants affecting 10 genes including FGFR3, COL1A1, COL1A2, COL2A1, TRIP11, SOX9, LBR, IFT172, FIG4, and DYNC2H15 were found to be associated with the SDs (Table 1). Among the 25 fetuses of SDs with genetic diagnosis, 52% (13/25) of fetuses contained variants in FGFR3, 16% (4/25) of fetuses contained variants in COL2A1, and 8% (2/25) fetuses contained variants in COL1A1.The other 4 fetuses carried variants in one of the 6 genes including TRIP11, SOX9, LBR, IFT172, FIG4, and DYNC2H1. These findings are in line with the high genetic variability of SDs, and also highlights that pathogenic variants in FGFR3 and collagen genes are the most common generic lesion for SDs.
Fibroblast growth factor receptor 3 (FGFR3) is one of four distinct membrane-spanning tyrosine kinases that serve as highaffinity receptors for a number of fibroblast growth factors and plays essential roles in skeletal development (Ornitz and Marie, 2015). Currently, variants in FGFR3 have been associated with at least 10 different skeletal disorders (Mortier et al., 2019), and the distribution of the known disease-associated FGFR3 germline mutations suggests there is no preferential type or location of mutations for a distinct disorder. In the present study, fetus 1 and 2 carried an identical variant c.742C > T (p.Arg248Cys) in the extracellular region, fetus 3 and 4 carried a variant c.746C > G (p.Ser249Cys) in the extracellular region, and fetus 8 carried mutation c.1118A > G (p.Tyr373Cys) in the transmembrane domain. These three mutation have been associated with thanatophoric dysplasia, type I (MIM_187600) (Brodie et al., 1999).The variant c.1144G > A (p.Gly382Arg) in fetus 5-7 and the variant c.1138G > A (p.Gly380Arg) in fetus 9-12 localize in the transmembrane domain, all of which are associated with achondroplasia (MIM_100800) (Hafner et al., 2006;Natacci et al., 2008). Fetus 13 inherited a heterozygous mutation c.1620C > G (p.Asn540Lys) in the TK-1 domain of FGFR3, which is a common cause of hypochondroplasia (OMIM_146000) (Mortier et al., 2000). Notably, a novel heterozygous COMP mutation c.1255-5C > T was also detected in fetus 13 and the mother, but not in the healthy family members (Figure 3). Variants in COMP has been linked to SDs (Mabuchi et al., 2003). Analysis through the ClinVar database indicates the COMP mutation c.1255-5C > T is benign, but whether this variant also contributed to the clinical phenotypes in the family of fetes 13 awaits further investigation.
Mutations in the type I collagen coding genes (COL1A1 and COL1A2) that affect collagen quantity or structure count to approximately 85% of osteogenesis imperfecta (OI) case (Marini et al., 2017;Tournis and Dede, 2018). Mutations that reduce the synthesis of normal type I collagen (quantitative defect) are usually associated with the milder OI, while mutations that alter the protein structure of collagen molecules (structural defect) mostly caused by substitutions of glycine by different amino acids in the helical domain of type 1 collagen are usually linked with more severe phenotypes (Marini et al., 2017). In the variants identified in the present study including COL1A1 variant c.1192G > T (p. Gly398Cys) and c.1373G > A (p.Gly458Glu), and COL1A2 variant c.2422G > T (p. Gly808Glu), a glycine in the helical domain of type 1 collagen was substituted by other amino acids. Thus, these variants may damage type 1 collagen functions and therefore cause SDs. COL2A1 encodes type II collagen that plays a role in the regulation of intramembranous and endochondral ossification. Heterozygous COL2A1 mutations are usually associates with a spectrum of dwarfism and skeletal malformation diseases (Zhang et al., 2020). In the present study, 2 novel likely pathogenic mutations c.1358G > T (p.Gly453Val) and c.3571G > A (p.Gly1191Arg), and 1 known pathogenic mutation c.1419 + 5G > T were identified in COL2A1 from fetus 17 to 19 (Table 1). Collectively, these findings support high genetic variability and also highlight the importance of molecular diagnosis for fetal SDs.

CONCLUSION
This study has proved that prenatal WES has a high diagnostic rate for SDs, which improves the clinical management of pregnancies and better inform family planning efforts. This study has also identified novel pathogenic variants for SDs, which broadens the mutation spectrum for this disorder and contributes to clinical consultation and subsequent pregnancy examination. As previously reported, VUS and incidental findings in prenatal genetic diagnosis by WES remains as counseling challenges for prenatal diagnosis for structural anomalies, including SDs. Moreover, large-scale prospective studies in fetuses with sonographically identified SDs will provide further information on the feasibility and potential impact of WES on prenatal counseling and pregnancy management.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below: National Center for Biotechnology Information (NCBI) BioProject, https://www.ncbi.nlm.nih.gov/ bioproject/, PRJNA747169.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the ethics committee of the Hunan Provincial Maternal and Child Health Care Hospital. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin. Written informed consent was obtained from the minor(s)' legal guardian/next of kin for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
HW, YP, and CT conceived the study. XH analyzed the results of prenatal ultrasound testing. SY, JP, JL, and JH analyzed the results of karyotyping and SNP array analyses. YP, XS, and CT analyzed the results of WES and prepared the manuscript. All authors contributed to the article and approved the submitted version.