Protein informatics combined with multiple data sources enriches the clinical characterization of novel TRPV4 variant causing an intermediate skeletal dysplasia

Abstract Background Transient receptor potential cation channel subfamily V member 4 (TRPV4) is an ion channel permeable to Ca2+ that is sensitive to physical, hormonal, and chemical stimuli. This protein is expressed in many cell types, including osteoclasts, chondrocytes, and sensory neurons. As such, pathogenic variants of this gene are associated with skeletal dysplasias and neuromuscular disorders. Pathogenesis of these phenotypes is not yet completely understood, but it is known that genotype–phenotype correlations for TRPV4 pathogenic variants often are not present. Methods Newly characterized, suspected pathogenic variant in TRPV4 was analyzed using protein informatics and personalized protein‐level molecular studies, genomic exome analysis, and clinical study. Results This statement is demonstrated in the family of our proband, a 47‐year‐old female having the novel c.2401A>G (p.K801E) variant of TRPV4. We discuss the common symptoms between the proband, her father, and her daughter, and compare her phenotype to known TRPV4‐associated skeletal dysplasias. Conclusions Protein informatics and molecular modeling are used to confirm the pathogenicity of the unique TRPV4 variant found in this family. Multiple data were combined in a comprehensive manner to give complete overall perspective on the patient disease and prognosis.

TRPV4-associated skeletal dysplasias. It has been noted that TRPV4 has high expression in sensory neurons, osteoclasts, osteoblasts, and chondrocytes . Of particular interest are the chondrocytes, cells that produce cartilage and compose early bones in the process of endochondral ossification. In wild-type mammalian fetuses, chondrocytes produce a cartilage matrix that is slowly replaced and lengthened by bone, a process that is continued after birth (Mackie, Ahmed, Tatarczuch, Chen, & Mirams, 2008). TRPV4 pathogenic variants have been experimentally proven to cause incomplete endochondral ossification by interfering with proper chondrocyte function Saitta et al., 2014;Weinstein, Tompson, Chen, Lee, & Cohn, 2014). It is proposed that Ca 2+ entry through mutant TRPV4 triggers overproduction of follistatin, an inhibitor of bone morphogenetic protein 2 (BMP2). BMP2 is necessary to promote production of several proteins that advance chondrocytes into the hypertrophic stage, a transition that precedes mature bone in endochondral ossification. With an increased follistatin concentration, chondrocytes are less sensitive to BMP2 signaling, meaning that bone is produced more slowly and less consistently than in the wild type (Leddy, McNulty, Guilak, & Liedtke, 2014). This proposed mechanism helps to explain how skeletal dysplasia can result from pathogenic variants of a gene that is expressed in multiple specialized cells.
The six skeletal dysplasias linked to TRPV4 pathogenic variants vary widely in severity. Familial digital arthropathy-brachydactyly is the mildest form, presenting only with progressive arthropathy and swelling of the hands and feet starting in early childhood. On the more severe end of the spectrum is metatropic dysplasia, which involves significant pelvic, spinal, and long bone abnormality and results in shortlimb short-stature dwarfism or premature death (Schindler et al., 1993). Our proband falls somewhere on the spectrum between these extremes presenting with a phenotype that shares aspects of spondyloepiphyseal dysplasia, Maroteaux type, and spondylometaphyseal dysplasia Kozlowski type. Both of these phenotypes are considered TRPV4-associated disorders of intermediate severity. The proband and her daughter present with the novel pathogenic TRPV4 variant NC_000012.12:c.2401A>G or NC_000012.12:p. Lys801Glu (referred to hereafter as p.K801E) and share some phenotypical features including short digits, abnormal nail beds, and short stature. We present the case of the proband and her family and confirm the pathogenicity of the newly discovered TRPV4 variant with molecular modeling techniques.

| CASE STUDY
Our proband, a 47-year-old female, initially presented to our clinic due to a lengthy history of skeletal abnormalities, bone pain, and osteoporosis. She has relatively short stature with a height of 5′2″. She reported that she had normal weight and height at birth with no notable bone abnormalities. As a child, she had scoliosis requiring bracing for correction. Since then, she has experienced "constant" bone and joint pain and her medical records indicated a history of anemia, vitamin D deficiency, and severe osteonecrosis (particularly in her hips). Physical examination of her fingers yielded multiple findings, including widened proximal interphalangeal joints, swan neck deformities, short broad digits, hyperplastic nails, and short nail beds (Figure 1a,b). Her hand joints were hypermobile, her feet were small, and her thorax was abnormally shaped (Figures 1c and 2a). We were also able to confirm her complaints of limited neck and back movement, as she had reported these areas were significant sources of pain. Overall, physical examination suggested a skeletal dysplasia as the source of the proband's symptoms. Skeletal survey and whole exome sequencing + mtDNA analysis were performed to determine the specific disorder. Patients' pedigree was determined for additional information ( Figure 3). Skeletal X-rays found multiple abnormalities in the proband, most notably in her spine, pelvis, and long bones. The proband's spine was notable for flattened vertebral bodies, narrowed disk spaces, and end plate sclerosis (Figure 2a). Her pelvis X-ray was abnormal ( Figure 2b). Her long bones exhibited metaphyseal abnormalities. The survey also provided skeletal confirmation for the shortened digits that had been apparent in the proband's physical examination. Additionally, the proband seemed to have diffuse arthropathy of the spine and large/small joints.
The proband's mtDNA analysis was negative, but whole exome sequencing uncovered a novel, likely pathogenic variant that explained her condition. The TRPV4 variant c.2401A>G (p.K801E) was found in a heterozygous state, consistent with the known autosomal-dominant transmission of TRPV4-associated disorders (Schindler et al., 1993). It is suspected that the proband received this variant from her father, who died at age 33. Per report, he had abnormal nails, abnormal hips, and short digits. The proband, in turn, passed this variant to her 22-year-old daughter, who has known phenotypical features include surgically corrected S-shaped scoliosis, brachydactyly, macrocephaly, and short stature at 4′11″. The differences in phenotype seen between these family members are unsurprising, as TRPV4 pathogenic variants do not exhibit clear genotype-phenotype correlations .
Given the skeletal and genetic findings in the proband, she ultimately fell on the disease spectrum somewhere between Maroteaux type and Kozlowski type spondyloepiphyseal dysplasia. At this point in time, her hip pain had increased enough to warrant surgical intervention, which was eventually pursued in spite of her significant osteoporosis. A right total hip arthroplasty was performed, and as of current date has not been complicated by fractures or new osseous lesions. The proband reported a significant improvement as compared to her preoperative level of pain and continues to follow with the clinic to monitor any further degenerative skeletal changes.

| Molecular modeling
The Transient receptor potential cation channel subfamily V member 4 (Trpv4 or TRPV4), called TRPV4 gene, is a nonselective calcium cation channel that regulates osmotic sensitivity and mechanical sensitivity in humans. The channel is regulated via a calmodulin-dependent mechanism with a negative feedback loop. It also plays a function in nonselective cation channel activation induced by 4-alpha-phorbol 12, 13-didecanoate as well as hypotonic stimulation in synoviocytes and also regulates production of IL-8 (Garcia-Elias, Lorenzo, Vicente, & Valverde, 2008;Strotmann, Schultz, & Plant, 2003). Lastly, the channel is inhibited through binding with phosphatidylinositol-4, 5-bisphosphate (Itoh et al., 2009). TRPV4 was taken from the NCBI Reference Accession Sequence: NP_067638: version NP_067638.3, which is encoded for the amino acid sequence: where the variant of unknown significance (VUS) is p.K801E. The NP_067638.3 sequence was used for computer-assisted modeling to evaluate the differences between the various X-ray structures for TRPV4 protein and the composite full-length (no omissions in structure) versus the VUS. Monte Carlo simulations were performed on the mutant to allow local regional changes for full-length 871 amino acids and when the p.K801E variant was introduced.
Using Monte Carlo, we were able to refine the multitude of X-ray structures (all partially missing fragments, side chains, loops, or rotamers), to generate a composite fulllength model for studies with simulations using our previous work with the Yasara SSP/PSSM method (Altschul et al., 1997;Hooft, Sander, Scharf, & Vriend, 1996;Hooft, Vriend, Sander, & Abola, 1996;King & Sternberg, 1996;Krieger et al., 2009;Qiu & Elber, 2006). The structure was relaxed to the YASARA/Amber force field using knowledge-based potentials within Yasara. The side chains and rotamers were adjusted with knowledge-based potentials, simulated annealing with explicit solvent, and small equilibration simulations using multiple parallel Yasara refinement protocols (Laskowski, Macarthur, Moss, & Thornton, 1993). Both the entire full-length structure and the VUS structure were modeled, filling in any gaps, or unresolved portions from the X-ray.
Refinement of the finalized model using either Schrodinger's LC-MOD Monte Carlo-based module or NAMD2 protocols was completed. All refinements started with YASARA generated initial refinement for variant p.K801E or the wild-type TRPV4 protein (Altschul et al., 1997;Krieger et al., 2009). The superposition and subsequent refinement of the overlapping regions yield a complete model for protein p.K801E or the wild type. The final structures were subjected to energy optimization with PR conjugate gradient with an R-dependent dielectric to ensure that the stereochemistry was in strict adherence to the force field.
Monte Carlo dynamics searching was completed on each model for conformational sampling, using methods previously described in the literature Caulfield & Devkota, 2012;Caulfield, Devkota, & Rollins, 2011;. Briefly, each TRPV4 or TRPV4 VUS was minimized with relaxed restraints using either Steepest Descent or Conjugate Gradient PR and then allowed to undergo the MC search criteria, as previously demonstrated Caulfield & Devkota, 2012;. The primary purpose of MC, in this scenario, is examining any conformational variability that may occur with different mutation, structural change, in or about, the region near to the site and possible effect on partner proteins complementary binding disruption with TRPV4 (Table 1, dysplasias and characteristics).

| RESULTS
For WT versus the variant p.K801E, we found the stability of the object from energetic calculations for G indicative of object stability shifts that correspond to worsening of the structural integrity from that of the wild type. In particular, for the tetramer, we find the change in object stability (folding state) and ability to partner (form tetramers). Stability of the object for TRPV4 (all 871 amino acids), as measured by G, is 649.73 kcal/mol*å 2 (162.43 per monomer), and per residue is 0.186 kcal/aa*mol*Å 2 . Whereas the p.K801E variant has a total stability G of 788.92 kcal/mol*Å 2 , per residue of 0.23 kcal/aa*mol*Å 2 shows an increase of nearly 24% in the per-residue energy for worsening of stability in the VUS over that of the wild type. If we look instead at direct interaction energy (favorable) for the tetramer complex instead of the individual object stability, we also see an interesting decrease in tetramerization favorability; where the wild-type, full-length tetramer has an G (interaction) of −120.13 kcal/ mol*Å 2 that is constructive, versus the p.K801E variant, which only has −59.81 kcal/mol*Å 2 indicating less robust inter-monomer binding. We have used this method in the past hypothesis generation and mutagenesis experiments Caulfield & Devkota, 2012;López-Vallejo et al., 2011;Reumers et al., 2005;Schymkowitz et al., 2005;Zhang et al., 2013).
This object stability designates changes in structure deleterious to configuration of the complex from proximate inspection (Figures 4 and 5). Additionally, we examined the local residues that may be useful to explain any change in function ( Figure 6). The molecular model for the full structure and its variant form were analyzed using our established techniques (Abdul-Hay et al., 2013;Ando et al., 2017;Caulfield & Devkota, 2012;Caulfield et al., , 2014Fiesel, Ando, et al., 2015;Fiesel, Caufield, et al., 2015;López-Vallejo et al., 2011;Puschmann et al., 2017;Zhang et al., 2013). The region around K801 is critical for various loops essential to the interface of the tetramer. Local residues within the 12Å cutoff near the variant site are critical to formation of inter-monomer interactions that form the pore complex for this cation channel protein, (Figures 4 and 5a,b). Additionally, we can see from the mapping of the surface, plus the loop changes in structure, and resulting effect on the tetramer stability, that the predicted tetramer of TRPV4 for the variant has a substantial loss of the interface compared to the wild type to the p.K801E variant and also the molecular surface shown in the zoomed-into insets in the figure panel give revealing loss of pore structure that holds the complex (Figures 5 and 6).
The p.K801E tetramer still forms a complex but is not as tightly held together and the pore becomes more disordered and has less regular structure than the wild-type protein (Figure 5c), which accounts for lowering its interaction potential with partner proteins and the loss of attractive G Caulfield & Devkota, 2012;López-Vallejo et al., 2011;Reumers et al., 2005;Schymkowitz et al., 2005;Zhang et al., 2013). Further use of protein-level knowledge as a diagnostic tool for clinician continues to provide new insights and help aid in the development of a database Harris et al., 2018;Richter et al., 2018;von Roemeling et al., 2018).

| DISCUSSION
The expression of TPRV4 in a multitude of specialized bone cells is the reason why such a large number of symptoms can result from a single mutation. Overproduction of follistatin, which interferes with endochondral ossification and is thought to cause skeletal dysplasia, only explains part of the proband's phenotype . She initially presented with bone pain and osteoporosis/osteonecrosis in addition to her skeletal abnormalities; her skeletal examination revealed widespread arthropathy as a source of her "constant" joint pain. Existing literature on TRPV4′s purpose in osteoclasts and cartilage helps explain these symptoms (Table 1).
Osteoclast differentiation is controlled by a Ca 2+ signaling pathway, where higher intracellular concentrations of the ion ultimately promote osteoclast-specific gene activity. As a regulator of Ca 2+ influx, TRPV4 can effectively control the number of osteoclasts produced in mammals. Pathogenic gain-of-function variants of TRPV4 increase Ca 2+ uptake, which ultimately raises the number of osteoclasts in an individual. This leads to higher levels of bone resorption than production, a resulting reduction in bone density, and even osteoporosis in some individuals Masuyama et al., 2012). Conversely, TRPV4 deficiency is known to abate bone loss usually caused by mechanical unloading, a phenomenon that occurs in inactive or bedridden individuals (Mizoguchi et al., 2008). This observation would suggest that a loss-of-function variant of TRPV4 is not behind the proband's osteopenia. Taken with the hypothesis that increased uptake of Ca 2+ is a generator of skeletal dysplasia, we can likely assume that NC_000012.12:c.2401A>G (p.K801E) is a gain-of-function variant.
The arthropathy observed with the proband challenges these above observations. Her arthropathy most closely resembles osteoarthritis, a phenotype that is associated with T A B L E 1 (Continued) pathogenic variants of TRPV4 in both experimental mouse models and human cases (Clark, Votta, Kumar, Liedtke, & Guilak, 2010;Lamandé et al., 2011). These results were seen with knockdown of TRPV4 expression (mice) or loss-of-function TRPV4 pathogenic variants (humans). Conversely, the skeletal dysplasia and osteoporosis also seen in our proband are known to be associated with gain-of-function TRPV4 variants Masuyama et al., 2012). This contrast suggests that further mechanisms maybe implicated in these manifestations as opposed to simply loss or gain of function. Due to this conflict of symptoms, we cannot say with confidence whether this TRPV4 variant causes loss of function or gain of function in TRPV4 protein.
Despite the potential TRPV4 variant conflict between the clinical features of the proband, they have all been seen before in TRPV4-associated disorders. We consider the proband TRPV4 structure that shows four monomer proteins that comprise the cation channel, as depicted in ribbons. Colors indicate the composite technique for our hybrid model technology to generate highly reliable structural models for VUS determination. This structural model is complete full-length sequence including residues and loops missing from X-ray structures. (b) Rotated view to look down into the channel pore, called the "bird's eye" view. The labeling for both models is indicated. Residues shown are rendered in either "licorice sticks" or van Der Waals (VdW) and using standard element coloring (O-red, N-blue, H-white, S-yellow) except for the carbon atoms that are colored to match the ribbon color a b to fit most closely with the spondylometaphyseal dysplasia, Kozlowski type phenotype, though there is some overlap with spondyloepiphyseal dysplasia, Maroteaux type. Both disorders include brachydactyly, platyspondyly, and metaphyseal abnormalities in their phenotypes, signs that were all seen in the proband. However, she does not have bowlegs or a waddling gait, features one might expect to see in these disorders (Schindler et al., 1993). Both also implicate some degree of short-trunk short-stature dwarfism, and while our proband stands at only 5′2″, she is taller than the expected heights for most patients having either disease. Autosomal dominant brachyolmia, a milder TRPV4-associated disorder, actually has the closest fit for the proband's height, with an expected adult range between 160 and 168 cm (or between 5′3″ and 5′6″) . The premature osteoarthritis seen in our proband is suggestive of spondylometaphyseal dysplasia, Kozlowski type, as is her childhood scoliosis. Conversely, her osteoporosis is known to be a symptom of spondyloepiphyseal dysplasia, Maroteaux type Schindler et al., 1993). The indications from the X-ray support the other findings with regard to the bone osteoporosis and dysplasia observed. The symptomatology of our proband cannot be clearly defined by a single TRPV4associated disorder; it seems that her case is comparable to a mosaic of 2-3 different disease phenotypes across the TRPV4-spectrum. Both the proband and her daughter tested positive for the c.2401A>G (p.K801E) variant of TRPV4, and while we were unable to test her father, it seems all but certain that he gave the proband this variant. The father's short digits, abnormal nails, and hip abnormalities were all phenotypic features seen in the proband; of these features, only short digits were also seen in her daughter. The father was also of mildly short stature, standing at 5′7″. The cases of this man's daughter and granddaughter combined with his phenotype suggest that he too had skeletal dysplasia ( Figure 3). If his hip issues were similar to those seen in the proband, it is possible that he had some form of arthropathy as well. When we investigate the confirmed family members with c.2401A>G (p.K801E), there is some overlap and discrepancy in phenotype. Both the proband and her daughter had scoliosis in early life that required correction, though it would appear the daughter's case was more severe as it required surgical action. It would be important for the daughter to undergo osteopenia evaluation, given a history of femur fracture at age 6. Bone and joint pain are exclusive to the proband, though her daughter might develop these symptoms later in life. The daughter also presents with macrocephaly, a feature that while not seen in the proband is not altogether unexpected given the variability of TRPV4 genotype-phenotype relationships .

| CONCLUSIONS
Our observation that a loss-of-function variant of TRPV4 is not behind the proband's osteopenia is determined and F I G U R E 6 Close-up views of the TRPV4 variant (p.K801E) surface and variant effect on loop structure and effect on tetrameric interface.
(a) Full-length model for the variant TRPV4 tetrameric structure colored by secondary structure and atom type. Interacting residues between the two proteins are shown. (b) The tetrameric structure for wild type (K801) is shown with same scheme. (c) Close-up view of the loops is shown for the wild type in this monomer-monomer interaction region. (d) Same view for the p.K801E is shown to illustrate the changes in the loop conformation and the effect on tetramer conformation. This VUS causes the pore structure to be altered for the cation channel, which also demonstrates changes in the protein Gibbs-free energy for stability a b c d supports the hypothesis that increased uptake of Ca 2+ is a generator of skeletal dysplasia, which taken together with the protein informatics and modeling, we find it plausibly likely that NC_000012.12:c.2401A>G (p.K801E) is a gainof-function variant of TRPV4. The structural models of these two variants concur with previous clinical data. TRPV4 is added into our repository for the database of all protein-level VUSs for future use to the clinical community to be made available as an online publicly available resource. Here, we demonstrate the fusion of protein informatics and modeling as another tool to complement the difficult task of variant classification.