Skip to main content

A polyclonal allelic expression assay for detecting regulatory effects of transcript variants

Abstract

We present an assay to experimentally test the regulatory effects of genetic variants within transcripts using CRISPR/Cas9 followed by targeted sequencing. We applied the assay to 32 premature stop-gained variants across the genome and in two Mendelian disease genes, 33 putative causal variants of eQTLs, and 62 control variants in HEK293T cells, replicating a subset of variants in HeLa cells. We detected significant effects in the expected direction (in 60% of variants), demonstrating the ability of the assay to capture regulatory effects of eQTL variants and nonsense-mediated decay triggered by premature stop-gained variants. The results suggest a utility for validating transcript-level effects of genetic variants.

Background

The interpretation of the functional effects of common and rare variants in the human population is a major objective in human genetics and genomics. Despite the success of mapping genetic associations to complex traits by genome-wide association studies (GWAS) and interpreting their effects on gene expression by expression quantitative trait loci (eQTL) studies [1,2,3,4,5], the causal variants at GWAS loci and eQTLs are usually unknown due to linkage disequilibrium (LD). Statistical fine-mapping methods [6,7,8] can help narrow down causal variants, but experimental validation of the performance of these methods is lacking. For rare variants, functional interpretation has distinct challenges even in the well-annotated coding regions. Rare disease studies often result in hundreds to thousands of potential disease-causing variants identified from whole-exome sequencing, and prioritization based on their functional effect is essential for research and clinical use [9].

Thus, there is a need for experimental methods to confirm the effects of common and rare variants. Methods such as massively parallel reporter assays (MPRAs) [10, 11], which couple regulatory sequences with an expression-correlated reporter, are high-throughput approaches for finding active regulatory variants outside of the gene body, and analogous methods exist for variants affecting splicing [12, 13]. However, the results of the assays show low concordance with eQTL data [10, 11], perhaps due to taking the variant out of its genomic context. Furthermore, MPRAs are not suited to testing variants within the transcript that can affect gene expression levels via post-transcriptional mechanisms, e.g., RNA stability. Regulatory variants are strongly enriched not only in promoters and enhancers, but also in UTRs and other transcript annotations [1], emphasizing the need for a method to validate them. Similar mechanisms also apply to a subset of rare disease variants, where stop-gained and frameshift variants can affect transcript abundance via nonsense-mediated decay (NMD) [14]. Stop-gained variants located 50–55 bp or more before the last exon junction often induce NMD, while variants located beyond this threshold are more likely to escape NMD and therefore produce truncated protein [15]. However, these predictions are not perfect [16]. Variants that trigger or escape NMD in the same gene can manifest in diseases with different symptoms or methods of inheritance [17], making it important to validate whether a given variant induces NMD.

CRISPR/Cas9 genome editing technology [18,19,20] has provided a means to introduce specific variants into the genome in order to validate their effects on expression in the native genomic context. However, editing one variant at a time, isolating hundreds of single-cell clones, genotyping and expanding clones, and measuring transcript abundance are hugely time-consuming and expensive processes. In addition to the resource cost of completing such an experiment, undetected large on-target mutations [21], off-target mutations, and other clone-specific abnormalities can create noise which requires many replicates of each desired genotype. A less labor-intensive genome editing approach analyzes allelic expression in the polyclonal edited cell population and has been used to validate the effects of specific rare variants [22] and all possible mutations in a particular exon using saturation mutagenesis [23].

In this study, we decided to develop and apply a similar polyclonal approach for medium-throughput testing of the expression level effects of eQTLs in transcribed regions and rare premature stop variants. We first tested rare premature stop variants with signs of NMD in the Genotype-Tissue Expression (GTEx) v8 data to validate the ability of our assay to detect effects on transcript abundance, and then applied the assay to fine-mapped eQTLs from GTEx. Finally, we assayed premature stop variants in two Mendelian disease genes, GLI3 and ROR2, to evaluate our ability to test NMD in a clinically useful context. Stop-gained variants towards the beginning of GLI3 are associated with Greig cephalopolysyndactyly, while variants towards the end of the gene are associated with the clinically distinct Pallister-Hall syndrome [24, 25]. Stop-gained variants towards the beginning of ROR2 are associated with the autosomal recessive Robinow syndrome, while variants towards the end of the gene are associated with autosomal dominant Brachydactyly type B1 [26, 27]. It has been hypothesized that in both genes, the disease manifestation is impacted by whether or not the variant triggers NMD. In such situations, experimental testing of NMD can be valuable for disease diagnosis and prognosis.

Methods

fgwas enrichment

First, we sought to establish the relevance of testing eQTL effects driven by variants within transcripts by analyzing the extent of cis-eQTL enrichment in functional elements of the genome. We used GTEx v6 fibroblast eQTL data and a diverse set of annotations: Gene annotations were obtained from GENCODE [28], and regulatory annotations (CTCF-binding site, enhancer, open chromatin region, promoter, promoter-flanking region, and TF binding site) were obtained from the Ensembl regulatory build release 80 [29]. Additional annotations include CADD variant consequence scores [30], SPIDEX machine learning-based prediction of splicing effects [31], experimentally validated miRNA-binding sites from Tarbase [32], 3′ UTR regulatory elements [33], and RNA-binding protein sites from CLIPdb [34]. Significant fibroblast eQTLs were analyzed for enrichment in these functional annotations using fgwas [35], with each annotation tested separately.

Assay design

The design of the assay is illustrated in Fig. 1a. In order to validate transcript regulatory variants’ allelic effects on transcript abundance, we utilized CRISPR/Cas9 genome editing with a gRNA specific to the locus of the variant of interest and a single-stranded DNA (ssDNA) template containing the alternative allele for homology-directed repair (HDR). For each variant of interest, we transfected the gRNA and ssDNA template into a well of inducible Cas9 human embryonic kidney cells (HEK293T). After editing, cells were harvested for gDNA and mRNA, followed by amplicon sequencing of the locus of interest in each. A regulatory effect of the variant is detected as a difference in the ratio of the alternative allele between gDNA and mRNA. This effect size is calculated as the log ratio of the alternative allele in cDNA over the ratio of the alternative allele in gDNA: log2(cDNA alt/ref/gDNA alt/ref), or the allelic fold change (aFC). Even though many cell lines are not diploid at all loci, this assay is detecting allelic events in cis where each homologous chromosome functions independently, and thus, variant effects are expected to be robust regardless of ploidy.

Fig. 1
figure 1

Polyclonal allelic expression assay to detect the effects of regulatory variants. a Assay schematic. Inducible Cas9 HEK293T cells undergo homologous recombination after transfection with the gRNA and ssDNA template in order to introduce the alternative allele to the cells. HeLa cells without inducible Cas9 were transfected with a Cas9 plasmid. Editing is followed by targeted sequencing of gDNA and mRNA to detect the ratio of alt/ref alleles in the polyclonal population of cells. b Table with the number of each type of control and putative regulatory variant edited with the assay in HEK293T cells. c Homologous recombination rate versus standard deviation for variants replicated 2–3 times with assay in HEK293T cells. The vertical line shows the 0.4% HDR cutoff which was used to filter variants for subsequent analysis. d Scatter plot showing the reproducibility of the effect size detected by the polyclonal allelic expression assay for two replicate experiments editing the same variants in HEK293T cells

Variant selection

In this study, we edited five types of variants: GTEx stop-gained, GTEx eQTL, disease gene stop-gained, non-eQTL synonymous control, and synthetic control variants (Fig. 1b and Additional file 1: Table S1).

Stop-gained variants from the general population were obtained from the GTEx v6 data release. Starting with all stop-gained variants that were singletons in GTEx v6, we used allele-specific expression (ASE) data from the fibroblast sample of the individual carrying the variant to select those that are likely triggering NMD. The selected variants have RNA-seq coverage of ≥ 20 reads, a reference ratio Ref/(Ref + Alt) > 0.7, and are located in a gene with > 0.5 RPKM in a published HEK293 RNA-seq dataset [36]. In HEK293 cells, 12,451 genes (61%) are expressed at the level of 0.5 > RPKM. Additionally, we required ASE data in at least 5 tissues and a first quartile of ASE across tissues of > 0.7 to select variants where NMD does not appear to be highly tissue-specific. Finally, we selected variants > 30 bp from the end of an exon for primer design. Six variants were used for editing.

eQTL variants were obtained from the GTEx v8 data release. Significant eQTL variants in fibroblasts were filtered for being within at least one protein-coding transcript, having a CAVIAR fine-mapping posterior probability of association > 0.8, an eGene with > 1 RPKM in HEK293 cells, and an effect size in the top quartile of the effect sizes of all associations (aFC > 0.30). The top 33 highest effect size variants with successful gRNA and primer design were chosen for editing.

Ten stop-gained variants for each of the disease genes GLI3 and ROR2 were created by changing a codon in the transcript to a stop codon. The stop codons were spaced 20 bp apart in both directions from the NMD cutoff point (55 bp upstream of last the exon-exon junction). The 6 disease variants tested were obtained from ClinVar [37], choosing disease-associated variants in the two genes on either side of the NMD threshold.

We selected 30 non-eQTL negative control variants from common synonymous variants in GTEx v8 data with an eQTL association p > 0.1 with the gene in which they reside. The templates for the 32 synthetic control variants were designed by introducing a nucleotide other than the reference or alternative allele at the stop-gained variant locus, which does not create a premature stop codon.

Cell culture

Genome editing was carried out in a doxycycline-inducible Cas9 HEK293T cell line, transduced with pCW-Cas9 plasmid (Addgene plasmid #50661 [38]), courtesy of the Sagi Shapira Lab. HEK293T cells were cultured in OptiMEM (Gibco) supplemented with 5% HyClone Cosmic Calf Serum (Fisher), 1% GlutaMAX (Gibco), 1% NaPyr (Corning), and 1% penicillin/streptomycin (Corning). The cells were passaged and maintained following the standard techniques in 5% CO2 and 95% air.

For replication in HeLa cells, cells were cultured in DMEM (Corning) with 4.5 g/L glucose supplemented with 10% heat-inactivated FBS (Gibco), 1% GlutaMAX (Gibco), 1% NaPyr (Corning), and 1% penicillin/streptomycin (Corning).

Genome editing

The protocol for the polyclonal editing assay can be found at https://doi.org/10.17504/protocols.io.7c6hize. gRNAs were designed with E-CRISP version 5.3 [39] using medium settings, with an NGG PAM, a 5′ G; excluding designs with more than 5 off-targets; and classifying off-targets as having up to 3 mismatches in the 5′ region of the gRNA. gRNAs were ordered as gBlocks gene fragments (IDT): a U6 promoter sequence followed by the specific gRNA and tracr sequence [40]. The gBlocks were amplified using the Q5 High-Fidelity 2X Master Mix (NEB) and gBlock amplification primers [40]. Homologous templates were designed by extracting the sequence 50 bp upstream and downstream of each variant and substituting the reference allele with the alternative allele. Stop-gained control templates have another nucleotide substituted in the variant position which does not create a stop codon. Homologous templates were synthesized as ultramers by IDT. If possible, primers which amplify both cDNA and gDNA were designed using the IDT primer quest, choosing those that cover the PCR target (region spanning the variant and DSB) with at least 15 bp between the PCR target and one primer and at least 60 bp to the other primer. Otherwise, cDNA- and gDNA-specific primers were designed using either the cDNA or gDNA sequence as the template. Nextera adapter sequences were appended to forward and reverse primer sequences as follows:

GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG+ForwardPrimerSequence

TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG+ReversePrimerSequence

Primers were ordered as standard oligos from IDT.

Twenty-four hours before transfection for CRISPR editing, iCas9 HEK293T cells were plated in 24-well plates and induced with 5 μg/mL of doxycycline, with a separate well for each targeted variant. Cells were transfected with 500 ng homologous template and 500 ng gRNA gblock using the Lipofectamine MessengerMAX transfection reagent. After 24 h, the transfection reagent was removed and replaced with new media. Cells were split after 4 days and 6 days, and DNA and RNA were extracted from the polyclonal edited cultures at 9 days. Seventy-five percent of the 24-well culture was harvested for RNA using IBI Isolate DNA/RNA Reagent according to the manufacturer’s instructions. Purified RNA was quantified by Nanodrop (Thermo Fisher). cDNA was synthesized with ~ 200 ng of purified RNA using 1/4 reactions of SuperScript IV VILO Master Mix with EZ DNase (Invitrogen). Another 10% of the cell culture was used for DNA extraction using 15 μL of QuickExtract (Lucigen). For the timecourse optimization experiment, mRNA and gDNA were extracted as above at days 4, 6, and 9.

For editing selected variants (Additional file 2: Table S1) in HeLa cells, to confer puromycin and Cas9 expression, we used transient transfection with pCC_01 plasmid (Addgene #139086) [41]. Cells were plated in 24-well plates and transfected at 90% confluency with 500–750 ng homologous template, 500 ng gRNA gblock, and 400 ng pCC_01 vector using Lipofectamine 3000 transfection reagent. After 24 h, transfection reagent was removed and replaced with fresh medium plus puromycin (1,5μg/mL). Cells were kept under puromycin selection for 72 h, until non-transfected control cells were completely eliminated. The medium was replaced at days 4 and 6. mRNA and gDNA were extracted from the edited cultures at 9 days as described above.

Library preparation

Amplicon libraries from cDNA and gDNA were created using either the same Nextera primers (if possible) or separate Nextera primers for cDNA and gDNA. One microliter of cDNA or gDNA was amplified using Q5 High-Fidelity 2X Master Mix (NEB). An indexing PCR was performed next using Nextera XT index kit primers (Illumina) and NEBNext High-Fidelity 2X PCR Master Mix (NEB) which resulted in dual barcoded amplicons with Illumina adapters. cDNA and gDNA libraries were mixed in equal volume and sequenced on the MiSeq using 150-bp paired-end reads. We obtained a median coverage of about 85,000 reads per sample.

Sequencing analysis

Fastqs generated from Illumina software were trimmed for adapter sequences and quality using trimmomatic. Reads were aligned to the gDNA or cDNA sequence specific for each amplicon and categorized as HDR, no edit, or NHEJ using EdiTyper [42]. Variants were eliminated if HDR in gDNA was greater than 30% (suggesting the cell line is in fact heterozygous for the variant). Samples were filtered out if they had fewer than 1000 reads covering the locus of interest. Additionally, samples were filtered out if they had an outlier NHEJ rate of greater than 80%, indicative of an alignment error. Reference and alternative allele counts were obtained from the EdiTyper results (no edit and HDR reads, respectively) for each cDNA and gDNA sample (Additional file 2: Table S2). The effect size for each variant was calculated as the log2((Alt/Ref in cDNA)/(Alt/Ref in gDNA)), or allelic fold change (aFC). An effect size of zero means the variant has no effect on transcript abundance.

Statistical analysis

Significance between the control variant distribution and the other experimental variant types was determined using a two-sided Wilcoxon rank sum test. An F test was utilized to detect a difference in the variance of aFC between non-eQTL control and synthetic control variants, and eQTL and control variants. For each individual regulatory variant, a p value was calculated from the z-score of the variant’s effect size based on the mean and standard deviation of the control distribution. The p values were then Bonferroni corrected, and variants with a corrected p value of less than 0.05 were considered significant.

eQTL effect size in GTEx

For the GTEx effect sizes for the eQTLs, we used the allelic fold change (aFC) estimates from the GTEx v8 data release [4, 43, 44]. For the eQTL effect size in each GTEx tissue, we used the aFC estimates calculated from eQTL data. For stop-gained variants, we calculated aFC as the log2 ratio of the alternative and reference allele counts in the RNA-seq data in GTEx. To analyze the variation of eQTL variants’ effects across GTEx individuals, we calculated the aFC for each eQTL variant across heterozygous individuals in GTEx as the ratio of allele counts in the gene body [43]. Samples were filtered for those with greater than 50 reads covering heterozygous sites in the gene.

Results

First, we assessed the right time point to harvest mRNA after transfection with CRISPR constructs. Since mRNA is likely to remain in the cell for hours to days after editing has occurred, we expect to see fewer edited mRNA molecules early after transfection. To find the optimal time point in HEK293T cells, we edited 14 control variants (in 14 different genes) which are not expected to have an effect on expression and harvested at three time points post-transfection: 4 days, 6 days, and 9 days. At 4 days, the edited allele is depleted in the mRNA (Additional file 4: Figure S1a). However, this effect is lessened after 6 days and gone by 9 days. Therefore, we used the 9-day time point for the assay in order to analyze only the mRNA which has been transcribed post-editing.

Next, we assessed how technical variation in editing efficiency or PCR amplification may affect the robustness of the assay, using data from HEK293T cells. We analyzed how the homology-directed repair (HDR) rate affects the standard deviation of variant effect size between editing replicates of 49 variants (2 replicates for 23 and 3 replicates for 26 variants). We found that very low HDR is associated with a higher standard deviation between replicates (Fig. 1c). Therefore, in subsequent analyses, we discarded any variant with an HDR rate of less than 0.4% as determined by the frequency of the alternative allele in the gDNA. The HDR rate varies greatly between loci but is very well correlated between replicates of the same variant (Spearman’s rho = 0.96, p = 2.7 × 10−14, Additional file 4: Figure S1b), suggesting that the results of the assay are not strongly influenced by PCR amplification bias or variation in transfection efficiency. The HDR rate was not significantly correlated to NHEJ rate (Spearman’s rho = 0.15, p = 0.27) as has been observed before [45], and neither HDR nor NHEJ rate was significantly correlated to two published predictors of gRNA efficiency [39, 46] (Additional file 4: Figure S2). There was a minimal correlation between the gene expression level and the standard deviation of the effect size between replicates, which was reduced with the HDR filter of 0.4% (Additional file 4: Figure S1d). Altogether, the effect sizes of the two replicates are well correlated (Spearman’s rho = 0.60, p = 1.3 × 10−3, Fig. 1d). These results indicate that the polyclonal allelic expression assay can be used to robustly detect the regulatory effects of genetic variants even at relatively low HDR rates and in lowly expressed genes.

In order to determine the optimal set of negative control variants, we compared the distribution of effect sizes of the synthetic control variants (new variants created in the same genes as the stop-gained variants) and non-eQTL variants (common synonymous variants where eQTL effects were tested in GTEx and not observed) in HEK293T cells. The synthetic control variants have several outlier variants with large effect sizes that are consistent in replicates (Additional file 4: Figure S1c). This suggests that a subset of synthetic variants affect transcript levels and are thus not ideal for negative controls. The non-eQTL control variants, however, have effect sizes consistently close to zero (median aFC = − 0.009), demonstrating the utility of population data in selecting non-functional negative control variants. The variance of the synthetic controls was significantly greater than the variance of the non-eQTL controls (1.1 versus 0.038; F test p = 2.4 × 10−8). The non-eQTL variants were thus utilized as the control distribution for comparison with the stop-gained and eQTL variants tested with the assay.

In order to analyze the effects of genetic variants on gene expression levels using the polyclonal allelic expression assay, we first analyzed the rare stop-gained variants from GTEx that are expected to trigger NMD. As a group, the stop-gained variants show the expected negative effect sizes as compared to the control distribution (Wilcoxon p = 9.8 × 10−5, Fig. 2a). Four of these variants individually deviate significantly from the control distribution (Bonferroni-corrected z test p < 0.05, Fig. 2a and Additional file 3: Table S3). These results demonstrate our ability to capture NMD effects with the assay.

Fig. 2
figure 2

Stop-gained and eQTL variants from GTEx show allele-specific regulatory effects on expression. a Effect size of non-eQTL control, eQTL, and stop-gained variants after editing with the polyclonal allelic expression assay. Triangular points mark variants whose effect sizes significantly deviate from the control distribution. b Correlation between the effect sizes of variants in GTEx and effect sizes from the polyclonal allelic expression assay. c Correlation between the effect sizes of variants in HEK293T cells and in HeLa cells from the polyclonal allelic expression assay. d eQTL effect size (aFC) in GTEx tissues for the 13 edited eQTL variants shown as boxplots, with lines indicating the median effect size in GTEx fibroblasts and in the assay. Asterisks mark variants which were significant in the assay in HEK293T cells

Next, we extended the assay to assess putatively causal eQTL variants within transcripts using GTEx fibroblast eQTLs, having established that transcript annotations are enriched in these eQTLs (Additional file 4: Figure S3). We chose fibroblasts because GTEx fibroblast transcriptome expression is highly correlated with that of HEK293 cells (rho = 0.68, p < 2.2 × 10−16, Additional file 4: Figure S4). The 33 eQTL variants chosen for editing are located within the transcript of the eGene with which they are associated and have a high posterior probability of causality based on CAVIAR fine mapping. After editing and QC filtering, 13 eQTL variants remained (Additional file 3: Table S3). The variance of the effect size of the eQTL variants was significantly higher than that of the control variants (0.49 versus 0.038; F test p = 1.7 × 10−5; Fig. 2a), which suggests that the edited eQTL variants as a whole have a greater regulatory effect than the edited control variants. Ten of the 13 variants have an effect in the same direction as the GTEx eQTL effect. Five of the eQTL variants are individually significantly different from the control distribution (Fig. 2a and Additional file 3: Table S3), and all five of these variants have an effect in the same direction as in GTEx. Additionally, there is a significant correlation between the effect size of the edited stop-gained, non-eQTL control, and eQTL variants and their effects in GTEx (Spearman’s rho = 0.60; p = 1.9 × 10−4, Fig. 2b), again indicating that the assay captures regulatory effects seen in the population.

In order to establish the robustness of the assay in a different cell line and analyze the potential cell type-specific effects, we replicated the assay for a total of 45 variants in HeLa cells, with 6 eQTL variants, 5 GTEx stop-gained variants, and 5 synonymous control variants passing QC (Additional file 3: Table S3). Due to the small number of control variants in this replication set, we analyzed variant effect sizes rather than distinguishing significant variants in HeLa cells. The effect sizes measured in HEK293T cells and HeLa cells were largely consistent (rho = 0.3, p = 0.258), with all stop variants showing the same direction of effect (Fig. 2c) and small effects of the control variants in HeLa cells. As in HEK293T cells, the HeLa effect sizes showed a correlation with GTEx effect sizes (Spearman’s rho = 0.42, p = 0.074; Additional file 4: Figure S5d), indicating that the assay is robust.

We expect that the lack of effect for some of the eQTL variants is due to the variants not actually being the causal regulatory variants of their association signals. Additionally, our cell line may not perfectly recapitulate the genetic regulatory effects of GTEx fibroblast samples. To investigate this, we looked at the variation in the effect size between GTEx tissues for each of the eQTL variants (Fig. 2d), as well as the variation of the effect sizes between HEK293 and HeLa cells (Fig. 2c, Additional file 4: Figure S5d). We also looked at inter-individual variation within fibroblast samples in GTEx, which may reflect more subtle cell type-specific genetic effects as well as the effects of other regulatory variants that the individuals may have. We measured the effect size in eQTL heterozygotes based on the allelic imbalance within the gene body (Additional file 4: Figure S5c), with eleven of the eQTL variants having sufficient data for this analysis. For all five significant variants in HEK293T cells, there is an agreement in direction between the polyclonal HEK293T aFC, median heterozygous aFC, and eQTL aFC. Some of the variations between the effect sizes observed in HEK293T and HeLa cells may be attributable not only to noise but also cell type-specific regulatory effects. Several of the other variants demonstrate a large range of effects in GTEx both across tissues and across individuals (Fig. 2d and Additional file 4: Figure S5c). The observed effect in the cell line, like an individual or tissue, is likely to fall somewhere in a spectrum of possible effects.

In order to apply our assay to the detection of nonsense-mediated decay triggered by disease-associated variants, we introduced stop-gained variants into two disease-associated genes: ROR2 and GLI3, primarily in HEK293T cells. Seven of the edited stop-gained variants are before the 55-bp threshold and were therefore expected to trigger NMD. Of these variants, all seven resulted in negative effect sizes, and the distribution of these variants was significantly different from that of both the four variants which were not expected to trigger NMD (Wilcoxon p = 6.1 × 10−3) and the non-eQTL control variants (Wilcoxon p = 5.8 × 10−4, Fig. 3a). When tested individually, six of the seven expected NMD variants are significantly different from the control distribution, indicating that we can sensitively detect NMD and NMD escape across the 55-bp boundary in these two genes (Fig. 3a and Additional file 3: Table S3).

Fig. 3
figure 3

The polyclonal assay effectively detects NMD in disease-associated genes. a Effect size in control variants, stop-gained variants after the NMD threshold, and stop-gained variants before the NMD threshold. Triangular points mark variants whose effect size significantly deviates from the control distribution. b The last two exons of NMD disease genes ROR2 and GLI3, showing the effect size (y-axis) and position in the transcript (x-axis) for each successfully edited variant. Disease-associated variants from ClinVar are labeled in red

In addition to the newly created stop-gained variants, we also included disease-causing stop-gained variants from ClinVar. The Arg442Ter mutation in ROR2 results in a stop codon right before the predicted NMD cutoff and is associated with the recessively inherited Robinow syndrome. We observe a significant negative effect of this variant (aFC = − 1.39, Bonferroni-corrected z test p = 2.2 × 10−11), which is consistent with NMD and the clinical manifestation of the disease (Fig. 3b). In contrast, the variant Trp749Ter is associated with dominant type B brachydactyly and falls after the NMD cutoff in the transcript. Our assay shows that Trp749Ter does not affect the expression level of ROR2 and therefore does not appear to be triggering NMD (aFC = − 0.17, corrected p = 1). The one disease variant tested in GLI3, Arg792Ter, falling immediately before the predicted border of NMD escape, shows evidence of triggering NMD with a negative effect size in the assay (aFC = − 1.02, corrected p = 3.6 × 10−6). This result is consistent with the clinical features of this variant, with Grieg cephalopolysyndactyly syndrome thought to be caused by haploinsufficiency in the gene GLI3. The results of editing stop-gained variants in these disease genes indicate that there is a sharp cutoff of NMD/NMD escape at the previously described 50–55-bp threshold and pinpoint the immediate molecular mechanism of NMD/NMD escape for these disease variants. For the three variants that were also tested in HeLa cells and passed QC, the NMD signal was consistent with that in HEK293T cells (Additional file 4: Figure S5d). Additionally, the results demonstrate the potential for utilizing this assay to assess whether a variant of clinical interest triggers NMD when it falls close to the threshold of NMD escape.

Discussion

In this study, we described a method utilizing CRISPR/Cas9 genome editing and targeted sequencing to validate regulatory variants without the need for isolating monoclonal cell lines. We demonstrated our ability to reliably detect the effects of stop-gained variants in the general population and in disease cases with the assay. The ability to experimentally assess the effect of potentially disease-causing stop-gained variants could lead to not only better understanding of the rules of NMD/NMD escape, but also more accurate diagnosis and prognosis. The American College of Medical Genomics recommends caution interpreting the pathogenicity of stop-gained or frameshift variants of unknown significance, especially in cases where the variant is in an exon which might be alternatively spliced, or close to the 3′ end of the transcript [47]. Even though RNA analysis from patients is increasingly used to support variant interpretation [27, 48, 49], establishing causality has been difficult since the lower expression of a mutant haplotype or gene could be driven by other genetic or environmental factors. Our approach can provide evidence that the introduction of the specific variant in question underlies transcript-level changes, thus reducing the ambiguity of variant effects. Furthermore, for genes where NMD/NMD escape is clinically relevant, saturation editing at the 50–55-bp border could build a high-resolution reference for variant interpretation, as it is currently largely unknown how sharp and variable this border is across genes.

This polyclonal assay has the ideal throughput for identifying causal variants from a list of a few to several dozen candidate variants discovered from a rare genetic study. It would be feasible to perform the polyclonal assay on a number of potential regulatory variants, sequencing mRNA and gDNA from the polyclonal culture, and then sort monoclonal cell lines from the same polyclonal culture for only the variants which demonstrate allele-specific regulatory activity. In this approach, the polyclonal assay narrows down the pool of variants to a reasonable number for in-depth follow-up with functional assays, protein quantification, etc. The straightforward nature of the assay makes it easily adoptable in any lab with tissue culture facilities and access to a sequencing instrument.

When we applied the polyclonal assay to eQTL variants, we detected increased variant effects on the expression levels as compared to controls, often in the same direction as the GTEx eQTL effect. Five of 13 variants had significant effects in HEK293T cells, all consistent with the GTEx eQTL data. This clearly demonstrates the ability of our assay to capture common regulatory variant effects. Some of the non-significant eQTL variants appear to have edited effect sizes consistent with GTEx in HEK293T and/or HeLa cells, but we lack the sensitivity to detect these small effects with confidence. In addition, some of the inconsistencies between the assay results and eQTL data are likely to originate from the eQTL data. Since we do not expect fine-mapping to always succeed in identifying the true causal variants at these loci, the undetected effects could represent these situations. Furthermore, with multiple eQTLs for the same gene being common [4], it is possible that eQTL effect sizes observed in populations reflect multiple regulatory variants in partial LD. Therefore, editing a single variant may not yield the same results as the full haplotype. When we looked at the aFC in heterozygous individuals for these variants in GTEx, we found a broad range of effect sizes, suggesting the presence of effects from multiple variants and potential modifiers that may not be captured by editing a single variant. Finally, assessing genetic regulatory effects even in closely matched cell lines does not necessarily capture effects measured in tissue samples. While this is likely to contribute to some of the differences, cis-eQTLs, especially in the transcribed region, are often highly robust across different tissues [4] and are expected to replicate in cell lines as well. We highlight that our approach maintains the genomic context of variants and native gene regulation. Thus, it does not suffer from the limitations of massively parallel approaches where discrepancies between eQTL and experimental data may be due to measuring genetic regulatory effects in artificial constructs [10, 11]. Altogether, more experimentation and further comparison of population and experimental results are required to fully understand the differences between experimental and population data.

Finally, we note that our assay is somewhat limited by HDR efficiency, which varies greatly between loci and led to filtering out a substantial number of targeted loci from our final analysis. Capturing the specific effect of the edited variant requires discarding any reads in the gDNA or cDNA which contain indels created through non-homologous end joining (NHEJ). Since NHEJ often dominates HDR in efficiency, this can result in low numbers of HDR reads. Research in improving the HDR rates in editing is ongoing [50,51,52], and likely HDR efficiency will be greatly improved in the future. Additionally, future improvements on base editor technology, which avoids the introduction of double-stranded breaks and therefore minimizes the risk of indels [53, 54], could also benefit this system and increase the sensitivity of the assay. Finally, while the assay can be applied to multiple cell lines to add robustness of the results and indicate cell type-specific variant effects, cell lines may not fully capture variant effects in vivo.

Conclusions

In summary, we have presented a method to validate the allele-specific effects of regulatory variants using CRISPR/Cas9 in human cell lines. When applied to eQTL variants, we see an increased regulatory effect over the control variants, suggesting we are capturing the allele-specific regulatory effects of these variants. Additionally, all of the significant eQTL variants have effects in the same direction as observed in GTEx, demonstrating the assay’s reliability in detecting eQTL effects. The assay is particularly robust in capturing variants triggering NMD across the genome and in rare disease genes, with potential applications in testing the effects of variants of unknown significance from rare variant studies.

Availability of data and materials

The sequencing dataset supporting the conclusions of the article is available in the figshare repository at https://doi.org/10.6084/m9.figshare.9883232 [55]. The laboratory protocol is available at https://doi.org/10.17504/protocols.io.7c6hize.

References

  1. Lappalainen T, Sammeth M, Friedländer MR, PAC ‘t H, Monlong J, Rivas MA, et al. Transcriptome and genome sequencing uncovers functional variation in humans. Nature. 2013;501(7468):506–11.

    CAS  PubMed  PubMed Central  Google Scholar 

  2. Grundberg E, Small KS, Hedman ÅK, Nica AC, Buil A, Keildson S, et al. Mapping cis- and trans-regulatory effects across multiple tissues in twins. Nat Genet. 2012;44(10):1084–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  3. Battle A, Mostafavi S, Zhu X, Potash JB, Weissman MM, McCormick C, et al. Characterizing the genetic basis of transcriptome diversity through RNA-sequencing of 922 individuals. Genome Res. 2014;24(1):14–24.

    CAS  PubMed  PubMed Central  Google Scholar 

  4. GTEx Consortium, Laboratory, Data Analysis &Coordinating Center (LDACC)—Analysis Working Group, Statistical Methods groups—Analysis Working Group, Enhancing GTEx (eGTEx) groups, NIH Common Fund, NIH/NCI, et al. Genetic effects on gene expression across human tissues. Nature. 2017;550(7675):204–13.

    PubMed Central  Google Scholar 

  5. Nicolae DL, Gamazon E, Zhang W, Duan S, Dolan ME, Cox NJ. Trait-associated SNPs are more likely to be eQTLs: annotation to enhance discovery from GWAS. Gibson G, editor. PLoS Genet. 2010;6(4):e1000888–10.

    PubMed  PubMed Central  Google Scholar 

  6. Hormozdiari F, Kostem E, Kang EY, Pasaniuc B, Eskin E. Identifying causal variants at loci with multiple signals of association. Genetics. 2014;198(2):497–508.

    CAS  PubMed  PubMed Central  Google Scholar 

  7. Wen X, Luca F, Pique-Regi R. Cross-population joint analysis of eQTLs: fine mapping and functional annotation. Gibson G, editor. PLoS Genet. 2015;11(4):e1005176.

    PubMed  PubMed Central  Google Scholar 

  8. Wellcome Trust Case Control Consortium, Maller JB, McVean G, Byrnes J, Vukcevic D, Palin K, et al. Bayesian refinement of association signals for 14 loci in 3 common diseases. Nat Genet. 2012;44(12):1294–301.

  9. Gilissen C, Hoischen A, Brunner HG, Veltman JA. Disease gene identification strategies for exome sequencing. Eur J Hum Genet. 2012;20(5):490–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  10. Tewhey R, Kotliar D, Park DS, Liu B, Winnicki S, Reilly SK, et al. Direct identification of hundreds of expression-modulating variants using a multiplexed reporter assay. Cell. 2016;165(6):1519–29.

    CAS  PubMed  PubMed Central  Google Scholar 

  11. van Arensbergen J, Pagie L, FitzPatrick VD, de Haas M, Baltissen MP, Comoglio F, et al. High-throughput identification of human SNPs affecting regulatory element activity. Nat Genet. 2019;526:68.

    Google Scholar 

  12. Ke S, Shang S, Kalachikov SM, Morozova I, Yu L, Russo JJ, et al. Quantitative evaluation of all hexamers as exonic splicing elements. Genome Res. 2011;21(8):1360–74.

    CAS  PubMed  PubMed Central  Google Scholar 

  13. Wang Y, Ma M, Xiao X, Wang Z. Intronic splicing enhancers, cognate splicing factors and context-dependent regulation rules. Nat Struct Mol Biol. 2012;19(10):1044–52.

    CAS  PubMed  PubMed Central  Google Scholar 

  14. Holbrook JA, Neu-Yilik G, Hentze MW, Kulozik AE. Nonsense-mediated decay approaches the clinic. Nat Genet. 2004;36(8):801–8.

    CAS  PubMed  Google Scholar 

  15. Nagy E, Maquat LE. A rule for termination-codon position within intron-containing genes: when nonsense affects RNA abundance. Trends Biochem Sci. 1998;23(6):198–9.

    CAS  PubMed  Google Scholar 

  16. Rivas MA, Pirinen M, Conrad DF, Lek M, Tsang EK, Karczewski KJ, et al. Human genomics. Effect of predicted protein-truncating genetic variants on the human transcriptome. Science. 2015;348(6235):666–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  17. Miller JN, Pearce DA. Nonsense-mediated decay in genetic disease: friend or foe? Mutat Res Rev. 2014;762:52–64.

    CAS  Google Scholar 

  18. Jinek M, Chylinski K, Fonfara I, Hauer M, Doudna JA, Charpentier E. A programmable dual-RNA-guided DNA endonuclease in adaptive bacterial immunity. Science. 2012;337(6096):816–21.

    CAS  PubMed  PubMed Central  Google Scholar 

  19. Cong L, Ran FA, Cox D, Lin S, Barretto R, Habib N, et al. Multiplex genome engineering using CRISPR/Cas systems. Science. 2013;339(6121):819–23.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. Shalem O, Sanjana NE, Hartenian E, Shi X, Scott DA, Mikkelson T, et al. Genome-scale CRISPR-Cas9 knockout screening in human cells. Science. 2014;343(6166):84–7.

    CAS  PubMed  Google Scholar 

  21. Kosicki M, Tomberg K, Bradley A. Repair of double-strand breaks induced by CRISPR-Cas9 leads to large deletions and complex rearrangements. Nat Biotechnol. 2018;36(8):765–71.

    CAS  PubMed  PubMed Central  Google Scholar 

  22. Li X, Kim Y, Tsang EK, Davis JR, Damani FN, Chiang C, et al. The impact of rare variation on gene expression across tissues. Nature. 2017;550(7675):239–43.

    PubMed  PubMed Central  Google Scholar 

  23. Findlay GM, Boyle EA, Hause RJ, Klein JC, Shendure J. Saturation editing of genomic regions by multiplex homology-directed repair. Nature. 2014;513(7516):120–3.

    CAS  PubMed  PubMed Central  Google Scholar 

  24. Furniss D, Critchley P, Giele H, Wilkie AOM. Nonsense-mediated decay and the molecular pathogenesis of mutations in SALL1 and GLI3. Am J Med Genet A. 2007;143A(24):3150–60.

    CAS  PubMed  Google Scholar 

  25. Johnston JJ, Olivos-Glander I, Killoran C, Elson E, Turner JT, Peters KF, et al. Molecular and clinical analyses of Greig cephalopolysyndactyly and Pallister-Hall syndromes: robust phenotype prediction from the type and position of GLI3 mutations. Am J Hum Genet. 2005;76(4):609–22.

    CAS  PubMed  PubMed Central  Google Scholar 

  26. Schwabe GC, Tinschert S, Buschow C, Meinecke P, Wolff G, Gillessen-Kaesbach G, et al. Distinct mutations in the receptor tyrosine kinase gene ROR2 cause brachydactyly type B. Am J Hum Genet. 2000;67(4):822–31.

    CAS  PubMed  PubMed Central  Google Scholar 

  27. Ben-Shachar S, Khajavi M, Withers MA, Shaw CA, van Bokhoven H, Brunner HG, et al. Dominant versus recessive traits conveyed by allelic mutations - to what extent is nonsense-mediated decay involved? Clin Genet. 2009;75(4):394–400.

    CAS  PubMed  Google Scholar 

  28. Harrow J, Frankish A, Gonzalez JM, Tapanari E, Diekhans M, Kokocinski F, et al. GENCODE: the reference human genome annotation for the ENCODE Project. Genome Res. 2012;22(9):1760–74.

    CAS  PubMed  PubMed Central  Google Scholar 

  29. Zerbino DR, Wilder SP, Johnson N, Juettemann T, Flicek PR. The ensembl regulatory build. Genome Biol. 2015;16(1):56.

    PubMed  PubMed Central  Google Scholar 

  30. Kircher M, Witten DM, Jain P, O’Roak BJ, Cooper GM, Shendure J. A general framework for estimating the relative pathogenicity of human genetic variants. Nat Genet. 2014;46(3):310–5.

    CAS  PubMed  PubMed Central  Google Scholar 

  31. Xiong HY, Alipanahi B, Lee LJ, Bretschneider H, Merico D, Yuen RKC, et al. The human splicing code reveals new insights into the genetic determinants of disease. Science. 2015;347(6218):1254806–6.

  32. Vergoulis T, Vlachos IS, Alexiou P, Georgakilas G, Maragkakis M, Reczko M, et al. TarBase 6.0: capturing the exponential growth of miRNA targets with experimental support. Nucleic Acids Res. 2012;40(Database issue):D222–9.

    CAS  PubMed  Google Scholar 

  33. Oikonomou P, Goodarzi H, Tavazoie S. Systematic identification of regulatory elements in conserved 3′ UTRs of human transcripts. Cell Rep. 2014;7(1):281–92.

    CAS  PubMed  PubMed Central  Google Scholar 

  34. Yang Y-CT, Di C, Hu B, Zhou M, Liu Y, Song N, et al. CLIPdb: a CLIP-seq database for protein-RNA interactions. BMC Genomics. 2015;16(1):51.

    CAS  PubMed  PubMed Central  Google Scholar 

  35. Pickrell JK. Joint analysis of functional genomic data and genome-wide association studies of 18 human traits. Am J Hum Genet. 2014;94(4):559–73.

    CAS  PubMed  PubMed Central  Google Scholar 

  36. Sultan M, Amstislavskiy V, Risch T, Schuette M, Dökel S, Ralser M, et al. Influence of RNA extraction methods and library selection schemes on RNA-seq data. BMC Genomics. 2014;15(1):675–13.

    PubMed  PubMed Central  Google Scholar 

  37. Landrum MJ, Lee JM, Benson M, Brown GR, Chao C, Chitipiralla S, et al. ClinVar: improving access to variant interpretations and supporting evidence. Nucleic Acids Res. 2018;46(D1):D1062–7.

    CAS  PubMed  Google Scholar 

  38. Wang T, Wei JJ, Sabatini DM, Lander ES. Genetic screens in human cells using the CRISPR-Cas9 system. Science. 2014;343(6166):80–4.

    CAS  PubMed  Google Scholar 

  39. Heigwer F, Kerr G, Boutros M. E-CRISP: fast CRISPR target site identification. Nature. 2014;11(2):122–3.

    CAS  Google Scholar 

  40. Arbab M, Srinivasan S, Hashimoto T, Geijsen N, Sherwood RI. Cloning-free CRISPR. Stem Cell Rep. 2015;5(5):908–17.

    CAS  Google Scholar 

  41. Legut M, Daniloski Z, Xue X, McKenzie D, Guo X, Wessels H-H, et al. High-throughput screens of PAM-flexible Cas9 variants for gene knockout and transcriptional modulation. Cell Rep. 2020;30(9):2859–2868.e5.

    CAS  PubMed  Google Scholar 

  42. Yahi A, Hoffman P, Brandt M, Mohammadi P, Tatonetti NP, Lappalainen T. EdiTyper: a high-throughput tool for analysis of targeted sequencing data from genome editing experiments. bioRxiv. 2020. https://doi.org/10.1101/2020.07.30.229088.

  43. Mohammadi P, Castel SE, Brown AA, Lappalainen T. Quantifying the regulatory effect size of cis-acting genetic variation using allelic fold change. Genome Res. 2017;27(11):1872–84.

    CAS  PubMed  PubMed Central  Google Scholar 

  44. Aguet F, Barbeira AN, Bonazzola R, Brown A, Castel SE, Jo B, et al. The GTEx Consortium atlas of genetic regulatory effects across human tissues. bioRxiv. 2019. https://doi.org/10.1101/787903.

  45. Miyaoka Y, Berman JR, Cooper SB, Mayerl SJ, Chan AH, Zhang B, et al. Systematic quantification of HDR and NHEJ reveals effects of locus, nuclease, and cell type on genome-editing. Sci Rep. 2016;6(1):23549–12.

    CAS  PubMed  PubMed Central  Google Scholar 

  46. Doench JG, Fusi N, Sullender M, Hegde M, Vaimberg EW, Donovan KF, et al. Optimized sgRNA design to maximize activity and minimize off-target effects of CRISPR-Cas9. Nat Biotechnol. 2016;34(2):18410.1101/78790391.

  47. Richards S, Aziz N, Bale S, Bick D, Das S, Gastier-Foster J, et al. Standards and guidelines for the interpretation of sequence variants: a joint consensus recommendation of the American College of Medical Genetics and Genomics and the Association for Molecular Pathology. Genet Med. 2015;17:405–23.

  48. Cummings BB, Marshall JL, Tukiainen T, Lek M, Donkervoort S, Foley AR, et al. Improving genetic diagnosis in Mendelian disease with transcriptome sequencing. Sci Transl Med. 2017;9(386):eaal5209.

    PubMed  PubMed Central  Google Scholar 

  49. Kremer LS, Bader DM, Mertes C, Kopajtich R, Pichler G, Iuso A, et al. Genetic diagnosis of Mendelian disorders via RNA sequencing. Nat Commun. 2017;8(1):15824.

    CAS  PubMed  PubMed Central  Google Scholar 

  50. Chu VT, Weber T, Wefers B, Wurst W, Sander S, Rajewsky K, et al. Increasing the efficiency of homology-directed repair for CRISPR-Cas9-induced precise gene editing in mammalian cells. Nat Biotechnol. 2015;33(5):543–8.

    CAS  PubMed  Google Scholar 

  51. Maruyama T, Dougan SK, Truttmann MC, Bilate AM, Ingram JR, Ploegh HL. Increasing the efficiency of precise genome editing with CRISPR-Cas9 by inhibition of nonhomologous end joining. Nat Biotechnol. 2015;33(5):538–42.

    CAS  PubMed  PubMed Central  Google Scholar 

  52. Aird EJ, Lovendahl KN, St Martin A, Harris RS, Gordon WR. Increasing Cas9-mediated homology-directed repair efficiency through covalent tethering of DNA repair template. Commun Biol. 2018;1(1):54.

    PubMed  PubMed Central  Google Scholar 

  53. Komor AC, Kim YB, Packer MS, Zuris JA, Liu DR. Programmable editing of a target base in genomic DNA without double-stranded DNA cleavage. Nature. 2016;533(7603):420–4.

    CAS  PubMed  PubMed Central  Google Scholar 

  54. Gaudelli NM, Komor AC, Rees HA, Packer MS, Badran AH, Bryson DI, et al. Programmable base editing of A•T to G•C in genomic DNA without DNA cleavage. Nature. 2017;551(7681):464–71.

    CAS  PubMed  PubMed Central  Google Scholar 

  55. Brandt, M, Ziosi, M, Gokden, A, Lappalainen, T. polyclonal_editing_fastq_files. figshare. 2020. https://doi.org/10.6084/m9.figshare.9883232.

Download references

Acknowledgements

We would like to thank Ana Vasileva for the advice and the GTEx Consortium, especially Farhad Hormozdiari and Francois Aguet, for the data.

Funding

This work was supported by R01MH106842 and R01GM122924.

Author information

Authors and Affiliations

Authors

Contributions

TL and MB conceived the idea and wrote the manuscript. MB filtered the variants, designed the constructs, and analyzed the results. MB, MZ, and AG carried out the editing and sequencing experiments. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Tuuli Lappalainen.

Ethics declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

T.L. is a scientific advisory board member of Variant Bio and Goldfinch Bio and owns stock in Variant Bio. M.B. is now an employee of Nature Communications, a journal from the same publisher; she did not have any access to or involvement with the editorial process at Genome Medicine. The remaining authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1.

Table S1. Information on the variants chosen for editing.

Additional file 2.

Table S2. Reference and alternate allele counts in cDNA and gDNA sequencing samples after editing variants.

Additional file 3.

Table S3. Effect sizes of edited variants which passed QC filtering steps.

Additional file 4.

Supplementary figure file containing Figure S1-S5.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Brandt, M., Gokden, A., Ziosi, M. et al. A polyclonal allelic expression assay for detecting regulatory effects of transcript variants. Genome Med 12, 79 (2020). https://doi.org/10.1186/s13073-020-00777-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13073-020-00777-8

Keywords