Whole-exome sequencing identified a novel variant in an Iranian patient affected by pycnodysostosis.

BACKGROUND: Whole-exome sequencing (WES) has emerged as a successful diagnostic tool in molecular genetics laboratories worldwide. In this study, we aimed to find the potential genetic cause of skeletal disease, a heterogeneous disease, revealing the obvious short stature phenotype. In an Iranian family, we used solo-WES in a suspected patient to decipher the potential genetic cause(s).
METHODS: A comprehensive clinical and genotyping examination was applied to suspect the disease of the patient. The solo clinical WES was exploited, and the derived data were filtered according to the standard pipelines. In order to validate the WES finding, the region harboring the candidate variant in the CTSK gene was amplified from genomic DNA and sequenced directly by Sanger sequencing.
RESULTS: Sequence analysis revealed a rare novel nonsense variant, p.(Trp320*); c.905G>A, in the CTSK gene (NM_000396.3). In silico analysis shed light on the contribution of the variant to the pathogenicity of pycnodysostosis. This variant was confirmed by Sanger sequencing and further clinical examinations of the patient confirmed the disease.
CONCLUSION: The present study shows a rare variant of the CTSK gene, which inherited as autosomal recessive, in an Iranian male patient with pycnodysostosis. Taken together, the novel nonsense CTSK variant meets the criteria of being likely pathogenic according to the American College of Medical Genetics and Genomics-the Association for Molecular Pathology (ACMG-AMP) variant interpretation guidelines.

INTRODUCTION

Rare syndromes are disorders affecting a scarce number of individuals in the world (Schieppati, Henter, Daina, & Aperia, 2008). The low number of patients or resources makes it difficult to detect the molecular cause of these conditions (Daneshjoo & Garshasbi, 2018). Instead, increasing knowledge in molecular biology in addition to using cutting‐edge technologies such as next‐generation sequencing (NGS) methods has enabled us to detect genetic defects causing some of the rare syndromes such as pycnodysostosis (PDO). The evident phenotype of PDO affecting the skeletal system is short stature (SS), and initially, based on this feature we selected the patient to detect the possible genetic cause(s).

We diagnosed a patient suspected of PDO as a rare autosomal recessive disease (OMIM 265800) (Figure 1a). Although the first case of this disease was described in 1923 by Sedano, Gorlin, and Anderson (1968), since 1962 to now, less than 200 cases have been detected or reported (Naeem, Sheikh, & Ahmad, 2009). Thus, this disease is categorized as a rare disease showing an obvious SS phenotype. The universal prevalence of PDO is estimated to be 1 to 1.7 per million with equal sex distribution (Xue et al., 2011). Besides SS, the most common features of PDO are increased bone density of long bones, pathological fractures with poor healing, open fontanels, stubby hands and feet with dystrophic nails, and to some extent the obtuse mandibular angle (Xue et al., 2011).

Figure 1

(a) This pedigree is comprised of four generations. The arrow appoints the proband of the family. The genetic status is shown as heterozygote: W/c.905G>A or homozygote: c.905G>A/c.905G>A; in this pedigree, white symbols: unaffected; red symbol: affected; squares: men; circles: females; parallel lines: consanguineous marriage, W: wild‐type allele. (b) Schematic genetic and protein maps of the CTSK gene (NM_000396.3). c.905G>A variant is located in exon 8, which encodes the mature domain of CTSK protein. Similar to other most papain‐like cysteine proteases, CTSK contains 329 amino acids that can be categorized in three distinct sections: a 15‐amino acid preregion, a 99‐amino acid proregion, and a 215‐amino acid mature active enzyme (Toral‐López et al., 2011). The low‐PH environment changes inactive CTSK to the active form by the removal of the N‐terminal proregion. (c) UCSC database used to show the conservation of specific nucleotide (G; highlighted as red) including the variant site in vertebrates, particularly in primates. The amino acid sequence of CTSK colored based on conservation scores by the ConSurf database. Scores ranged from 1 to 9, where a score of 9 represented a highly conserved residue. ConSurf demonstrates evolutionary conservation profiles for proteins of known/unknown structure according to the phylogenetic relations between homologous sequences as well as amino acid's structural and functional importance. (d) The 3D structure of CTSK is shown. The picture was rendered with PyMOL (v.0.99rc6). The original site of Trp302 is emphasized by a highlighted zone and locally zoomed

Various mutations in CTSK gene have been reported in PDO patients which most of these mutations occur in the mature domain of CTSK protein (Xue et al., 2011) (Figure 1b); almost 70% of the mutations have been identified in the mature domain of CTSK, 24.24% in the proregion, and 6.06% in the preregion (Xue et al., 2011). The CTSK gene (NM_000396.3) spans around 12 kb and contains 8 exons (Figure 1b). The encoded protein, CTSK or cathepsin K, is a member of the papain‐like cysteine protease family, and there is a high similarity between this protein and other cathepsins such as S and L (Gelb, Shi, Chapman, & Desnick, 1996). CTSK is highly expressed in osteoclasts (Xue et al., 2011), and its mRNA is detectable in macrophages and bone marrow‐derived dendritic cells (Asagiri et al., 2008; Honey & Rudensky, 2003). Additionally, cathepsin K plays a vital role in osteoclast‐mediated bone resorption by degenerating the bone matrix proteins, such as type I collagen, osteopontin, and osteonectin (Hou et al., 1999; Mujawar et al., 2009). In a nutshell, impairment of CTSK‐mediated osteoclast apoptosis/senescence may also be responsible for the higher bone density, which is a hallmark feature in patients affected by PDO (Chen et al., 2007).

In the present study, the main object was to resolve the genetic diagnosis of SS phenotype in an Iranian male patient with normal parents. By applying whole‐exome sequencing (WES), a novel variant, p.(Trp320*), in the CTSK gene was identified. The physical assessments and medical evaluations confirmed the PDO due to the novel identified variant in this family. According to the American College of Medical Genetics and Genomics‐the Association for Molecular Pathology (ACMG‐AMP) variant interpretation guidelines (Biesecker & Harrison, 2018; Oza et al., 2018), this variant could be classified as a likely pathogenic variant, albeit functional studies to further confirm the pathogenicity of the variant in appropriate animal models will be recommended.

MATERIALS AND METHODS

Editorial policies and ethical considerations

The study protocol was approved by the local medical ethics committee of Tarbiat Modares University, Tehran, Iran. A written informed consent was obtained from all subjects before applying enrollment in this study. The family also provided written informed consent for publication of their pertinent images included in this paper. All of the patient's clinical information and the medical histories were collected at the Dr. Azimy Genetic Center, Arak, Iran.

Subject and clinical investigations

We enrolled 3 members of the family in our study (Figure 1a). The consanguineous family, originating from Arak province of Iran, was suspected to pycnodysostosis. The family was ascertained for the present study, and the affected subject underwent meticulous medical records including a comprehensive physical examination, bone density testing, and radiography tests. Some of the important clinical indices used in this study are summarized in Table 1.

Table 1

Clinical features of the patient in the family

Features	IV.1
Age at initial assessment/age at molecular assessment	17/20
Gender	Male
Father age at conception	31
Weight at birth	Below 5th percentile (3.1 ± 0.01 kg)
Height at birth	Below 5th percentile (45.2 ± 0.10 cm)
Head circumference at birth	36.3 ± 0.10 cm
Height at assessment	137 ± 0.10 cm
Weight at assessment	34 ± 0.01 kg
Sexual maturity rating	Stage IV
Short stature (<150 cm)	Observed
Increase in bone density	Observed
Open fontanels and sutures	Open sutures of anterior fontanel and closed posterior fontanel
Frontal and parietal bossing	Observed
Fractures	Developed easy fractures for four times
Obtuse mandibular angle	Observed
Short fingers and hypodontia	Observed
Stubby hands and feet with osteolysis of the distal phalanges	Observed
Prominent eyes with bluish sclerae	Observed
Hypoplastic maxilla	Observed
Grooved palate	Observed
Scoliosis	Scoliosis was evident in the thoracic region
Dysplastic nails	Dysplastic nails were evident in this patient
Nonpneumatized paranasal sinuses	Not applicable
Prominent nose	Observed
Macrocephaly	Observed
Asymmetric skull	Observed

DNA extraction

To apply genetic tests, around 10 ml of blood samples was taken from the patient and his parents, and then, genomic DNA was isolated from the samples by Roche DNA Extraction Kit (Cat. No. 11814770001, Roche Life Science). Thence, DNA concentrations were measured by Thermo Scientific™ Nanodrop 2000 (Thermo Fisher Scientific).

Whole‐exome sequencing

About 1 μg of genomic DNA sample of the patient (IV.I) was subjected to high‐throughput sequencing. RNA capture baits against approximately 60 Mb of the Human Exome (targeting > 99% of regions in CCDS, RefSeq, and Gencode databases) were used to enrich regions of interest from fragmented genomic DNA with Agilent's SureSelect Human All Exon V6 kit. The generated library was sequenced on an Illumina platform to obtain an average coverage depth of approximately 100×. Typically, around 97% of the targeted bases were covered >10×.

Bioinformatics analysis

An end‐to‐end in‐house bioinformatics pipeline including base calling, alignment of reads to GRCh37/hg19 genome assembly, primary filtering of low‐quality reads and probable artifacts, and subsequent annotation of variants was applied. The reads were aligned to the reference genome (hg19/NCBI37.1) with SNP & Variation Suite version 8.0 (SVS v8.0) and DNASTAR Lasergene12 (DNASTAR Inc.). Small indel detection was used with the Unified Genotyper tool from GATK tools in Galaxy online database (http://www.usegalaxy.org). The missense, nonsense, silent, and indel variant rates were estimated by Galaxy online tool and finally were confirmed by Ivariantguide® (Chaudhry & Tainsky, 2019).

Regarding the pedigree and inheritance mode, we assumed that the variant(s) might be inherited as homozygous manner and we hereby excluded the heterozygous variants, and then, several filtering steps were applied to prioritize all variants: (a) Variants in dbSNP132 (https://www.ncbi.nlm.nih.gov/projects/SNP), 1000 Genomes Project (Siva, 2008), and GnomAD (Karczewski & Francioli, 2017) with allele frequencies more than 1% were excluded. (b) The rest of the variants underwent further exclusion in the Exome Sequencing Project (ESP) (http://evs.gs.washington.edu/EVS) and Exome Aggregation Consortium (ExAC) database. (c) The intragenic, intronic, untranslated regions (UTRs), and synonymous variants were excluded from later analysis. (d) The SIFT (Ng & Henikoff, 2003), Provean (Choi & Chan, 2015), and MutationTaster (Schwarz, Rödelsperger, Schuelke, & Seelow, 2010) were used to predict the pathogenicity/prioritization of candidate variants (Table 2). To be on the safe side, the filtering was applied with assuming the heterozygous mode of inheritance, but this did not lead to candidate variant(s) in this family.

Table 2

Several online databases that used to confirm the pathogenicity of the variant in CTSK

Variant	Gene	Zygosity			MutationTaster	EXAC	SIFT	1K Genome	Iranome	PROVEAN	ENTPRISE‐X
		Index	Mother	Father
p.(Trp320*)	CTSK	Hom.	Het.	Het.	Disease causing	N.R	Damaging	N.R	N.R	Deleterious	Deleterious

Abbreviations: Het, Heterozygote; Hom, Homozygote; N.R, Not Reported.

To narrow down the list of candidate variants as much as possible, we used Phenolyzer (Yang, Robinson, & Wang, 2015), Face2Gene (https://www.face2gene.com/), and Varcards (Li et al., 2017). All suspected pathogenic variants were double‐checked in HGMD® (Stenson et al., 2003) and ClinVar (Landrum et al., 2015).

ConSurf server (Glaser et al., 2003) and UCSC database (Karolchik et al., 2003) were applied to provide an evolutionary conservation profile for CTSK protein and DNA sequence, respectively, to better predict the potential disrupting effects of the variant (Figure 1c and d). For further consideration, the frequency of the variant was checked out on Iranome as the local database (Fattahi et al., 2019) (Figure 2a). All information related to the in silico predictions such as allele frequency, cosegregation results, and pathological predictions is summarized in Table 2.

Figure 2

(a) Schematic representation of filtering strategies exploited in this research. For more investigation, the filtering steps evaluated by regard to this fact that the disease could be engendered by autosomal dominant; however, we could detect no relevant variant according to this supposition. (b) Sequence chromatogram shows a homozygous state of the nucleotide sequence of c.905G>A. (c) Radiography shows the large head circumstance, obtuse mandibular angle, and scoliosis, which was evident in the patient. (d) Short fingers revealed by X‐ray and the dysplastic flat nails were evident in the patient (For detailed, please refer to Figure S1)

Segregation analysis

Sanger sequencing in forward and reverse directions was performed to validate the candidate variant, and then, segregation analysis was performed in the family. The primers were designed by Primer3.0 (http://bioinfo.ut.ee/primer3-0.4.0) Web‐based server (Table S1). The lack of SNPs in the genomic region corresponding to the 3′ ends of primers was inspected by looking through the dbSNP database. The specificity of primers was checked by the in silico‐PCR tool in UCSC genome browser (https://genome.ucsc.edu/cgi-bin/hgPcr) and Primer blast on NCBI genome browser (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). Finally, polymerase chain reaction (PCR) was performed in standard conditions, and then, PCR products were sequenced by ABI 730XL, using the conventional capillary system. Chromatograms were analyzed by GenomeCompiler online tool (http://www.genomecompiler.com) and Mutation Surveyor v.3.24 (Softgenetics) to identify the alternation. Variants were annotated based on the standards of the Human Genome Variation Society (HGVS) nomenclature (den Dunnen et al., 2016).

Three‐dimensional structure modeling

The identified variant, p.(Trp320*), is located in the mature domain of cathepsin K (Figure 1e). The protein families and domains were analyzed using ScanProsite (Gattiker, Gasteiger, & Bairoch, 2002), and sequence alignments were recruited by using ClustalW (http://www.ebi.ac.uk/clustalw). The protein structures and possible effects of the novel variant on cathepsin K were analyzed by PyMOL (DeLano, 2002) after building the PDB structure file based on the Phyre2 (Kelley, Mezulis, Yates, Wass, & Sternberg, 2015) and SWISS‐PROT (Bairoch & Boeckmann, 1992). The RAMPAGE online tool (https://www.mordred.bioc.cam.ac.uk/~rapper/rampage.php) was applied to check out the detailed residue‐by‐residue stereochemical quality on the basis of a Ramachandran plot to know whether the modeled structure was acceptable.

RESULTS

Clinical findings

The patient, 18‐year‐old male, showed typical clinical features of the hereditary PDO, for example, abnormal skeletal system, scoliosis, SS with skeletal dysplasia, increased density of bones, open sutures of anterior fontanel and closed posterior fontanel, history of easy fractures for four times, frontal bossing, micrognathia, asymmetric skull and macrocephaly, short fingers with dysplastic flat nails, midface retrusion, and prominent nose. In addition, dental abnormalities including severe crowding, poor oral hygiene, periodontal problems, delayed exfoliation of primary teeth, an eruption of permanent teeth, enamel hypoplasia, obliteration of pulp chambers, and hypercementosis were evident. The patient's radiographs exhibited an obtuse mandibular angle, the general increase in bone density, and open fontanels and sutures (Figure 2c).

In the assessment period, Tanner Stages, also known as Sexual Maturity Rating, was examined as stage IV. The history of short stature and growth retardation noted at 3 years of age, while at birth time weight (3.1 ± 0.01 kg; reference in an Iranian population: 3.26 ± 0.60 kg) and height (45.2 ± 0.10 cm; reference in an Iranian population: 51.9 ± 3.59 cm) were below 5th percentile in comparison with the normal values (Heydari, Emamghoreishi, & Amini, 2009) and head circumference (36.3 ± 0.1 cm; reference in an Iranian population: 35.13 ± 1.45 cm (Esmaeili & Esmaeili, 2015)), was above 95th percentile. Further bone evaluations showed increased bone density with the preservation of growth plate and sign of old fractures. Other important detected clinical features are summarized in Table 1. In addition to obvious dysplastic nails, radiographic examination of the index case verified the increased bone density and short fingers in comparison with the average values (Figure 2d).

Regarding the pedigree provided, the parents were consanguineous. Further laboratory analysis affirmed normal values for leukocyte and thrombocyte counts, plasma phosphate, and alkaline phosphatase. Additionally, no obvious phenotypic abnormality was observed in the parents.

Genetic analysis

To elucidate the underlying genetic cause(s), genomic DNA was obtained from the patient and analyzed by solo clinical WES, and finally, p.(Trp320*) variant was confirmed by Sanger sequencing (Figure 2b).

In the patient, the detected SNVs and deletion/insertions were analyzed by several filtering steps. In total, 91,395 variants were found by WES in the proband after alignment and SNV calling. By applying several exclusion processes such as base quality filtering, and frequency in dbSNP132, 1000 Genomes Project, ESP, GnomAD, and ExAc databases, 51 variants were identified. These variants then prioritized based on the patient's phenotypes.

Eventually, regarding the phenotypes by utilizing CentoMD (Trujillano), ClinVar, phenolyzer (Yang et al., 2015), Face2Gene, Human phenotype ontology (Robinson et al., 2008), and Varcards, only one relevant variant was identified. Consequently, the samples from the available members of the family were subjected to Sanger sequencing and segregation of the candidate nonsense variant of the CTSK gene was confirmed (Figure 2b). Finally, we classified the variant based on ACMG‐AMP guidelines (http://wintervar.wglab.org) into a likely pathogenic variant.

To caveat, besides considering autosomal recessive inheritance mode for this family, we considered autosomal dominant inheritance, but this analysis did not result in any candidate variant. All details of filtering steps are accessible in Table 3, and all statistics regarding the variants are given in Table S2.

Table 3

The number of variants identified in patients through WES

Patient	IV.I
Total variants	91,395
Variants after base quality filtering	77,909
Homozygous variants	31,383
Nonsynonymous/indel/splice site variants	8,928
Novel variants (dbSNP132/1000GP queried)	51
Genes with plausible disease association	6
Phenotype analysis	1

In silico predictions

We used RAMPAGE online tool to show the detailed residue‐by‐residue stereochemical quality of CTSK in accordance with the Ramachandran plot. The structural model of CTSK indicated almost 98.1% of residues in the most favored regions, around 0.7% of residues in allowed regions, and only 0.2% of residues in the outlier regions, which suggested that the modeled structure of this protein was acceptable.

Various in silico predictor databases such as SIFT, PolyPhen‐2, MutationTaster, Provean, and ENTPRISE‐X (Zhou, Gao, & Skolnick, 2018) were used to evaluate the possible pathogenicity of the variant. All detailed results are described in Table 2.

DISCUSSION

Bone development is a complex process in which a balance between bone formation and resorption is delicately maintained (Razmara et al., 2019). In fact, osteoblasts and also osteoclasts have an important hand in this process, and any perturbation influencing this procedure causes multitudinous genetic bone diseases ranging from osteoporosis to osteopetrosis (Tanaka, Nakayamada, & Okada, 2005). To date, many genes have been identified related to skeletal abnormalities (Table S3), although the molecular and cellular basis of a great number of skeletal genetic disorders is still uncertain and this knowledge gap can be filled by using advanced technologies such as NGS, which has importantly expedited the detection of responsible genetic causes. Among NGS techniques, WES has been developed into a robust and cost‐effective tool to identify the genetic cause in rare diseases (Razmara & Garshasbi, 2018). Using WES, we succeeded to identify a novel variant in the CTSK gene in an Iranian male patient presenting SS phenotype. We also provided evidence supporting the causative role of homozygous p.(Trp320*) in index case using in silico predictions. Nevertheless, the messy reality of unknown significant variants is not going away anytime soon, much as clinicians, geneticists, and patients might wish it would; in this study by applying different stages of filtering, 51 unknown significant variants were detected in the patient from which 6 genes were plausible showing an association with reported skeletal diseases. To narrow down as much as possible, we used the clinical data of the patient in various databases, for example, Face2gene (Table 3). This helped us to decrease the number of candidate variants to only one variant, thence confirmed by Sanger sequencing. The c.905G>A variant in the CTSK gene segregated with the condition in the family.

Cathepsin K deficiency does not influence on the function of osteoclast‐mediated extracellular acidification (Chen et al., 2007), whereas CTSK mutations have been detected to halter the ability of osteoclasts to degrade collagen rather than demineralize the extracellular matrix. Furthermore, cathepsin K may also function as a potential regulator of apoptosis and senescence, controlling osteoclast numbers in vivo (Chen et al., 2007). Studies on the Ctsk −/− mice have revealed that the number of chondroclasts and osteoclasts is increased in joint tissues of the mice, whereas osteoclast numbers are increased in the mouse bones (Soki et al., 2018). Similarly, some scrapes of evidence showed that the deletion of CTSK in osteoclasts enhances bone formation in vivo (Lotinun et al., 2013). These data can justify the increased bone density in the patients affected by pycnodysostosis.

The clinical presentation of the patient in the family was originally described as SS and solo‐WES detected a novel variant in the CTSK gene. The hallmark features of this disease such as increased bone density, frontal bossing, micrognathia, prominent nose, short stature, delayed abnormal tooth eruption, and fragility fractures were observed in the index case. Indeed, one of the two original descriptions of PDO in 1962 named it an osteopetrosis variant (Pangrazio et al., 2014). However, PDO is usually a progressive but relatively benign condition. The applied clinical assessments confirmed that the patient is affected by pycnodysostosis, not osteopetrosis. Evaluation of the radiographs available for the patient showed the obtuse mandibular angle on a craniolateral view and absence of acroosteolysis of the hands (Figure 2c). Moreover, it seems that haploinsufficiency, a dominant phenotype in diploid organisms that are heterozygous for a loss‐of‐function allele, is not responsible for the impairment of CTSK, because none of the parents showed pertinent clinical symptoms. Thus, not only did this help us to exclude the osteopetrosis but also it provided cogent evidence of pycnodysostosis in the index case.

The p.(Trp320*) variant may lead to the production of a truncated protein with loss of downstream functional domain in the cell; as a result, cells should forestall this process by nonsense‐mediated decay response (NMD response), which is increasingly appreciated as one of the central mechanisms of RNA surveillance, with a great role in the physiological control of gene expression. Chen et al. showed that a deficiency of cathepsin K prevents inflammation and bone erosion in rheumatoid arthritis and periodontitis and reveals its shared osteoimmune role (Hao et al., 2015). Based on the identified variant, we propose NMD or loss of downstream functional domain as judicial mechanisms of pathogenesis of p.(Trp320*) variant in the patient, but more studies need to unearth the exact molecular mechanisms that contribute to the pathogenicity.

Additionally, it has been shown that cathepsin K plays a role in hormone activation or degradation (Tepel, Bromme, Herzog, & Brix, 2000), glucose metabolism (Yang et al., 2008), and pathogenesis of obesity (Xiao et al., 2006). Thus, the identification of novel functional variants can broaden the horizons toward understanding the mechanisms in which cathepsin K plays roles.

To conclude, our results indicate that the novel nonsense variant, c.905G>A; p.(Trp320*), in the CTSK gene (NM_000396.3) might be the genetic cause of pycnodysostosis. The putative variant meets the criteria of being pathogenic, but we firmly advise applying functional analysis with animal models to inspect the distinctive pathological roles of this variant.

CONCLUSION

In summary, we report the CTSK variant p.(Trp320*), which is associated with PDO in a male Iranian patient. We hope that this identification may also yield new insights into the mechanisms of human bone disease, especially in relation to metabolic alterations occurring during this process. We also suggest doing functional analysis by applying appropriate animal models to decipher the mechanism of pathogenesis of the variant prior to using in genetic counseling.

CONFLICT OF INTEREST

The authors stated no conflict of interest.

AUTHORS' CONTRIBUTION

Conceived and designed the experiments: M.G, E.R, M.A.D. Conducted the experiments: M.G, M.A.D, and E.R. Analyzed and interpreted the data: E.R, M.A.D, H.A, and E.E.G. Contributed reagents/materials/analysis tools: E.R. Wrote the paper: E.R, M.G, A.R.B, M.GH, and M.D. Designed the figures: E.R. All authors read and approved the final manuscript.

Click here for additional data file.

Click here for additional data file.