Novel Pathogenetic Variants in PTHLH and TRPS1 Genes Causing Syndromic Brachydactyly.

Skeletal disorders, including both isolated and syndromic brachydactyly type E, derive from genetic defects affecting the fine tuning of the network of pathways involved in skeletogenesis and growth-plate development. Alterations of different genes of this network may result in overlapping phenotypes, as exemplified by disorders due to the impairment of the parathyroid hormone/parathyroid hormone-related protein pathway, and obtaining a correct diagnosis is sometimes challenging without a genetic confirmation. Five patients with Albright's hereditary osteodystrophy (AHO)-like skeletal malformations without a clear clinical diagnosis were analyzed by whole-exome sequencing (WES) and novel potentially pathogenic variants in parathyroid hormone like hormone (PTHLH) (BDE with short stature [BDE2]) and TRPS1 (tricho-rhino-phalangeal syndrome [TRPS]) were discovered. The pathogenic impact of these variants was confirmed by in vitro functional studies. This study expands the spectrum of genetic defects associated with BDE2 and TRPS and demonstrates the pathogenicity of TRPS1 missense variants located outside both the nuclear localization signal and the GATA ((A/T)GATA(A/G)-binding zinc-containing domain) and Ikaros-like binding domains. Unfortunately, we could not find distinctive phenotypic features that might have led to an earlier clinical diagnosis, further highlighting the high degree of overlap among skeletal syndromes associated with brachydactyly and AHO-like features, and the need for a close interdisciplinary workout in these rare patients.

Introduction

The vertebrate skeleton derives from the condensation of precursor mesenchymal stem cells that, after the patterning phase, differentiate to osteoblasts and chondrocytes. Genetic diseases of skeletal development and growth can derive from molecular alterations in different genes of the same pathway or in genes part of pathways involved in regulatory networks controlling skeletogenesis, like bone and cartilage cells patterning and differentiation, and in the growth‐plate development.( )

Brachydactyly (BD), a family of developmental disorders characterized by shortening of hands, feet or both, was classified on the anatomic and genetic basis into five groups and was added to the international Nosology and Classification of Genetic Skeletal Disorders in 2001.( )

Brachydactyly type E (BDE) is a rare BD subtype, presenting with the shortening of metacarpals/metatarsals, that can be found as an isolated trait or as part of several genetic syndromes.( ) Isolated BDE (OMIM #113300) derives from autosomal dominant HOXD13 alterations, although most BDEs occur in association with syndromic diseases, such as the Albright hereditary osteodystrophy (AHO), Turner syndrome, and others. Dysfunction of different components of the same network/module may result in overlapping phenotypes.( , )

Syndromes including BDE as a pathognomonic trait can roughly be divided according to the presence (pseudohypoparathyroidism [PHP], acrodysostosis type 1 [ACRDYS1]) or absence of multihormonal resistance (hypertension with BD syndrome [HTNB], BDE with short stature parathyroid hormone‐like hormone [PTHLH] type, brachydactyly‐mental retardation syndrome [BDMR], acrodysostosis type 2 [ACRDYS2], tricho‐rhino‐phalangeal syndrome [TRPS], and Turner syndrome [TS]).( )

The existence of pathways controlling specific cellular and tissue functions, leading to related and clinically overlapping genetic diseases, is exemplified by disorders due to the impairment of the parathyroid hormone (PTH)/parathyroid hormone‐related protein (PTHrP) signaling pathway. Defects affecting different members of this pathway determine several previously mentioned syndromic disorders, such as PHP, ACRDYS1 and ACRDYS2, and HTNB. According to a recent proposal, genetic diseases caused by defects of PTH/PTHrp signaling elements were named inactivating PTH/PTHrP signaling disorders (iPPSD), followed by numbers for specific subtypes according to the underlying molecular defect (iPPSD1, loss‐of‐function mutations in PTH1R; iPPSD2, loss‐of‐function mutations in GNAS; iPPSD3, methylation defects at one or more GNAS DMRs; iPPSD4, PRKAR1A mutations; iPPSD5, PDE4D mutations; iPPSD6, PDE3A mutations; iPPSDx, no molecular defects identified).( , ) Thus, iPPSDs represent an updated description and classification of heterogeneous, but still similar disorders, previously identified as pseudohypoparathyroidism and related disorders (including pseudopseudohypoparathyroidism, progressive osseous heteroplasia, and acrodysostosis) (Supplementary Table S1). In the absence of a molecular diagnosis, performing a correct identification of the different subtypes may be rather difficult in some cases because of interfamilial and intrafamilial variability.( , )

Two additional genes should be considered for establishing a correct differential diagnosis in patients presenting BDE within an “AHO‐like phenotype,” because both determine syndromic BDE associated with additional skeletal and endocrinological phenotypes: the PTHLH (OMIM*168470) and the TRPS1 (OMIM*604386) genes associated with “BDE with short stature” syndrome (BDE2; OMIM#613382) and tricho‐rhino‐phalangeal syndrome (TRPS1; OMIM #190350, and TRPS3; OMIM #190351), respectively.( , , , )

PTHLH encodes for PTHrP, a hormone mainly secreted by chondrocytes, perichondrial cells, and osteoblasts with many biological actions throughout life on cartilage, mammary gland, and dental development, on central nervous system activity, and skin and hair follicles. It binds to PTHR1 activating the Gsα‐cAMP‐PKA‐PDE4D pathway.( , , , ) Haploinsufficiency of PTHLH consequent to loss‐of‐function heterozygous mutations leads to premature differentiation of chondrocytes and abnormal bone size, with cone‐shaped epiphyses and their premature fusion.( ) The deriving syndrome, reported in 2010 by Klopocki and colleagues,( ) consists of patients with BDE and, in most cases, short stature. Additionally, the clinical presentation may include delayed dental eruption and/or oligodontia, hypoplastic nails, learning disabilities, altered mammary gland development, mild craniofacial dysmorphism, and delayed bone age.( , , ) To date, 12 mutations have been described, located along the whole coding region encoding for bioactive peptides obtained from the PTHLH prohormone.

Autosomal dominant mutations in the TRPS1 gene were described for the first time in 2000 by Momeni and colleagues,( ) who demonstrated their association with tricho‐rhino‐phalangeal syndrome, characterized by a triad of hair (sparse and slow‐growing hair), craniofacial (bulbous nose, long flat philtrum, thin upper vermilion border, protruding ears), and skeletal abnormalities (brachydactyly with phalangeal cone‐shaped epiphyses, hip malformations, short stature due to progressive postnatal growth retardation).( , , ) Since then, great variability in clinical findings has been reported even within the same family, and several additional skeletal alterations and systemic manifestations were observed in some patients.( , ) In particular, TRPS type I and type III subtypes, showing variable severity of skeletal abnormalities, were defined, which differed regarding the type of TRPS1 causative mutation. Type I derived from gene deletions and frameshift mutations, whereas the most severe type III from missense mutations in the GATA ((A/T)GATA(A/G)‐binding zinc‐containing domain) DNA‐binding zinc‐finger domain. Recently, the discovery of genetic alterations outside the GATA and the Ikaros‐type zinc‐finger encoding regions demonstrated the presence of additional functional key elements of TRPS1, like the nuclear localization signals (NLS), and the need to further investigate other protein domains.( , , ) To date, only four disease‐causing mutations in exon 3 were reported in the literature (one frameshift and three missense) (Supplementary Table S2). Nonsense and frameshift variants led to nonsense‐mediated mRNA decay and a dominant negative effect, respectively, whereas missense alterations were supposed to alter the interactions with RUNX2 and induce a milder phenotype.( , , )

The present work aimed to establish, by the use of whole‐exome sequencing (WES) and in vitro functional studies, a correct clinical and genetic diagnosis in syndromic patients with AHO‐like skeletal malformations. After a review of available clinical data we tried to identify, if present, those phenotypic manifestations useful for an earlier diagnostic definition of mutated patients. Moreover, to our knowledge, this is the first study demonstrating the pathogenicity of TRPS1 missense variants located outside both the NLS and the GATA and Ikaros‐like binding domains.

Patients and Methods

Patients

The present study includes five patients, referred to our laboratory after the observation of BDE and a complex phenotype including AHO features associated with resistance to PTH (n = 1) and without hormone resistances (n = 4). Patients were investigated for the presence of major (resistance to PTH, ectopic ossifications, and BDE) and minor (resistance to TSH, additional hormone resistances, motor and/or cognitive impairment, intrauterine and/or postnatal growth retardation, obesity/overweight and flat nasal bridge and/or maxillary hypoplasia and/or round face) criteria for iPPSD, but genetic testing failed to identify any molecular defect at iPPSD‐associated genes.

Clinical details of investigated patients are summarized in Supplementary Table S3. The study was performed in compliance with relevant legislation and institutional guidelines and was approved by the IRCCS Fondazione Cà Granda Ospedale Maggiore Policlinico Institutional Committee. All subjects involved in the study subscribed the informed consent for genetic and epigenetic studies.

Exome sequencing and bioinformatic analysis

Previous investigations included the direct sequencing of GNAS, PRKAR1A, and PDE4D coding regions, and PDE3A exon 4 (ENSEMBL reference sequence ENST00000354359, ENST00000589228, ENST00000502484 and ENST00000359062, respectively), and MS‐MLPA of the GNAS locus (ME031 GNAS probemix; MRC‐Holland, Amsterdam, The Netherlands) and they did not reveal (epi)genetic abnormalities. Moreover, we excluded 2q37 deletions associated with the AHO‐like syndrome, known as brachydactyly‐mental retardation syndrome (BMDR), by MLPA (P264‐Human Telomere‐9 probemix; MRC‐Holland, Amsterdam, The Netherlands).( , )

Libraries for WES were prepared using 50 ng of patients' DNAs and the SureSelectQXT Clinical Research Exome V2 kit (Agilent Technologies, Santa Clara, CA, USA), then 1.3 pmol of libraries were sequenced by paired‐end high‐throughput NGS on the NextSeq550 platform (Illumina, San Diego, CA, USA). Raw sequencing run outputs .bcl files were converted into FASTQ files and demultiplexed by the Illumina bcl2fastq software. The quality control showed the production of 2×150 nucleotides paired‐end read sequences with a mean region coverage depth of 165× and a coverage uniformity of 91%. The BaseSpace BWA Enrichment v2.1 App (Illumina, San Diego, CA, USA) was used to process data (alignment through BWA‐MEM v0.7.7 on GRCh37/hg19 and Picard v1.79) and obtain VCF files (GATK HaplotypeCaller v1.6). The eVai, Expert Variant Interpreter v0.6 tool (enGenome Srl), allowed the prioritization of potentially deleterious of single nucleotide variants (SNVs) and small insertions and deletions (indels) on the base of the functional relevance (i.e., missense, nonsense, indels, and splice‐site changes), taking into account both dominant and recessive inheritance models. Selected damaging rare variants were confirmed by direct sequencing and we verified their presence in publicly available databases, The Human Gene Mutation Database (HGMD®), Leiden Open Variation Database (LOVD), and ClinVar, and in the literature (Fig. 1). Mutation nomenclature followed the Human Genome Variation Society (HGVS) guidelines, based on the nucleotide and protein numbering of the Locus Reference Genomic (LRG) sequence format adopted by LOVD.( )

Fig. 1

(A,B) Electropherograms of the confirmatory direct sequencing of genetic variants discovered by WES. (C) Schematic representation of the PTHLH isoform NM_198965.1, its proteolytic processing pattern into bioactive peptides and known mutations mapping in coding regions. Prohormone: 177 amino acids. Bioactive peptides: mature N‐terminal (having PTH‐like and growth regulatory activities) (amino acids 37–70), mid‐region (regulating calcium transport and cell proliferation and containing the nuclear localization signal, NLS) (amino acids 103–130) and C‐terminal secretory (modulating osteoclast activity) (amino acids 143–146) peptides. (D) Schematic representation of the TRPS1 isoform NM_14112.2, protein domains and known mutations in coding regions. Protein: 1294 amino acids. Protein domains of the TRPS1 transcription factor predicted by SMART: at the N‐terminal, region involved in the nuclear translocation, 7 zinc finger motifs (amino acids 235–260/346–371/447–472/536–567/627–650/679–702/705–728); one DNA‐binding zinc finger GATA‐type motif (amino acids 903–953), surrounded by the putative NLSs; at the C‐terminal, region involved in the homodimer formation, and 2 zinc finger IKAROS‐like (amino acids 1228–1250/1256–1280). Known coding mutations are reported in the upper part (mutations found in our patients reported in bold; @ highlights recurrent mutations). NLS = nuclear localization signal; SMART = simple modular architecture research tool; WES = whole‐exome sequencing.

Cell lines for in vitro functional studies of novel TRPS1 variants

The functional effect of novel TRPS1 variants, located outside the already confirmed disease‐associated protein domains, was determined using the human hepatocellular carcinoma cell line HepG2 (LGC Standards, Sesto San Giovanni, Italy), selected after database search on Protein Atlas, which confirmed to be a natural TRPS1 knockout model, and in the human mesenchymal stem cell (MSC) L88/5 cell line (kindly provided by Dr. Karin Thalmeier), that are TRPS1‐expressing cells used to evaluate the effect of TRPS1 overexpression and of the mutant heterozygous condition in precursors cells of both osteoblasts and chondrocytes. HepG2 offered the chance to work in the absence of the confounding role of a basal TRPS1 expression, which could affect the effect of mutated TRPS1 on SOX9 and STAT3 transcription repression (Fig. S1). The L88/5, on the other hand, allowed us to investigate the potential dominant‐negative effect of our novel TRPS1 variants, located outside known disease‐associated protein domains, and their effect on RUNX2 expression, a fundamental differentiation factor involved in both bone and cartilage development.( )

Cells culture conditions and transfection

HepG2 were cultured on collagen‐coated plates in IMDM+Glutamax® supplemented with 20% FBS, 1% non‐essential amino acids, and 1% penicillin/streptomycin (Sigma‐Aldrich, Milan, Italy). L88/5 were cultivated in RPMI‐1640 supplemented with 10% FBS and 1% of penicillin/streptomycin (Sigma‐Aldrich, Milan, Italy).

Mutant plasmids were obtained from 10 ng of wild‐type pCMV6‐AC‐GFP encoding the full length TRPS1 transcript (RG215856: TRPS1 NM_014112 Human Tagged ORF Clone; OriGene, Rockville, MD, USA) by site‐directed mutagenesis using the QuickChange II XL Site‐Directed Mutagenesis kit (Agilent Technologies, Cedar Creek, TX, USA) and 125 ng of mutation‐specific primers to introduce our novel c.343C>T and the c.422C>T variants, the c.2894G>A variant in exon 7 (mutational hot spot located in the NLS2 domain, used as mutated positive control) and the already described variants c.47G>A and c.203A>C (known variants were renamed according to the reference sequence NM_014112; primer sequences available upon request).( , , ) Transformed ampicillin‐resistant E. coli colonies were recovered from agar plates and expanded over night and plasmids were purified by the Plasmid Plus Maxi Kit (Qiagen, Valencia, CA, USA) and tested by Sanger sequencing. Twenty thousand (20,000) cells/cm2 of HepG2 and 15,000 cells/cm2 of L88/5 were grown (each condition was performed in triplicate) until 80% confluence, transfected with 1 μg of wild‐type and mutated pCMV6‐AC‐GFP TRPS1 expression vectors by the Mega Tran 2.0 transfection reagent (OriGene, Rockville, MD, USA). As control, only Mega Tran 2.0 was applied. The observed transfection efficiency was 74% in agreement with the manufacturer's expected ranging from 51% up to 79%.

Gene expression in transfected cells

RNAs were collected, using the Trizol protocol followed by treatment with Dnase, after 48 hours in both HepG2 and L88/5 cell lines. cDNAs obtained by using the PrimeScript™ RT Master Mix Perfect Real Time (Takara Bio USA, Inc., San Jose, CA, USA) were analyzed for TRPS1, SOX9, STAT3, RUNX2, RNASEP, and GAPDH gene expression by qPCR (QuantStudio3; Applied Biosystems, Foster City, CA, USA) using the ΔΔCt method (primers sequence and genomic location available upon request). Data analysis by the Expression Suite Software v.1.3 (Applied Biosystems) allowed to determine the relative fold expression change (relative quantification [RQ] reported as log2 fold changes) of TRPS1, SOX9, STAT3, and RUNX2 in transfected cells, using mock cells as reference (Figs. 2 B, 3).

Fig. 2

Gene expression analysis in transfected the HepG2 cell line and nuclear and cytoplasmatic localization of TRPS1 GFP‐tagged wt and mutated proteins. (A) By immunocytochemistry the cellular localization of TRPS1 GFP‐tagged (green) transcription factors after 48 hours transfection of pCMV6‐AC‐GFP TRPS1 plasmids in HepG2. Nuclei were counterstained by DAPI (blue). The graph summarizes the quantification of TRPS GFP‐tagged by ImageJ software. Scale bars = 50 μm. Values obtained by the TRPS1 localization analysis of mutants were statistically significant (p < 0.001, Student’s t test). (B) Gene expression analysis of TRPS1 target genes, SOX9 and STAT3 in HepG2 pCMV6‐AC‐GFP TRPS1 transfected. The y‐axis reports the binary logarithm (ln2) RQ in mRNA expression versus mock, used as calibrator reference. Statistical analysis: ordinary one‐way ANOVA. SOX9: F = 95.75, p <0.0001; multiple comparisons p values: wt = <0.0001, c.47G>A = 0.1434, c.203A>C = 0.8592, c.343C>T* = 0.0032, c.422C>T* = 0.0002, c.2894G>A = <0.001. STAT3: F = 804.5, p = <0.0001; multiple comparisons p values: wt = <0.0001, c.47G>A = >0.9999, c.203A>C = 0.0008, c.343C>T* = <0.0001, c.422C>T* = 0.0009, c.2894G>A = 0.9956. RQ = relative fold expression change.

Fig. 3

RUNX2 gene and protein expression analysis in the L88/5 cell line by immunocytochemistry. (A) RUNX2 expression by rhodamine immunostaining (red) in L88/5 (36 hours after pCMV6‐AC‐GFP TRPS1 plasmids transfection) and nuclear counterstaining by DAPI (blue). The graph resumes the staining quantification by ImageJ tool. Scale bars = 50 μm. RUNX2 immunocytochemistry values in mutants were statistically significant, but the one for the c.343C>T variant (p < 0.001, Student's t test). (B) Gene expression analysis of RUNX2 in wild‐type and pCMV6‐AC‐GFP TRPS1 transfected L88/5 cells. The y‐axis reports the binary logarithm (ln2) RQ in mRNA expression versus mock, used as calibrator reference. Statistical analysis: ordinary one‐way ANOVA. RUNX2: F = 220.9, p = <0.0001; multiple comparisons p values: wt = <0.0001, c.47G>A = 0.0368, c.203A>C = 0.0846, c.343C>T* = <0.0001, c.422C>T* = 0.0011, c.2894G>A = <0.0001. RQ = relative fold expression change.

Cellular localization and dominant‐negative effect of new TRPS1 variants

To determine the molecular mechanism underlying the defective functioning due to novel sequence alterations of TRPS1, we investigated the ability of TRPS1 green fluorescent protein (GFP)‐tagged transcription factors to be translocated into the nucleus in HepG2 cells and the dominant‐negative effect on the wild‐type allele in L88/5 cells.

For HepG2 nuclear counterstaining, cells were fixed 48 hours posttransfection in cold acetone (5 minutes at room temperature [RT]), permeabilized with 0.3% of Triton X‐100 (Sigma‐Aldrich) in PBS (15 minutes at RT), then incubated with 4′,6‐diamidino‐2‐phenylindole (DAPI) (Sigma‐Aldrich, Milan, Italy) (1 hour at) RT (Fig. 2). For L88/5 staining, cells were fixed 48 hours posttransfection in cold acetone (5 minutes at RT), permeabilized with 0.3% of Triton X‐100/PBS (15 minutes at RT), then incubated with 1% bovine serum blocking solution (30 minutes at RT). Immunocytochemistry (IF) was performed using the rat monoclonal immunoglobulin G (IgG) RUNX2 primary antibody (MAB2006; R&D systems, Minneapolis, MN, USA) and the Alexa Fluor 546 Goat anti‐rat IgG secondary antibody, while nuclei were stained with DAPI (Fig. 3). Images were acquired by Zeiss AxioObserver microscope equipped with high‐resolution digital video camera (AxioCam; Carl Zeiss Microscopy, Inc., Dublin, CA, USA) and Apotome system for structured illumination and recorded by AxioVision software 4.8. Quantifications of TRPS1 cellular localization (n = 6/group), calculated whole cell surface area versus nucleus area, and of RUNX2 expression (n = 9/group) were determined by the Image J software tool (NIH, Bethesda, MD, USA; https://imagej.nih.gov/ij/) (Figs. 2 and 3).

Results

Identification of new genetic variants and bioinformatic analysis

WES identified three different PTHLH and two TRPS1 mutations in five patients (Fig. 1). More than 70,000 variants were identified in each patient, which were filtered and prioritized according to pathogenicity scores computed using a subset of 18 out of the 28 American College of Medical Genetics and Genomics (ACMG)/Association for Molecular Pathology (AMP) criteria for variant pathogenicity assessment by the eVai tool.( ) The eVai benchmarking analysis of ACMG/AMP preclassification sensitivity (pathogenic/likely pathogenic prediction) and specificity (benign/likely benign prediction) were determined as 78% and 94%, respectively.

Genetic variants in the PTHLH gene (patient 1: c.2T>C, p.Met1?; patient 2: c.166C>T, p.Arg56X; patient 3: c.299del, p.Glu100Glyfs*5) showed high pathogenicity scores (PS = 6 and 7). The ACMG/AMP criteria evaluated these variants as “null variants” (missense, nonsense, frameshift, canonical ±1 or 2 splice sites, initiation codon, single or multiexon deletion) in a gene with probability of being LoF intolerant (pLI) ≥0.9 (score predicting the pLI for dominant inheritance conditions, developed by the Exome Aggregation Consortium [ExAC]). They were absent in controls from population databases (Genome Aggregation Database [gnomAD], ExAC, and 1000 Genomes Project [1000GP]) and presented different functional, conservational, and evolutionary scores (PaPI [pseudo amino acid composition to score human protein‐coding variants] ≥0.5, DANN [deleterious annotation of genetic variants using neural networks] ≥0.9 and dbscSNV [database developed for functional prediction and annotation of all potential non‐synonymous single‐nucleotide variants] ≥0.9) predicting damaging variants. Both the p.R56X variant, located in the region encoding the PTH‐like domain, and the p.E100Gfs*5 variant, inside the mid‐region domain containing the NLS, are predicted to affect the correct translation and lead to a truncated and nonfunctional protein product. Even though we did not check PTHLH mRNA levels for the unavailability of a patient sample, it is known that nonsense variants result in truncated, incomplete, and usually nonfunctional protein products. Many organisms, including humans, degrade mRNAs containing nonsense variants before translation by a mechanism called nonsense‐mediated mRNA decay (NMD), in order to avoid the creation nonfunctional, possibly even harmful, polypeptide products. The p.M1? variant, affecting the first methionine, probably led to the complete loss of translation from one allele (Fig. 1). The potential outcome of this variant was further investigated by online bioinformatic tools producing predictions of translation start in vertebrate (ATGpr [https://atgpr.dbcls.jp], DNA Functional Site Miner [http://dnafsminer.bic.nus.edu.sg], NetStart [https://services.healthtech.dtu.dk], and PreTIS [http://service.bioinformatik.uni-saarland.de/pretis]). All identified in‐frame translation initiation sites (TIS) were found in the 5′ untranslated region (UTR) or the 3′UTR regions and, according to the location of stop codons, the expected open reading frames (ORFs) could encode for small polypeptides and they did not include codons for PTHLH functional domains.

Custom filtering of WES data from patients 4 and 5 succeeded in finding two rare variants in TRPS1 (patient 4: c.343C>T, p.His115Tyr; patient 5: c.422C>T, p.Pro141Leu), that presented a pathogenicity score of 2. According to the ACMG/AMP criteria, TRPS1 molecular alterations were considered pathogenetic being missense variants in more than half of the transcripts of a gene with a low rate of benign missense variation according to ExAC (Z score >0) and in which missense variants are a common mechanism of disease, with more than 10 pathogenic variants reported in ClinVar, of which at least half were missense. Moreover, a deleterious effect on the gene product was predicted by PaPI or DANN or dbscSNV scores.

Most of the missense TRPS1 mutations described up to now are located in the GATA and Ikaros‐like DNA‐binding zinc‐finger domains, and only few pathogenetic variants at other exons encoding for additional known protein domains were reported (Fig. 1). Only four TRPS1 variants at exon 3, different from those identified in our patients, were reported (three missense and one frameshift), different from those identified in our patients, further supporting the predicted causative effect of our p.H115Y and p.P141L variants, although located outside protein domains already associated with the pathological phenotype development (Supplementary Table S2 and Supplementary References).

In vitro functional studies

Functional investigations started confirming the absence of basal TRPS1 expression and the concomitant presence of SOX9 and STAT3 target genes in the HepG2 cell line (Fig. S1). Mutant TRPS1 transcription factors lost their physiological repressor activity on target genes compared to the wild‐type (SOX9: c.47G>A RQ = 1.199 ± 0.273, c.203A>C RQ = 0.882 ± 0.336, c.343C>T RQ = 1.417 ± 0.433, c.422C>T RQ = 1.556 ± 0.129, c.2894G>A RQ = 3.167 ± 0.831 versus mock RQ = 0.644 ± 0.171 and WT RQ = −1.773 ± 0.493; STAT3: c.47G>A RQ = 1.456 ± 0.305, c.203A>C RQ = 0.733 ± 0.044, c.343C>T RQ = 3.232 ± 0.055, c.422C>T RQ = 0.821 ± 0.072, c.2894G>A RQ = 1.522 ± 0.164 versus mock RQ = 1.453 ± 0.332 and WT RQ = ‐6.264 ± 0.396), supporting a pathogenic role (Fig. 2).

We then investigated the molecular mechanism underlying the loss‐of‐function transcriptional effect of both novel c.343C>T and c.422C>T and previously reported c.47G>A, c.203A>C, and c.2894G>A TRPS1 variants and we found that, in the HepG2 cell line, they were responsible for a defective nuclear translocation of the mutated proteins. In particular, thanks to the GFP tag of plasmid‐derived proteins, we observed a prevalent nuclear localization of the wild‐type (nuclear versus cytoplasmatic = 70% versus 30%) against a prevalent cytoplasmatic localization of mutated TRPS1 (nuclear versus cytoplasmatic = c.47G>A 22% versus 78%, c.203A>C 17% versus 83%, c.343C>T 28% versus 72%, c.422C>T 12% versus 88% and c.2894G>A 20% versus 80%) (Fig. 2). Finally, we tested whether these variants determined a dominant negative effect on the protein product deriving from the wild‐type allele; ie, the physiological condition observed in heterozygous TRPS patients. To this scope, we used the L88/5 cell line. We transfected L88/5 cells and the subsequent gene expression and immunocytochemistry analysis confirmed the loss of TRPS1 repressory effect on RUNX2. In particular, mutant plasmids were not able to downregulate RUNX2 mRNA expression (mock RQ = 0.755 ± 0.239 versus WT RQ = −1.547 ± 0.844 c.47G>A RQ = 1.567 ± 0.307 versus c.203A>C RQ = 1.459 ± 0.294 versus c.343C>T RQ = 2.761 ± 0.223 versus c.422C>T RQ = 1.755 ± 0.364 versus c.2894G>A RQ = 8.803 ± 0.276), even though the most effective loss of regulatory effect on RUNX2 was exerted at protein level, as demonstrated by the statistically significant increase of IF positivity (mock about 20%, wild‐type <10%, mutants about 30% to 50%) (Fig. 3).

Patients' phenotype

Patients reported in the present study were, at first examination, suspected to be affected by iPPSD because of BDE associated with other AHO features and, in one case, PTH resistance.

The PTHLH mutated patients' phenotype included: brachydactyly type E in all patients (shortening of 3rd, 4th, and 5th metacarpals and, additionally, of 3rd and 4th metatarsals and 2nd, 3rd, 4th, and 5th phalanges), associated with short stature and delayed dental eruption and/or numerical teeth abnormalities (patients 2 and 3). Additionally, we observed rhizomelic shortening of limbs and skull vault thickening in patient 1, cafè‐au‐lait spots in patient 2, and insulin resistance, obesity, posterior plagiocephaly, Chiari malformation II, pityriasis rosea, and striae rubrae in patient 3 (Table 1). In the paternal branch of patient 1’s family other cases of short stature, including her brother who also had brachydactyly, were reported. The family history of patient 2 included: (A) maternal grandfather, diabetes mellitus type 2 (DMT2) and colon cancer; (B) paternal grandfather, liver cancer; and (C) father, paternal aunt, paternal grandmother, short stature and apparent brachydactyly. Finally, anamnesic notes of patient 3 were DMT2 (paternal grandfather); plagiocephaly at birth and Chiari II malformations (mother).

Table 1

Clinical Presentation of Patients at the Moment of the Clinical Diagnosis

Parameter	Patient 1 PTHLH	Patient 2 PTHLH	Patient 3 PTHLH	Patient 4 TRPS1	Patient 5 TRPS1
Sex	F	F	M	F	F
Age at clinical diagnosis (years)	34	10	17	13	13
Major iPPSD criteria
rPTH (PTH; Ca; P)	No (55.8 pg/mL; 9.2 mg/dL; 3.6 mg/dL)	No (25.4 pg/mL; 10.5 mg/dL; 3.8 mg/dL)	No (45 pg/mL; 10 mg/dL; 3.85 mg/dL)	Yes (212 pg/mL; 6.3 mg/dL; 4 mg/dL)	NA
EO	No	No	No	No	No
Br	Yes (brachymetacarpia, most pronounced in IV and V; brachymetatarsia, particularly pronounced in IV, shortening of the middle phalanges)	Yes (brachymetacarpia; dysmorphic 2nd phalanges II and V; brachymetatarsia III; brachyphalangy of toes)	Yes (brachymetacarpia left IV and V and right III, IV and V; bilateral brachyphalangy 1st II, III and IV)	No	Yes (metacarpals and 3rd phalanges)
Minor iPPSD criteria
rTSH	No	No	No	NA	Subclinical hypothyroidism (TSH levels NA)
Additional HR	No	No	No	No	No
M/C impairment	No	No	No	No	Yes (IQ45‐IQP46‐IQV55; ADD; ODD)
Short stature	Yes	Yes	Yes	No	Yes (H NA)
OB/OW	No	No	Yes (W=88.6kg BMI=32.7kg/m2)	No	Yes (W NA)
DF	No	No	No	No	No
Additional clinical features	Rhizomelic shortening of the upper limbs (arm span = 135 cm); lower limb shortening; skull vault thickening	Delayed eruption of definitive molar teeth; cafè‐au‐lait spots	Posterior plagiocephaly (surgery); Chiari malformation II (surgery); pityriasis rosea; striae rubrae; supernumerary tooth; insulin resistance	No	Advanced bone age; motor stereotypies; hypertrophic pyloric stenosis; incomplete elbow joint extension; transitional vertebra at the lumbosacral junction; levoscoliosis

ADD = attention deficit disorder; Br = brachydactyly; Ca = serum calcium (reference values = 9‐10.5 mg/dL); DF = dysmorphic facies; EO = ectopic ossification; F = female; H = height; HR = hormone resistances; IQ = intelligence quotient; IQp = performance intelligence quotient; IQv = verbal intelligence quotient; M = male; M/C = motor and/or cognitive; NA = data not available; OB/OW = obesity/overweight; ODD = oppositional defiant disorder; P = serum phosphate (2.8‐4.5 mg/dL); PTHLH = parathyroid hormone like hormone; rPTH = PTH resistance (reference values=10‐65 pg/mL); rTSH = TSH resistance; TRPS1 = tricho‐rhino‐phalangeal syndrome 1; W = weight.

The phenotype of patients bearing TRPS1 mutations included PTH resistance in patient 4, although, due to the lack of an informative post‐therapeutical follow‐up, we were not able to exclude that such resistance derived from vitamin D insufficiency, and brachydactyly, short stature, intellectual disability, obesity, subclinical hypothyroidism, advanced bone age with cone‐shaped epiphysis, abnormal joint mobility, scoliosis, a transitional vertebra at the lumbosacral junction, motor stereotypies and the hypertrophic pyloric stenosis in patient 5 (Table 1). We retrieved few data on the family history of patient 4: positivity for ischemic heart disease at a young age was reported and her sister was affected by partial epilepsy, coronary heart disease (CHD) with aortic subvalvular stenosis and learning difficulties.

Parental biological samples were not available for genetic testing, so we could not determine whether the detected variants were inherited or de novo and, consequently, their inheritance pattern.

Discussion

At first examination patients reported in the present study, except patient 4, presented with BDE as the main clinical finding with the addition of some features typical of the AHO phenotype; thus, the first clinical hypothesis was that they were affected by AHO. In particular, patient 4 was proposed to suffer from pseudohypoparathyroidism because of the detection of mild PTH resistance. Genetic analysis did not reveal any (epi)genetic alteration at GNAS, PRKAR1A, PDE4D, and PDE3A genes nor 2q37 deletions.

After WES, patients were rediagnosed as affected by BDE with short stature or TRPS, keeping in mind the known significant variability and incomplete penetrance of BDE2 and TRPS‐related symptoms. Four out of five of the discovered PTHLH and TRPS1 variants were novel to the literature, under the notion that most mutations are private. Only the nonsense PTHLH p.Arg56X in patient 2, first described in 2016 by Jamsheer and colleagues,( ) is now considered the only known recurrent pathogenic variant.( )

The PTHLH gene is translated into a precursor prohormone containing the signal peptide and the bioactive peptides: the N‐terminal domain, with PTH‐like and growth regulatory activities, the mid‐region containing the NLS, that regulates calcium transport and cell proliferation and the C‐terminal domain, which modulates the osteoclast activity. The mutations found in patients 1 (p.Met1?), 2 (p.R56X), and 3 (p.E100Gfs*5) affected the translation of the whole protein, the PTH‐like domain and the NLS‐containing domain, respectively, with the last two leading to the production of truncated products.

All PTHLH mutated patients presented BDE and short stature, which are clinical features considered hallmarks of BDE2. Moreover, patients 2 and 3 presented delayed dental eruption and/or numerical tooth abnormalities, that are frequently but not obligatory observed in BDE2.( , , , ) The additional features reported in patient 3 (ie, insulin resistance, obesity, posterior plagiocephaly, Chiari malformation II, pityriasis rosea, and striae rubrae) were not previously clearly associated with BDE2.

TRPS1 is a tissue‐specific transcriptional repressor with nine zinc‐finger domains, which regulates differentiation, proliferation, and apoptosis in cartilage, kidneys, and hair follicles. To date, most reported missense mutations are located in GATA DNA‐binding and Ikaros‐like zinc‐finger domains, and only two mutations in exon 3, outside the regions encoding for known disease‐associated functional domains, were reported (Supplementary Table S2 and Supplementary References). TRPS1 mutations discovered in our patients (p.H115Y and p.P141L) were located in exon 3, containing a Kozak consensus ATG translation start site and the first C2H2‐type zinc‐finger domain, and computational pathogenicity analysis predicted a damaging effect. Besides, the p.P141L variant was located near the highly conserved phosphoserine residue at position 140, further supporting a possible damaging effect. The pathogenetic role of p.H115Y and p.P141L amino acid substitutions in the disease development was confirmed by in vitro functional studies, demonstrating the loss of transcriptional repressor activity on SOX9 and STAT3 target genes compared to the wild‐type TRPS1. In particular, our investigations demonstrated that novel TRPS1 variants here described affected the nuclear translocation of the TRPS1 transcription factor and that mutated TRPS1 exerts a dominant negative effect on residual wild‐type protein in heterozygous cells. Note that the regulatory effect of TRPS1 on RUNX2 is mainly indirectly exerted by cellular mechanisms able to reduce the RUNX2 protein lifespan rather than its gene repression.

At the first clinical evaluation, patients bearing TRPS1 mutations were suspected to have pseudohypoparathyroidism and AHO because of the presence of PTH resistance in patient 4 and of BD, short stature, intellectual disability, obesity and subclinical hypothyroidism in patient 5. Patient 4 also showed low vitamin D levels, but the association of hypocalcemia with markedly elevated PTH serum levels (>200 pg/mL) is highly evocative of hormone resistance, although, unfortunately, no data after vitamin D supplementation were available. Even though resistance to the action of PTH and obesity are not typical TRPS features, single affected TRPS cases have already been reported in the literature. Both Böhles and Ott( ) and Pereda and colleagues( ) described subjects with partial resistance to PTH without associated hypocalcemia and obesity, although this last finding may be a confusing factor due to pandemic obesity. Unfortunately, few information about patient 4 phenotype was available, whereas patient 5 also presented other TRPS‐associated signs, such as advanced bone age with cone‐shaped epiphysis, abnormal joint mobility and scoliosis, that helped in confirming the pathogenicity of the novel TRPS1 variant. In conclusion, this work allowed to establish a correct genetic diagnosis, followed by an adequate follow‐up and genetic counseling, in five patients with syndromic BDE and to discover four novel mutations in PTHLH and TRPS1 genes, further expanding the spectrum of genetic defects associated with two extremely rare conditions such as BDE2 and TRPS. To our knowledge, we first experimentally demonstrated by a functional study that TRPS1 missense variants located outside known fundamental functional protein domains can alter the normal functioning of the TRPS1 transcription factor and associate with mild‐to‐moderate TRPS phenotypes. No genotype‐phenotype correlations were observed, both considering the location and the type of the genetic variant. Unfortunately, we could not find distinctive phenotypic features that might have helped, before the use of broad genetic analysis by NGS, in early identification of BDE2 and TRPS, further highlighting the high degree of overlap among the different (but similar) syndromes associated with brachydactyly and other AHO‐like features, and again stressing the need for a close interdisciplinary workout in rare patients.

Author Contributions

Francesca Marta Elli: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing; Deborah Mattinzoli: Data curation, Formal analysis, Methodology, Resources, Validation, Writing – review & editing; Camilla Lucca: Data curation, Investigation, Methodology, Resources, Writing – review & editing; Matteo Piu: Formal analysis, Investigation; Maria A. Maffini: Formal analysis, Investigation, Writing – review & editing; Jole Costanza: Data curation, Formal analysis, Writing – review & editing; Laura Fontana: Data curation, Writing – review & editing; Carlo Santaniello: Formal analysis, Investigation, Resources; Concetta Forino: Formal analysis, Investigation, Writing – review & editing; Donatella Milani: Formal analysis, Investigation, Resources, Writing – review & editing; Maria Teresa Bonati: Formal analysis, Investigation, Resources, Writing – review & editing; Andrea Secco: Formal analysis, Investigation, Writing – review & editing; Roberto Gastaldi: Formal analysis, Investigation, Writing – review & editing; Piergiorgio Messa: Resources, Supervision, Writing – review & editing; Monica Miozzo: Resources, Supervision, Writing – review & editing; Maura Arosio: Supervision, Writing – review & editing.

Conflict of Interests

All the authors declare the absence of any potential conflicts of interest.

Web Resources

OMIM, http://www.omim.org/; ENSEMBL, https://www.ensembl.org; HGMD®, http://www.hgmd.cf.ac.uk/; LOVD, at https://databases.lovd.nl/; ClinVar, https://www.ncbi.nlm.nih.gov/clinvar/; eVai, https://www.engenome.com/; Illumina, https://support.illumina.com/downloads/; Protein Atlas, https://www.proteinatlas.org/ENSG00000104447‐TRPS1/cell; Expression Suite, https://www.thermofisher.com/it/en/home/technical‐resources/software‐downloads/expressionsuite‐software.html; Image J, https://imagej.net; ACMG/AMP guidelines, https://www.amp.org/clinical‐practice/practice‐guidelines/; SMART, http://smart.embl.de; PreTis, https://service.bioinformatik.uni‐saarland/pretis/background.php; NetStart, www.cbs.dtu.dk/services/NetStart/; DNAFSminer, https://dnafsminer.bic.nus.edu.sg; ATGpr, https://atgpr.dbcls.jp.

Peer Review

The peer review history for this article is available at https://publons.com/publon/10.1002/jbmr.4490.

Appendix S1 Supporting Information

Click here for additional data file.