Interpreting ciliopathy-associated missense variants of uncertain significance (VUS) in an animal model

Purpose Better methods are required to interpret the pathogenicity of disease-associated variants of uncertain significance (VUS), which cannot be actioned clinically. In this study, we explore the use of a tractable animal model (Caenorhabditis elegans) for in vivo interpretation of missense VUS alleles of TMEM67, a cilia gene associated with ciliopathies. Methods CRISPR/Cas9 gene editing was used to generate homozygous knock-in C. elegans worm strains carrying TMEM67 patient variants engineered into the orthologous gene (mks-3). Quantitative phenotypic assays of sensory cilia structure and function measured if the variants affect mks-3 gene function. Results from worms were validated by a genetic complementation assay in a human TMEM67 knock-out hTERT-RPE1 cell line that tests a TMEM67 signaling function. Results Assays in C. elegans accurately distinguished between known benign (Asp359Glu, Thr360Ala) and known pathogenic (Glu361Ter, Gln376Pro) variants. Analysis of eight missense VUS generated evidence that three are benign (Cys173Arg, Thr176Ile, Gly979Arg) and five are pathogenic (Cys170Tyr, His782Arg, Gly786Glu, His790Arg, Ser961Tyr). Conclusion Efficient genome editing and quantitative functional assays in C. elegans make it a tractable in vivo animal model that allows rapid, cost-effective interpretation of ciliopathy-associated missense VUS alleles.


INTRODUCTION
Exome and genome sequencing have revolutionised our ability to identify the genetic causes of disease. Missense variants (single codon altered to encode a different amino acid) are the most numerous class of protein-altering variants 1 but only a subset are associated with disease 2 . Based on disease features and patterns of inheritance, identified variants are classified as benign, likely benign, uncertain significance, likely pathogenic or pathogenic (as defined by the Association for Clinical Genomic Science; ACGS) 3 . For novel or previously uncharacterized variants, the only evidence available to assess their pathogenicity is population allele frequency and analysis by in silico tools (e.g. SIFT/PolyPhen/CADD ), which are not sufficient to meet the threshold for a 'likely pathogenic' classification according to best practices established by ACGS 3

. A VUS (variant of uncertain significance)
classification is made when there is insufficient evidence to conclude on pathogenicity 3,4 .
Currently, most missense variants are classified as VUS (116,215/198,080=58.7%, accessed from ClinVar 5 September 2021). Since VUS classifications cannot be used for molecular diagnosis of disease, they cannot be acted upon clinically. Therefore, a VUS classification can delay or prohibit accurate disease management or genetic counselling, and prevents patients from accessing gene-specific therapies and clinical trials 6,7 . Given the wealth of genomic data available, there is a pressing clinical need to definitively reclassify VUS as benign or pathogenic. Effective experimental strategies for the functional interpretation of VUS are required to address this problem 8 . With the emergence of advanced genetics tools such as CRISPR-Cas9 gene editing, non-rodent model organisms such as zebrafish, Drosophila and C. elegans are emerging as robust in vivo experimental platforms for determining variant pathogenicity [9][10][11] .
Ciliopathies are a heterogenous group of at least 25 inherited disorders with clinically overlapping phenotypes, caused by pathogenic variants in almost 200 genes 12 . Ciliopathies are caused by defects in cilia, which are microtubule-based organelles, typically 2-20 microns long, that extend from the surfaces of most cell types. Motile cilia propel cells through a fluid or push fluid across a tissue surface. Primary cilia act as cellular "antennae" 13 , transducing a wide variety of extrinsic chemical and physical (eg. light, odorants) signals into the cell 14 .
Primary cilia are especially important for coordinating cell-cell communication signaling pathways (eg. Shh, Wnt, PDGF-a) essential for development and homeostasis 15 . Ciliopathies affect many organ systems, causing a broad range of clinical phenotypes of varied severity and penetrance that include cystic kidneys, retinal dystrophy, bone abnormalities, organ laterality defects, respiratory tract defects, infertility, obesity, neurodevelopmental defects and cognitive impairment 16 .
TMEM67 is a transmembrane protein that localizes specifically to the proximal 0.2-1.0 m µ region of the ciliary axoneme called the transition zone 17,18  elegans is a tractable model system that can provide evidence of variant pathogenicity.

C. elegans maintenance
All C. elegans strains in this study were maintained at 20ºC or 15ºC on nematode growth medium seeded with OP50 E. coli using standard techniques 22 . All worm strains are listed in Table S2.

CRISPR/Cas9 to engineer mks-3 mutants in C. elegans
CRISPR protocols were performed as previously described 23 in nphp-4(tm925) using an unc-58 co-CRISPR strategy. crRNA are listed in Table S3. The CRISPR efficiency (defined as the percent of F1 pools that were positive for the edit by PCR) varied from 1-35% with an average 15%. Accuracy of the engineered variants was confirmed with Sanger sequencing.
unc-58 was also sequenced and unintended unc-58 mutations 24 were outcrossed. PCR and sequencing primers are listed in Table S4.

C. elegans quantitative phenotyping assays
Dye filling, roaming, and chemotaxis assays were performed as previously described 25

C. elegans wide-field imaging and quantification of fluorescence
Young adult hermaphrodites were immobilised on 4% agarose pads in 40 mM tetramisole (Sigma). Images were acquired with a 100x (1.40 NA) oil objective on an upright Leica DM5000B epifluorescence microscope and captured with an Andor iXon+ camera. Image analysis was performed with FIJI/ImageJ (NIH). MKS-3::GFP fluorescence was quantified as previously described 23 .

Transmission electron microscopy
Young adult hermaphrodites were processed as previously described 25

PCR and sequence validation of crispant cell-lines
To extract DNA from colonies within the 96-well plates, cells were washed with 1x

Whole cell extract preparation and western immunoblotting
Whole cell extracts containing total soluble proteins were prepared from hTERT-RPE1 cells and visualized using a ChemiDoc MP imaging system (BioRad Inc., Hercules, CA, USA).
Ratios of active phosphorylated ROR2 : unphosphorylated ROR2 isoforms were calculated by quantitating band intensity using ImageLab 5.2.1 software (BioRad Inc.) for three biological replicates, as described previously 26 .

Statistical analyses
All

Selection of TMEM67 variants for analysis
TMEM67 variants were selected using two criteria: (i) conservation of the mutated amino acid in the worm orthologue (MKS-3) ( Figure S1 ), and (ii) presence of an adjacent Cas9 PAM site in the C. elegans genome to facilitate CRISPR gene editing. Using these criteria, we modelled eight missense TMEM67 VUS in C. elegans mks-3 ( Figure 1A, Table S1 ). We also included two known benign and two pathogenic variants as controls.

Effect of VUSs on transition zone localization of MKS-3 and ultrastructure
To provide further insight into the damaging or benign nature of the TMEM67 VUS alleles, we examined the effect of the variants on the subcellular localisation of MKS-3. In C. elegans sensory neurons, transmembrane MKS-3 localizes to the ciliary transition zone (TZ), which corresponds to the most proximal 1 μm of the ciliary axoneme 18 . Fusion PCR was used to generate linear mks-3::gfp fragments ( Figure S3A ) containing variants of interest.
We also assessed the effect of the TMEM67 VUS variants on the ultrastructure of cilia and TZs using transmission electron microscopy (TEM). Specifically, we analyzed two alleles: VUS1 (predicted pathogenic) and VUS2 (predicted benign). Pathogenic2 was also examined as a control. Wild-type amphidal pores contain 10 rod-shaped ciliary axonemes emanating  ( Figure 3D ). The VUS2-containing strain displays a similar TZ phenotype to the positive control, whereas the VUS1-containing strain shows enhanced TZ membrane-MT detachment that is similar to the mks-3(∆); nphp-4(∆) negative control and Pathogenic2 ( Figure 3D ). Altogether, the ultrastructure data confirms the pathogenic and benign nature of VUS1 and VUS2, respectively.

In vitro genetic complementation assay of TMEM67 VUS function in human cell culture
To further validate our findings from C. elegans , we utilized an in vitro human cell culture-based assay of TMEM67 function. Previously, we demonstrated that TMEM67 is required for phosphorylation of the ROR2 co-receptor and subsequent activation of non-canonical Wnt signalling 26 ( Figure 4A ). Here, we developed a genetic complementation approach using a validated hTERT-RPE1 crispant cell-line that has biallelic null mutations in TMEM67 ( Figure S4 ). In the absence of TMEM67, phosphorylation of ROR2 was not stimulated by exogenous treatment with the non-canonical ligand Wnt5a. Transient transfection with full-length wild-type TMEM67 fully rescued ROR2 phosphorylation following Wnt5a treatment. These responses allowed us to develop a series of assays to determine the relative effects of VUS on biological function. In this assay, transfection of Benign1 allowed 204.7% induction of phospho-ROR2 levels by Wnt5a relative to control ( Figure 4B ). In contrast, Pathogenic2 did not rescue biological function (90.2% induction) ( Figure 4B ). Comparison of all VUS, normalized to wild-type TMEM67 responses across three independent biological replicates, enabled us to interpret VUS1, VUS4, VUS5, VUS6 and VUS8 as pathogenic, and VUS3 and VUS7 as benign ( Figure 4C ). VUS2 was not tested in this assay.

DISCUSSION
Ciliopathies are multisystem disorders that affect many organs including kidneys, liver, and retina. While the organ systems affected by cilia dysfunction are not present in C. elegans , the basic biology of primary cilia is conserved. Despite undoubted context-specific distinctions, cilia proteins are functionally conserved 31,32 , therefore allowing us to model missense variants of these proteins in worms.
In this study, we exploited efficient genome editing and quantitative phenotypic analysis in C. elegans to determine the pathogenicity of TMEM67 variants. This approach accurately classified known pathogenic and known benign variants. We also generated a pathogenicity prediction for all eight missense VUS alleles analyzed. Three VUS were phenotypically benign (VUS2(Cys173Arg), VUS3(Thr176Ile),VUS7(Gly979Arg)) and five were phenotypically pathogenic (VUS1(Cys170Tyr), VUS4(His782Arg), VUS5(Gly786Glu), VUS6(His790Arg), VUS8(Ser961Tyr)). These pathogenic predictions were validated in a human cell culture assay of TMEM67 function. We propose that experiments in C. elegans can interpret the pathogenicity of VUS and provide evidence towards their reclassification as benign or pathogenic.
Several in silico algorithms have been developed to predict the pathogenicity of missense variants. However, their accuracy is inconsistent [33][34][35] . Indeed, we tested five commonly used in silico tools (MISTIC, SIFT, Poly-Phen, CADD, and REVEL) and found that they return deleterious/damaging predictions for all 8 of the VUS examined in this project ( Figure S5 ).
The only exceptions are VUS1/3/8, where 1 or 2 of the tools returned non-deleterious or intermediate scores ( Figure S5 ). Therefore, the algorithm predictions do not align with our observations in C. elegans and cell culture experiments. We conclude that in vivo modelling of missense variants in C. elegans more accurately predicts results in human cells than currently available prediction algorithms.
The quantitative assays employed in this study are suitable for high throughput analysis. Live animal fluorescence-activated cell sorting can be used for high throughput dye filling assays 36 and automated worm tracking can be used for high throughput roaming and chemotaxis assays 37 . Additionally, machine learning can streamline the analysis of complex datasets to predict VUS pathogenicity 38 . One limitation to modelling patient variants in endogenous C. elegans genes is the conservation of amino acid residues. Although many cilia genes have clear orthologs in C. elegans, the amino acid conservation varies from 30-60%. C. elegans genes can be "humanized" to remove this limitation [38][39][40] . In summary, this study highlights that C. elegans is a practical model for variant interpretation of ciliary genes. Analysis of ciliopathy-associated VUS in C. elegans is accurate, quick, affordable, and easily interpretable. While this study focused on TMEM67, we anticipate that VUS alleles of any conserved ciliary genes can be modelled and characterized in C. elegans .
The authors declare no competing interests.

Figure 1. TMEM67 variants analyzed in this study
A) Twelve variants in TMEM67/mks-3 were generated by CRISPR-Cas9 gene editing and characterized in C. elegans . The schematic shows the relative positions of the variants along the length of the proteins. The domains are conserved between humans and worms, but the worm protein lacks the N-terminal signal peptide. B) RaptorX protein structure and domain organization predictions for the full-length TMEM67 and MKS-3 proteins. RaptorX is a deep learning algorithm that predicts secondary and tertiary structures of proteins without close homologues or known structures in the protein data bank. Ribbon diagrams of proteins are rainbow-coloured (red at N-terminus to dark blue at C-terminus) with variants indicated (magenta -known pathogenic; yellow -VUS; green -known benign) within the predicted protein domains. Note that transmembrane helices 1-5 (green to cyan) are followed by a single coiled-coil domain (light blue) and then two C-terminal transmembrane helices 6-7 (dark blue).

Figure 2. Quantitative phenotyping of cilia-dependent phenotypes in C. elegans
All assays were performed blind with at least three independent biological replicates. A) Lipophilic dye filling assay of the four phasmid (tail) neurons. The number of cell bodies which uptake dye was counted (values range from 0-4). The bar graph indicates the proportion of the population with dye uptake in 0 to 4 phasmid neurons. The number of worms is shown in brackets. Statistical significance according to Kruskal-Wallis followed by Schaich-Hammerle post hoc test. B) Assessment of worm roaming behaviour normalised to wild-type. A single young adult hermaphrodite was placed on a food-rich plate for 20 hours and the roaming activity was quantified. The number of worms is shown in brackets. Box plots indicate the maximum and minimum values (bars), median, lower quartile, and upper quartile. Statistical significance according to Kruskal-Wallis followed by Dunn's post hoc test. C) Quantification of worm chemotaxis towards benzaldehyde after 60 minutes. Assay is performed on a population of 50-300 worms. The number of assays is shown in brackets.Statistical significance according to ANOVA followed by Tukey's post hoc test. Box plots indicate the maximum and minimum values (bars), median, lower quartile, and upper quartile. D) Integration of the phenotypic results from the three quantitative assays into one value. Averages from each assay were normalized to the nphp-4 control (with a maximum score of 1.0 for each assay) and summed. Values were ranked from highest (benign) to lowest (pathogenic). C) The mean number of cilia observed at the mid-pore region. The number of pores is shown in brackets. D) TZ membrane-microtubule detachment phenotypes of differing severity: mild -TZ microtubules are slightly detached from the membrane, which is expanded to a small degree; moderate -many TZ microtubules are clearly detached from the membrane, which is extensively expanded on one side of the axoneme; severe -most or all of the TZ microtubules are fully detached from the membrane, which is highly expanded, indicating that the TZ is ectopically docked in the periciliary membrane compartment (PCMC) cytoplasm. The number of TZs is shown in brackets.