Abstract
To maximize the impact of precision medicine approaches, it is critical to identify genetic variants underlying disease and to accurately quantify their functional effects. A gene exemplifying the challenge of variant interpretation is the von Hippel–Lindautumor suppressor (VHL). VHL encodes an E3 ubiquitin ligase that regulates the cellular response to hypoxia. Germline pathogenic variants in VHL predispose patients to tumors including clear cell renal cell carcinoma (ccRCC) and pheochromocytoma, and somatic VHL mutations are frequently observed in sporadic renal cancer. Here we optimize and apply saturation genome editing to assay nearly all possible single-nucleotide variants (SNVs) across VHL’s coding sequence. To delineate mechanisms, we quantify mRNA dosage effects and compare functional effects in isogenic cell lines. Function scores for 2,268 VHL SNVs identify a core set of pathogenic alleles driving ccRCC with perfect accuracy, inform differential risk across tumor types and reveal new mechanisms by which variants impact function. These results have immediate utility for classifying VHL variants encountered clinically and illustrate how precise functional measurements can resolve pleiotropic and dosage-dependent genotype–phenotype relationships across complete genes.
Main
Delineating rare genetic variants underlying disease phenotypes remains a major challenge in human genetics. For the majority of cancer-associated genes, more variants of uncertain significance (VUS) have been reported than variants whose phenotypic effects are known1,2,3,4. In the context of both germline testing and tumor profiling, VUS represent a missed opportunity to improve patient care through precision medicine approaches.
Most variants observed are too rare to enable statistically robust genotype–phenotype associations. Computational models of variant effects have improved due to greater availability of training data and the use of machine learning5,6,7,8,9,10. However, such models are not accurate enough to dictate clinical decisions without additional evidence11. Mechanistic knowledge of variants in tumor suppressor genes can inform which individuals will benefit from preventative measures and guide therapeutic selection12,13,14. There is, however, a scarcity of functional data available for linking variants to phenotypes15.
The von Hippel–Lindau (VHL) tumor suppressor is a 213-amino acid protein encoded on chromosome 3p that functions as an E3 ubiquitin ligase in complex with cullin-2, elongins C and B (ELOC and ELOB) and ring-box 1 (ref. 16). In normoxic conditions, VHL ubiquitinates the ɑ-subunit of hypoxia-inducible factor (HIF), targeting HIF for proteasomal degradation. In hypoxic conditions, HIF is protected from VHL-mediated degradation and signals to promote glycolysis and angiogenesis. Loss of VHL function due to mutation can lead to constitutive HIF activity and tumor development17.
Somatic VHL mutations are frequently observed in renal cell carcinomas (RCCs), most commonly clear cell RCC (ccRCC). During ccRCC evolution, chromosome 3p deletion typically precedes a loss-of-function (LoF) mutation to the remaining VHL allele, resulting in increased HIF activity18. VHL mutations have been observed in other types of RCC and extrarenal cancers, but their functional significance is less certain.
Pathogenic germline variants in VHL predispose patients to different neoplasias in an autosomal dominant manner, a condition known as VHL disease19. Affected patients have varying susceptibilities to different tumors, including ccRCC, pheochromocytoma and hemangioblastoma. The risk of each tumor depends largely on the specific germline variant. Classically, type 1 VHL disease variants lead to complete LoF and include whole-gene deletions, nonsense variants and frameshifting insertions and deletions (InDels). Type 1 variants predispose patients to ccRCC, hemangioblastoma and other neoplasms in a HIF-dependent manner. Type 2 variants, in contrast, are associated with high pheochromocytoma risk and are most often missense variants. Attempts have been made to subclassify type 2 VHL disease further—type 2C disease is marked by pheochromocytomas only, type 2A disease includes hemangioblastomas and other benign tumors and type 2B disease further includes ccRCC20.
These clinical classifications have helped explain patterns of tumors in families, yet a complete molecular accounting of how different mutations confer distinct pathologies has remained elusive. In a curated database of VHL mutations21, many variants have been associated with both type 1 and type 2 disease. Other variants have been implicated in recessive diseases, such as congenital polycythemia or germline VHL deficiency22,23. Although protein-truncating variants typically cause type 1 disease, in rare instances, patients with nonsense variants have presented with type 2 disease24. Such discrepancies highlight the challenge of developing individualized surveillance and therapy plans for patients without functional evidence of variants’ effects.
In addition to variants known to cause VHL disease, there are over 800 VUS in VHL reported in ClinVar3. It is unknown what fraction of these variants cause disease, and likewise how many variants yet to be reported may prove pathogenic. While genetic evidence may be converging on a near-complete set of pathogenic VHL alleles21, substantially more variants are predicted by computational models to be deleterious than have been linked to disease25. It is unclear whether such variants have yet to be encountered due to their rarity, whether they are incompatible with life or whether they are truly benign.
For patients harboring VHL variants whose phenotypic effects are unknown, well-calibrated functional data may prove useful in aiding diagnosis and management. The recently demonstrated efficacy of HIF2ɑ inhibitors for preventing ccRCC progression14 suggests quantifying variants’ precise effects on HIF regulation may prove valuable for guiding therapeutic selection. More broadly, VHL serves as a powerful model to assess the extent to which functional data can recapitulate genotype–phenotype relationships in humans, owing to the fact that it is frequently mutated in ccRCC2 and the extensive knowledge regarding phenotypic effects of germline variants21.
Here we systematically measure the functional consequences of VHL variants across the complete gene by using saturation genome editing (SGE). In total, we scored 2,268 single-nucleotide variants (SNVs) for HIF-dependent effects on cellular fitness, defining LoF variants underlying ccRCC development with 100% accuracy. Our assay captures clinically meaningful differences in the degree of functional impairment among pathogenic alleles and reveals new mechanisms explaining genotype–phenotype associations, suggesting a role in improving diagnostic and therapeutic precision.
Results
An optimized SGE assay to precisely score VHL variants
To develop a high-throughput assay for VHL variants, we assessed the effect of VHL loss across cell lines. Except for kidney-derived lines, CRISPR-induced knockout of VHL almost uniformly reduces cell fitness (Fig. 1a)26,27. To investigate VHL’s essentiality in the haploid human line HAP1, InDels in exon 2 were generated with CRISPR and sequenced at multiple timepoints. Robust depletion of InDels over time confirmed the essentiality of VHL for normal HAP1 proliferation (Fig. 1b). The strong selection against InDels was eliminated by prior knockout of HIF1A (Fig. 1b), indicating VHL loss confers a HIF-dependent growth defect.
SGE is a method by which all possible SNVs in a genomic region are assayed in multiplex using CRISPR–Cas9 editing28. When SGE is performed in HAP1, a single variant is engineered per haploid cell, allowing variants’ effects on growth to be quantified by next-generation sequencing (NGS). Seven SGE libraries were made to tile the coding sequence of VHL, as well as exon-proximal regions of introns and a region deep within intron 1 covering the 5′ end of a reported pseudoexon (Fig. 1c)29,30. Each library consisted of all possible SNVs in a region cloned into vectors with homology arms to facilitate genomic integration (Fig. 1d).
To measure more subtle growth effects, SGE experiments were performed using a highly optimized protocol modified from published work31 to feature improved transfection efficiency, a longer time course and addition of 10-deacetyl-baccatin-III (DAB) to maintain haploidy32 (Fig. 1e,f, Extended Data Figs. 1 and 2 and Supplementary Note). Amplicon sequencing of genomic DNA (gDNA) collected on days 6 and 20 was used to calculate a ‘function score’ for each SNV reflective of cellular fitness (Methods). SNVs with significantly reduced function scores (that is, ‘depleted’ SNVs) were defined for each SGE region by applying a false discovery rate (FDR) of 0.01. After stringent quality filtering, function scores for n = 2,268 SNVs were obtained, comprising 85.4% of all possible SNVs in SGE regions (Supplementary Table 1).
Mapping LoF variants
The majority of variants designed but not scored map to the GC-rich 5′ region of exon 1, where the rate of editing was lowest (Extended Data Fig. 3a). In contrast to other coding regions, the 5′ region of exon 1 lacked significantly depleted SNVs, including nonsense variants (Extended Data Fig. 3b,c). This is consistent with a previously characterized alternative translation initiation site at p.M54 producing a fully functional VHL isoform33, and corroborates the lack of pathogenic variants proximal to p.M54 in ClinVar (Fig. 1c).
Moving forward, we restricted analyses to n = 2,200 SNVs scored with high reproducibility by excluding SNVs assayed in the 5′-most SGE region (Fig. 2a,b). Among n = 115 remaining SNVs upstream of p.M54, none were significantly depleted. Likewise, no SNVs assayed in the 3′-untranslated region (UTR) or deep within intron 1 scored significantly (Fig. 2b,c). In contrast, between p.M54 and p.R200, all but one nonsense variant was depleted (n = 43, median score = −2.4), as were all canonical splice site SNVs (n = 24, median score = −2.3; Fig. 2b,c). Most missense variants scored neutrally, although 22.4% were depleted.
To better resolve mechanisms of functional impairment, we derived n = 1,626 ‘RNA scores’ for coding variants by performing targeted RNA-sequencing on day 6 and day 20 samples (Extended Data Fig. 4 and Supplementary Table 1). RNA scores reflect SNVs’ effects on full-length VHL mRNA levels.
Comparison of RNA scores to function scores reveals only large reductions in mRNA confidently predict LoF (Fig. 3a). Indeed, 16 of 17 SNVs with RNA scores below −3.0 were significantly depleted, reflecting a minimum mRNA dosage required for normal growth (Extended Data Fig. 4c). While RNA scores across timepoints were highly correlated, variants depleted in mRNA on day 6 tended to be less depleted on day 20 (Extended Data Fig. 5a–e), suggesting selection for cells expressing sufficient VHL mRNA. Indeed, only 6 of 17 variants with day 6 RNA scores below −3.0 had day 20 RNA scores below −3.0….
