Effects of Pathogenic Variations in the Human Rhodopsin Gene (hRHO) on the Predicted Accessibility for a Lead Candidate Ribozyme

Abstract

Purpose: The mutation-independent strategy for hammerhead ribozyme (hhRz) or RNA interference (RNAi)-based gene therapeutics to treat autosomal dominant diseases is predicated on the hypothesis that a single therapeutic would equivalently suppress all/most of the diverse mutant mRNAs in patients with the disease phenotype. However, the hypothesis has not been formally tested. We address this through a comprehensive bioinformatics study of how mutations affect target mRNA structure accessibility for a single lead hhRz therapeutic (725GUC↓), designed against human rod rhodopsin mRNA (hRHO), for patients with hRHO mutations that cause autosomal dominant retinitis pigmentosa.

Methods: A total of 199 in silico coding region mutations (missense, nonsense, insert, deletion, indel) were made in hRHO mRNA based on Human Gene Mutation Database and Database of Single Nucleotide Polymorphisms. Each mRNA was folded with MFold, SFold, and OligoWalk algorithms and subjected to a bioinformatics model called multiparameter prediction of RNA accessibility. Predicted accessibility of each mutant over both a broad local region and the explicit lead ribozyme annealing site were compared quantitatively to wild-type hRHO mRNA.

Results: Accessibility of the 725GUC↓ site is sensitive to some mutations. For single nucleotide missense mutations, proximity of the mutation to the hhRz annealing site increases the impact on predicted accessibility, but some distant mutations also influence accessibility.

Conclusions: A mutation-independent strategy appears viable in this specific context but certain mutations could significantly influence ribozyme or RNAi efficacy through impact on accessibility at the target annealing site/region. This possibility must be considered in applications of this gene therapy strategy.

Among the many known disease genes (current number: 22; RetNet, https://sph.uth.edu/retnet/sum-dis.htm) that cause autosomal dominant retinitis pigmentosa (adRP), the first gene identified and source of much ongoing investigation is disease-promoting mutations in the human rod rhodopsin mRNA (hRHO) gene (Online Mendelian Inheritance of Man: 180380).^1,2 Rhodopsin, the visual pigment of the rod photoreceptor, constitutes up to 90% of the total protein in the rod outer segment where phototransduction is initiated. A functional RHO protein depends on correct folding, post-translational modifications, cell trafficking, and interaction with 11-cis-retinal to form a ground state visual pigment that can be activated by light. Mutant RHO proteins frequently do not properly fold into functional conformations, blocking mature visual pigment formation.^3,4 This can activate cell stress machinery to promote apoptotic photoreceptor cell death.⁵ This results in the severe, progressive vision impairment that characterizes adRP. The rationale for gene therapy to treat adRP and other autosomal dominant retinal degenerative diseases is to develop therapeutic agents capable of reducing the levels of toxic mutant proteins in the cell, allowing cells to maintain viability and (ideally) function to support vision longer into the life of the affected patient. As gene therapy studies to treat dominant retinal degenerations, and in particular hRHO mutations, are moving toward clinical trial, this study is both relevant and timely. 

The retina has often been selected as a focus for the development of post-transcriptional gene silencing (PTGS) agents. Among these agents are ribozymes, small catalytic RNAs that on binding a target mRNA, undergo a conformational change, and cleave the target at NUH sites (N = any nucleotide, U = U, H = C, U, A). This leads to increased degradation of targeted mRNA and decreased translation (expression) of the protein encoded by it. Although ribozymes are catalytic as RNA elements, short hairpin RNAs (shRNAs) anneal to target mRNAs and recruit host cell proteins to cleave the target mRNA. As gene therapy agents, hammerhead ribozyme (hhRzs) and shRNAs act in trans and interact with a target mRNA by a second-order biochemical reaction based on concentration-dependent biophysical collision frequency. The hhRz or shRNA must colocalize with target mRNA and generally be in concentrative excess to anneal and then cleave the target. 

Challenges arise while developing PTGS agents against target mRNAs, reflected in the very small number of agents that successfully passed clinical trials and are in use.⁶ A major barrier and bottleneck to the development of PTGS agents is the dense structure of any target mRNA, because the first limiting step in PTGS is the ability of any agent to anneal to the target mRNA, which is severely limited by structural inaccessibility.⁷ Most potential hhRz or shRNA cleavage sites in a target mRNA are inaccessible due to dense and stable secondary folding, tertiary folding, protein coating, and dynamical fluctuations.^5–10 Successful annealing requires the target region to be largely single stranded under physiologic conditions. Such sites in target mRNAs are rare. 

Our focus is on hhRzs as therapeutics because they are likely to have greater specificity than shRNAs, which have extensive off-target effects that are potentially toxic in the therapeutic context.^11–14 Because of this, hRHO has been targeted in multiple gene-based therapy development studies attempting to suppress its expression in the context of adRP with the use of ribozymes.^5,10,15–19 

Due to the increasingly vast number of mutations that may be found in a given autosomal dominant disease gene, strategies have been developed to target mutant genes or their mRNAs that are mutation independent.^6,15,16 Rather than attempt to make a new therapeutic for each individual mutation, which is both fiscally and biologically impractical (i.e., most single nucleotide [nt] mutations are buried in the mRNA structure and inaccessible), a single therapeutic is sought that can be used to target and suppress all or most known mutant mRNAs/proteins in a given autosomal dominant disease gene. For instance, if one identifies a highly accessible site in a target mRNA that can be attacked by a ribozyme or other PTGS agents (e.g., Rz725GUC↓), then that therapeutic might be useful to treat most/all known mutations as long as the mutations do not profoundly affect the accessibility of the target site to therapeutic RNA annealing. Because such a therapeutic will also suppress the wild-type (WT) mRNA, the knockdown agent must be combined with an allelic variant expression construct that encodes the WT protein through an mRNA that cannot be cleaved by the given therapeutic (e.g., silent mutation for valine of in-frame cleavable 725GUC↓ to noncleavable 725GUG). We call this strategy “combined knockdown and reconstitute” (CKDRT).⁶ 

hRHO mRNA has multiple pathogenic and tabulated variants thought to be pathogenic (The Human Gene Mutation Database [HGMD]; http://www.hgmd.cf.ac.uk).^20–22 No study to date has used the rich repertoire of known hRHO mutations, or in fact diverse mutations in any gene, to address the question of the viability of a CKDRT approach to treat adRP as a model autosomal dominant retinal degeneration. In this study, we use an approach to map mRNA target accessibility (developed in-house), which we call multiparameter prediction of mRNA accessibility (mppRNA). We used mppRNA to successfully predict potential ribozyme cleavage sites in hRHO in silico, which has identified lead hhRzs candidates with intracellular efficacy.^5,10 MFold is used to predict the most stable secondary structure of a target mRNA through a free energy minimization approach. The single structure with the minimum folding energy (MFE) and an ensemble of structures within a designated degree of higher free energy (less stable) in the folding neighborhood are identified.²³ SFold uses a Boltzmann-weighted sampling of all folding space to estimate the true probability of single stranded-ness at each nucleotide position.²⁴ Another algorithm, OligoWalk (OW), available in RNAStructure, uses the output of MFold and calculates the ribozyme annealing and local structure breakage energies with the concept that a more accessible site is likely to exist as a single-stranded, unstructured form with less energy.²⁵ We take the vectorial outputs of MFold, SFold, and OW and combine them in a simple multiplicative model (set of intersection) to achieve an accessibility map along the target mRNA that has proven versatility.^5,10 

Our lead hhRz candidate against hRHO, Rz725GUC↓ (Fig. 1), targets the GUC↓ cleavage motif at nucleotide 725, which resides in a local region of high-weighted accessibility in hRHO. This lead candidate was identified by first using mppRNA to map accessibility followed by a high throughput screen (HTS) to test a multitude of hhRz expression constructs in cellulo.¹⁰ mppRNA has been successfully used with hRHO to screen for available target sites for ribozyme cleavage.^5,10 Our studies illustrate there is a relationship between mppRNA predicted accessibility and the ability for a ribozyme to cleave in vitro and in cellulo. 

Figure 1

 

Lead candidate ribozyme. Hammerhead ribozyme developed by Yau et al. to target hRHO hhRz cleavage site 725GUC↓. The lead candidate is composed of 7-nt upstream and downstream antisense flanking regions and a 4-bp length helix II capped by a GAAA tetraloop. This ribozyme was selected among 32 tested ribozymes for cleavage in cellulo via secreted alkaline phosphatase assay as previously described.¹⁰

The human hRHO gene is subject to hundreds of pathogenic changes, including single nucleotide missense/nonsense; small insertions, deletions, or indels and gross insertions or deletions. Most of these mutations result in adRP.^2,26 It is established that a mutation within the target binding site of a ribozyme can decrease its efficacy by either altering its rate of binding or release of products from the antisense flanking regions or directly altering its enzymatic ability to cleave at the target if the mutation affects the NUH site directly.^27–30 Here, we took advantage of our predictive mppRNA model to explore the possibility that pathogenic missense mutations, nonsense mutations, insertions, or deletions can lead to a sufficient structural change in a target mRNA that could perturb accessibility for our lead PTGS agent (Rz725GUC↓). Because prospective patients who would be receiving a therapy of this sort will have varied mutations in the gene being suppressed, understanding the effect of mutations on the ability of the single agent to knock down levels of mutant mRNA is highly important for the eventual clinical application of PTGS agents for adRP or other autosomal dominant diseases. 

Materials and Methods

Mutants and Ribozyme

The in silico hRHO mRNA transcript used was a 1532-nt mRNA originally cloned and sequenced by Nathans and Hogness (GenBank Entry NM_000539.3).²⁶ The transcript has a 95-nt 5′untranslated region (UTR), full coding region (1044 nt), stop codon (UAA), and a 394-nt 3′UTR that contains the first dominant poly-A signal (nt 1506) and continues 26 nt past the polyA signal to nt 1532 at the approximate site where 3′ cleavage and polyA addition would naturally occur. PolyA tails (∼200 to 250 nt) have no structure and therefore are not included in predictive folding experiments. This 1532-nt sequence is the dominant form of polyadenylated hRHO mRNA made in vivo.²⁶ The dominant mature hRHO mRNA is approximately 1.7 to 1.8 kb due to unstructured polyadenylation at the 3′ end. Nondominant transcripts exist but were not used for in silico analysis in this study. All hRHO mutations were made in silico in the 1532-nt dominant form of polyadenylated hRHO mRNA. 

The lead ribozyme candidate for an mutation-independent (MI) CDKRT lead gene therapy strategy targets the GUC↓ site at position 725 in the coding region of the hRHO transcript and is composed of 7-nt antisense flanking regions on each side of the cleavage (nonhybridizing) nt (C of GUC↓), and a 4-bp-length Stem II helix capped by a GAAA tetraloop (Fig. 1).¹⁰ This agent was the best-performing hhRz among 32 tested ribozymes for cleavage in cellulo in a HTS in which the hhRzs were embedded in an engineered high-flexibility region of the adenoviral viral-associated RNA I (VAI) scaffold RNA, demonstrated as a 60% mRNA knockdown measured by secreted alkaline phosphatase reporter assay.¹⁰ The NUH↓ sites selected for screening in the prior studies were chosen based on their predicted accessibility (rank ordered) in the hRHO mppRNA map. The 725 region had the highest level of predicted accessibility, which correlated with the level of experimental target suppression. 

A total of 199 hRHO mutations in the coding region of the mRNA (nts 96 to 1142) were surveyed in this study and were found in the HGMD, National Center for Biotechnology Information's (NCBI) Database of Single Nucleotide Polymorphisms (dbSNP), with one novel mutation characterized from an adRP patient seen in Buffalo, NY (Sullivan et al., unpublished observations, 2016).^20,22All known single nucleotide missense and nonsense mutations found on HGMD that were pathogenic for adRP were used in this study (Fig. 2). Notably, we evaluated six missense mutations at the known hotspot amino acid residue P347, which all cause severe early-onset adRP and the P23H mutation responsible for the most common form of hRHO adRP in the United States.^2,31–33 A smaller number (15) of nonpathogenic single nucleotide polymorphisms (SNPs) within 200 nts of the target nucleotide 725 were also sampled from NCBI's dbSNP that did not appear on the HGMD list to evaluate the impact of local or regional variations on PTGS target accessibility. Pathogenic small insertions, small deletions, small indels, and large insertions or deletions were also explored, totaling 24. Insertions, deletions, and indels were generated from those listed at HGMD.^34,35 

Figure 2

 

Distribution of known hRHO mutations. An illustration of the types and locations of 199 mutations surveyed for predicted cleavage accessibility via mppRNA analysis. All known mutations are located in the coding region of the mRNA. The types of mutations are color coded. Sites with multiple surveyed mutations are represented as having an increased frequency. For mutants affecting multiple nucleotides, the first nucleotide affected was depicted. Mutations were selected from HGMD and NCBI's dbSNP. Mutations including single nucleotide point mutations, small and large insertions, deletions, or indels were sampled. The common P23H single nucleotide missense mutation responsible for most cases of adRP in the United States, as well as six single nucleotide missense mutations causing mutations at amino acid P347, a severe form of adRP, were surveyed.

mppRNA Approach to Target mRNA Accessibility Mapping

A mutant hRHO sequence was created for each of the mutations from the 1532-nt WT hRHO mRNA using the information in the databases and the source papers and folded using three secondary structure predictive algorithms. All foldings were done under standard conditions commonly used for these algorithms (37°C, 1 M salt). Briefly, the first algorithm, MFold (http://mfold.rna.albany.edu/), which examines the frequency of being single stranded in all structures in the acquired ensemble, uses energy minimization to predict the MFE structure and an ensemble of structures of higher (less stable) free energies within a set criterion range of energy difference from the MFE structure (parameters used: 10% window of free energy from the MFE structure, difference window size = 3, up to 100 structures in the ensemble).²³ The MFE secondary structural fold generated from MFold for WT hRHO shows the dense secondary structure of a typical target mRNA for PTGS therapeutic development (Fig. 3). The output of MFold used for mppRNA is the normalized ssCount frequency vector (Fig. 4A), which represents the percentage of predicted structures where each nucleotide is single stranded. This is a biased probability estimate of accessibility only in the structural neighborhood of the MFE. The second algorithm, SFold (http://sfold.wadsworth.org/), uses a Boltzmann-weighted sampling of the total structural state ensemble and outputs a direct true estimate of the probability of being single stranded along the target string over the entire folding space.²⁴ The output of SFold used for mppRNA is the sStrand vector (Fig. 4B). Another output of SFold shows how the sample of structures (1000/sample) clusters into discrete distributions with repeatable probability densities. A method of displaying such distributions is through multidimensional scaling (MDS), which in effect collapses the three-dimensional (3D) map of location in the space (base pairing differences) versus free energy into a plane with unitless dimensionless axes that displace the distribution of structures.^36,37 MDS is a way of looking at how single nucleotide polymorphisms can affect RNA structure.^38–42 The third algorithm requires the RNAStructure package (http://rna.urmc.rochester.edu/RNAstructure.html). RNAStructure has a version of MFold that generates an ensemble of structures. The OW algorithm uses the MFold output and calculates the local target unfolding energy (target break energy) for overlapping windows (of arbitrary size) of RNA over the entire length of mRNA.²⁵ The 15-nt window was used in OW to represent the span of annealing for a standard hhRz (including the Rz725GUC↓) with 7-nt antisense flanks on either side of the cleavage nucleotide (C↓; Fig. 1), which does not base pair with the hhRz. The OW output (.rep file) displays the free energy (kcal/mol) of the unfolding energy for each overlapping 15-nt segment (1-nt steps) along the entire length of the mRNA (Fig. 4C); it is already a filtered version because of the stipulated window size for a standard symmetrical hhRz. We conduct a linear transform to convert this energy vector into a unitless positive probabilistic scale by adding a scalar constant to all centered points and then normalizing the output by the same constant. The constant chosen was 41.6 kcal/mol or the absolute value of the total free energy (−41.6 kcal/mol) of the RNA sequence 5′-GGGGGGGGGGGGGGGAAAAACCCCCCCCCCCCCCC-3′ (calculated via RNAStructure). This sequence folds into only one predicted low energy structure, a stem loop of 15 high energy base pairs (G≡C), connected by a small single stranded loop of five A residues. The total energy of this structure was used for normalization of the OW output, as it represents the maximum folding free energy of a 15-nt stable segment of folded mRNA (such a sequence does not occur in hRHO, but this is used because it is a normalization parameter that can extend beyond this particular target in future comparisons). In doing this normalization, we convert the free energy scale into a unitless scale between 0 and 1, which is an additional probability estimator of accessibility (local high free energy regions are likely to be single stranded). 

Figure 3

 

Example of 2D fold of WT hRHO mRNA generated using MFold. (A) Full 1532-nt hRHO mRNA. The MFE structure is shown. (B) Expanse of the 725 region. Blue arrow indicates the 725GUC↓ cleavage site.

Figure 4

 

Independent accessibility maps. (A) MFold, (B) SFold, and (C) OligoWalk outputs. An mppRNA map is generated by multiplying probabilities of each output for each nucleotide position. Large span representing the local regional integration (blue bar) (640 to 764 nt) and hhRz anneal site representing the ribozyme site integration (red bar) (718 to 732 nt) are indicated.

The raw MFold and SFold outputs are then smoothed by a 15-point averaging window to substantially decrease probability of a 0 accessibility score at one or a few nucleotides in one vector output from being projected into the final mppRNA multiplicative product vector output. The 15-point smoothing also allows the MFold and SFold outputs to represent, nucleotide for nucleotide, the accessibility of the same 15-nt segment as the OW output along the entire length of the mRNA. The 15-nt smoothed vector outputs of each of these three algorithms (0 ≤ X ≤ 1) are then multiplied together, with current assumption of equivalent weighting among the output vectors, and an mppRNA map is generated that represents the intersection access probability at each nt residue (Fig. 5). mppRNA seeks regions that are accessible as assessed by all three algorithms (MFold, SFold, and OW) (intersection set of the three probability estimators). Note the broad local region of accessibility (640 to 764 nts) that includes the targeting site for Rz725GUC↓ (718 to 732 nts). 

Figure 5

 

hRHO mppRNA full map. (A) Full-length mppRNA map of WT hRHO, using vector outputs from MFold, SFold, and OligoWalk. (B) Expanse of the local 725 region. The 725GUC↓ cleavage site is located in a region of high predicted accessibility, identified as a mesa or elevated plateau in the mppRNA map. Blue bar (near abscissa) represents the local regional integration from nts 640 to 764. Red bar represents the ribozyme binding site integration from nts 718 to 732.

By integrating the area under the curve of the 15-nt span of the mppRNA map centered at the site of cleavage, one can generate a quantitative scalar value that represents the overall predicted accessibility for the explicit target region, referred to here as the ribozyme site integration (718 to 732 nt). mppRNA was performed on WT hRHO and each mutant, and the site integration was determined for each mutant. For mutants that cause the numbering to shift, the area around the residue that corresponds to nt 725 in WT hRHO was integrated. A larger integration of the entire peak representing the broad area of increased accessibility, referred to here as the local integration (640 to 764 nt), was performed to assess the impact of the mutations on a broader area of the mRNA secondary structure. Because it has been shown that efficacy of ribozymes against hRHO correlates to the areas of peaks in the mppRNA map, it is possible to make judgements about altered efficacy of our lead hhRz among various mutant genotypes.^5,10 

Quantitative Analysis

Statistical and graphical analysis was conducted in Origin (Versions 6.1, 8.1, and 2016; MicroCal, Northampton, MA, USA) using ANOVA and post hoc t-tests for between-condition comparisons with a previously chosen significance level of 0.05. Tests for uniformity of variance and normal distribution of data were conducted (uniformity of variance: Levene's, Shapiro's, Wilk's; Gaussian distribution: Anderson-Darling, Shapiro-Wilk, Kolmogovov-Smirnov). Post hoc parametric Student t-test and Welch's t-test were conducted. Welch's t-test does not depend on equivalence of variance and is robust to non-Gaussian distributions. Degrees of freedom (df) are shown where pertinent. 

Results

Single Nucleotide Missense/Nonsense Mutations

The two integrations in both a global (total mRNA) and a local mppRNA map of the WT hRHO mRNA are displayed in Figure 5. The global mppRNA map (Fig. 5A) shows that accessibility varies considerably throughout hRHO mRNA. The most accessible region (peak with the greatest integration weight) of the mppRNA map is broad and covers the span of nucleotides 640 to 764, which constitutes the local integration region that is illustrated by the blue bar underneath the expanded region of the global map (Fig. 5B). Within the local area is the actual 15-nt target annealing site (718 to 732 nts) in hRHO mRNA for our lead hhRz (Rz725GUC↓; Fig. 5B, red bar). Both the local and annealing site integral values from the mppRNA map were analyzed in the context of mutational effects on accessibility for a single lead hhRz candidate therapeutic (Rz725GUC↓). 

A total of 176 in silico transcripts of single nucleotide point hRHO mutations were analyzed for impact on accessibility within the local and 725 hhRz annealing site spans in the mppRNA maps. The data for the MFold, SFold, OW, and mppRNA maps for the WT and each hRHO mutant are provided (Supplementary Table S1). Qualitatively, the region of accessibility around the 725 cleavage site was maintained in most single nucleotide mutations (Supplementary Table S2). Most single nucleotide mutations exerted little effect on either local or site accessibility relative to WT hRHO mRNA (17 of 176 [9.7%] showed a statistically significant difference in site accessibility by t-test; one showed a statistically significant difference in the local accessibility by t-test [U755G, relative integration 0.84803, t(248) = −2.28418; P = 0.02321]; Supplementary Table S2). The distribution of predicted accessibility for each of the 176 single nucleotide missense/nonsense mutations is shown (Fig. 6). Accessibility within or around the 725 hhRz target site could theoretically increase or decrease due to a proximate or remote mutation in the hRHO mRNA. The local accessibility of the mutants remained within ±15% of WT and a line with nonsignificant slope and y-intercept approximating 1 could be readily fit through the data (Fig. 6A). Sequence regions where mutations show the largest impact on 725 accessibility are within the actual hhRz binding site and those flanking an insulator sequence from nucleotides 500 to 600. Mutations within the insulator region show low impact on accessibility in the 725 region. An area at the tail end of the coding region (nts 1116 to 1129) also has a high number of mutations with some impact on 725 accessibility; however, none of these single nucleotide mutations exerted a statistically significant shift in accessibility in the local integration analysis (Supplementary Table S2). In contrast, the hhRz annealing site accessibility around 725 was substantially impacted by mutational variation, with a range from +28.2% to −67.9% (Fig. 6B; + values indicate enhanced accessibility and − values indicate decreased accessibility). Not surprisingly, the largest areas of change are mutations immediately within and bracketing the annealing sequence of the Rz725GUC↓ cleavage site. There are no known mutations at the 725 cleavage nucleotide reported to date. Several mutations (G720A, G723U, and A727U) within the 725 annealing site region exert a statistically significant difference on the integrated site accessibility values. This effect may be nucleotide-specific because an alternative mutation, A727G, does not exert a significant change in site accessibility. There were significant changes in annealing site accessibility from single nucleotide mutations surrounding the hhRz annealing span (U715A, U715G, C738A, A741U, and U742G) and from single nucleotide mutations bracketing the insulator region (U469G, C488G, C498U, U639A, U654G, G655A, G658A, and C673U). There is a consistent relationship between the proximity of the hRHO mutation and its impact on accessibility at the site of hhRz cleavage (Table 1; Fig. 6). Other hRHO mutants that promoted statistically significant change in accessibility at the 725 annealing site include, surprisingly, the single nucleotide mutation U355A, which is remote from the 725 region, and the large insertion mutation, 1032ins150del8, which may be expected to largely alter overall secondary structure of the hRHO mRNA fold. 

Figure 6

 

Distribution of accessibility relative to WT hRHO. (A) Local integration and (B) site integration. Relative accessibilities were found by normalizing the area under the curve (integral) for each mutant by the area under the curve for the WT mppRNA over the appropriate span of sequence (mutant accessibility/WT accessibility). Local integration shows minimal influence of mutation on structure. A line (red) could be fit through the data with essentially zero slope and a y-intercept of 1 (slope: 2.9776E-6 ± 1.10249E-5; intercept: 0.98886 ± 0.00773). Site integration indicates a substantially greater impact of mutation. An insulator sequence is observed from approximately nts 500 to 600, flanked by regions of high influence on accessibility at the 725 region. Another region of high accessibility is found at the end of the coding region (1114 to 1140; coding region is 96 to 1142).

Table 1

 

Accessibility of Selected Missense Mutations in hRHO mRNA

 

Table 1 contains noteworthy examples of single nucleotide mutants that have different effects on the predicted accessibility of the mRNA, chosen based on the magnitude of changes in accessibility, shape of the local/site mppRNA map, statistical significance, and clinical relevance. Some mutants had considerable differences of predicted accessibility in either one or both regions of integration seemingly dependent on the mutation's proximity to the 725 residue. Maps of the single nucleotide mutations with the largest relative changes from WT hRHO are presented (Fig. 7). Mutation A664G (Asp190Gly) shows the second highest local accessibility (112.4% relative to WT, t[df = 248] = 1.4189; P = 0.15718). The site accessibility also increases but is also not statistically significant (116.9% relative to WT, t[28] = 1.42455; P = 0.16534). A664G is located in the large region of susceptibility just after the insulating region. Likewise, mutation U755G (Phe220Leu) shows the largest significant decrease (84.8% of WT, t[248] = −2.28418; P = 0.02321) in accessibility of the local integration. The site accessibility change shows a nonsignificant increase (114.9% of WT, t[28] = 1.64737; P = 0.11066). The distant mutation, G1116A (Glu341Lys), has the greatest but statistically nonsignificant increase on the local integration (113.3% of WT, t[248] = 1.64317; P = 0.10162). The site integration for G1116A shows a nonsignificant increase in accessibility (114.5% of WT, t[28] = 1.05567; P = 0.30014). The site accessibilities of the mutants were subject to larger changes than the local accessibilities relative to WT. The largest statistically significant variations occurred when the mutation was in close proximity to the site of ribozyme annealing and cleavage. Two mutations illustrate this: G720A (Val209Met), which has the highest site accessibility of all mutants (128.2% of WT, t[28] = 2.39529; P = 0.02354), and G723U (Val210Phe), which has the lowest of all (32.1% of WT, t[28] = −7.89601; P = 1.34E-8). Although these mutations are within the antisense binding region of our ribozyme, we acknowledge the predicted accessibility measured here does not address the potential impact of mutation on the thermodynamics and kinetic behavior of 725 hhRz due to a mismatch of base pairing in helices I and III.⁴³ These examples likely are representative of potential unreported mutations. 

Figure 7

 

Select mutants showing different results. All 199 surveyed mutations had mppRNA maps generated, and area under the curve integrations were performed over the local regional (blue bar) and ribozyme site (red bar) regions. In these select mutants, qualitative and quantitative shifts in the regional/local structure of the mppRNA map are observed in the vicinity of the 725 hhRz targeting site. This is an indication for changes in stability over the target region. Different mutants were compared with the same region of the WT mRNA mppRNA map. (A) WT hRHO. (B) A664G. (C) G720A. (D) G723U. (E) U755G. (F) G1116A.

The 725 target site in hRHO resides approximately at 50% of the length of the mRNA. We explored adRP mutations near the (remote) 5′ and 3′ ends of the coding region. The mutation C163A causes Pro23His, the most common mutation causing adRP in the United States. The predicted accessibilities, both local and site, for this mutation did not differ significantly from that of WT hRHO (Table 1; Fig. 8A) (local accessibility, 102.7% of WT, t[248] = 0.35738; P = 0.72111; site accessibility, 107.2% of WT, t[28] = 0.57014; P = 0.57313). Multiple mutations at codon 347, including Pro347Ser, are responsible for a severe early-onset and rapidly progressive form of adRP. These P347X mutations (P347A, P347L, P347Q, P347R, P347S, and P347T) occur through changes at nts 1134 and 1135 as shown in Table 1. The mutation Pro347Ser (C1134U) showed the largest decrease in relative site accessibility of the six mutations (84.3% of WT, t[28] = −1.80676; P = 0.08156; Fig. 8B). The mutant C1134A responsible for Pro347Thr also showed a nonsignificant decrease (88% of WT, t[28] = −1.32941; P = 0.19444). Other P347X mutants did not show changes in relative site accessibility substantially different from hRHO. The effects of P347X mutations on local accessibility quantification show values within 3% of WT that were not significant. At distances of more than 400 nts, different single nucleotide mutations can promote some changes in the accessible architecture of a targeted region, but none have proven statistically significant. 

Figure 8

 

mppRNA maps of select clinically significant hRHO mutants. (A) Local regional map of the WT hRHO mRNA. (B) Single nucleotide missense mutation C163A, representing the common adRP mutant P23H. An insignificant increase in site accessibility of 7.2% is noted (t[28] = 0.57014; P = 0.57313), and the local map architecture is maintained compared with WT. (C) A map illustrating the mutation C1134U, responsible for the P347S mutation. A decrease of 15.7% is noted for the site integration (t[28] = −1.80676; P = 0.08156). Local regional (blue bar) and ribozyme site (red bar) regions are indicated by the bars.

To further evaluate quantitative accessibility differences between WT and mutant hRHO mRNAs, the integrated vector weights from the site and local regions surrounding the 725GUC↓ cleavage site were compared using a two-way t-test (Supplementary Table S2). Of the 199 mutations, 18 (9.0%) lead to statistically significant changes in the site 725 accessibility (17 missense mutations and one large insertion [1032ins150del8]). For the local integration around 725, there are only two mutations (U755G and1032ins150del8, 1.0%), which promote a statistically significant change in accessibility. 

Because only a relatively small subset of mutations could affect structural accessibility in hRHO targeting by our lead 725 hhRz PTGS agent, this is evidence in support of the MI-CKDRT therapeutic strategy for this construct. 

To explore the extent to which the site accessibility is related to the local accessibility, a linear regression was performed on the accessibility predicted for each mutant by each method and plotted as site versus local integration. A significant positive correlation was observed (R = 0.51531, P < 0.0001), suggesting that the site and local accessibility around nt 725 follow similar trends when affected by a given mutation (Fig. 9). This is consistent with the idea that RNA folds locally and that the annealing site of the 725 hhRz is a subset (nested within) of the larger local regional accessible region. 

Figure 9

 

Comparison of 15- vs. 125-nt integrations. A linear correlation (R = 0.51531, P < 0.0001) is noticed between changes in the local integration versus the ribozyme site integration when each is normalized to the integral for WT hRHO for that region. This suggests that, although related (the ribozyme site sequence is nested within the local regional sequence), there is some independence between changes in the accessibility on a ribozyme binding site versus a regional local scale.

We further investigated the class of mutants that exert statistically significant effects and the class of mutants that do not exert statistically significant effects for their impact on the 725 region (Fig. 10). The mppRNA map of the 725 region of the WT hRHO mRNA is shown for comparison (Fig. 10A). The average mppRNA map in the 725 region of all of the hRHO mutants that exert statistically significant effects demonstrates substantial differences relative to the WT (Fig. 10B). The average mppRNA map in the 725 region of all of the hRHO mutants that do not exert statistical significant effects shows similarity to the WT hRHO mRNA map (Fig. 10C). The average site integrals for both classes were compared with WT and the mean of only the statistically significant class was different from WT (Fig. 10D). Similarly, the average local integrals for the two classes of mutants were compared with WT, but the means were not significantly different (Fig. 10E). It is clear that the individual effects of certain mutants on the annealing site accessibility are represented as an averaged class effect. 

Figure 10

 

Comparison of local structural mRNA impact by mutant classes. We compare the effect on structure of the class of mutants that exert statistically significant impact on the 725 region versus the larger class of mutants that do not exert statistically significant impact, in comparison to the WT mRNA. (A) Local mppRNA map of the WT hRHO. (B) Local average mppRNA map of all hRHO mutants that exerted statistically significant effects on local/regional mppRNA maps around the 725 target site. (C) Local average mppRNA map of all hRHO mutants that did not have statistical effect on 725 region accessibility. (D) Quantitative statistical comparison of accessibility in the immediate hhRz binding site region (nt 718 to 732) for WT versus class of statistically significant mutants (SS) and the class of nonstatistically significant mutants (Non-SS) (mean ± SEM, P values; WT mean = 0.19342 ± 0.01346; SS mean = 0.13653 ± 0.00865; Non-SS mean = 0.19413 ± 0.01404) (SS: t[28] = 3.55571; P = 0.00136, significant; Non-SS: t[28] = −0.03625; P = 0.97134). Asterisk indicates statistical significance. (E) Statistical comparison of accessibility in the local region of hhRz binding (nt 640 to 764) for WT versus mutants with and without significant impact (WT mean = 0.1982 ± 0.0105; SS mean = 0.18192 ± 0.00968; Non-SS mean = 0.19745 ± 0.01026) (SS: t[248] = 1.1402; P = 0.2553; Non-SS: t[248] = 0.05082; P = 0.95951).

To evaluate the potential impact of a remote mutation that had a significant impact of the site accessibility we further investigated the U355A mutation. This mutation occurs in a region of substantial stable secondary structure and within the 356 GUC↓ cleavage site that does not suppress with hhRzs.⁵ We first explored the MFE structure generated by MFold for the U355A mutant versus the WT mRNA (Fig. 11A). An impact of U355A is seen on the local structure relative to the WT mRNA but does not extend into the 725 region of the mutant mRNA, which shows identical local structure to the WT mRNA. Evaluating the ensemble of possible structures the U355A mutation has an effect on the region around the 725 cleavage site and the mppRNA site integral shows significantly different accessibility relative to WT mRNA but the local integration did not (Fig. 11B). Reasoning that the effect of the mutation on remote accessibility could be reflected in the ensemble of all potential structures, we used SFold to investigate for propagated differences in structure due to the mutation (Fig. 11C). Using the multidimensional scaling (MDS) analysis of SFold, one can see that the WT hRHO mRNA demonstrates two clustered distributions of structures and the single MFE structure resides in only one of those clusters. The U355A mutation promotes a change in the probabilistic distribution of structures, which is consistent with a larger scale effect of the mutation. 

Figure 11

 

Impact of a single mutation (U355A, Val87Asp) on global mRNA structure. (A) Local regions of MFE structures of the hRHO WT mRNA and U355A mRNA. The MFE structures had similar free energies and global structures (WT: −530.89 kcal/mol; U355A: −535.48 kcal/mol; Supplementary Fig. S1). The mutation (red arrows) has a local impact on structure in U355A versus WT mRNA in the same region but does not impact the local structure around the 725 hhRz target cleavage site (blue arrows in U355A and WT mRNAs). (B) mppRNA comparison of the WT hRHO mRNA map and the U355A map in the region around the 725 target site. The site integration was significantly different for the U355A mutation relative to the WT hRHO mRNA (t[28] = −2.31025; P = 0.02846), whereas the local/regional integration was not (t[248] = −1.2809; P = 0.20143). (C) Multidimensional scaling plots generated by the SFold RNA statistical module for the WT hRHO and mutant U355A mRNAs. WT hRHO mRNA shows two clusters of sampled structures, and the MFE structure is located in the cluster with the highest probability (WT: two clusters; probabilities: 0.543* [red dots], 0.457 [green dots]; *cluster with MFE structure; second trial: probabilities 0.537, 0.463 [data not shown]). U355A hRHO mRNA also shows two clusters like the WT hRHO mRNA, but the distributions are substantially different (U355A: two clusters; first trial: probabilities: 0.847 [red dots], 0.153* [green dots]; second trial: probabilities: 0.861 and 0.139 [data not shown]). Note that the low probability cluster (green dots) now containing the MFE is substantially diffuse compared with the lower probability cluster for the WT mRNA. Each MDS analysis was repeated twice to demonstrate the repeatability of clusters in the sampling of 1000 structures for both the WT and U355A mRNAs. These outcomes suggest that the U355A single nucleotide mutation has a more global effect on the entire set of possible structures engaged by the hRHO mRNA.

Figure 11

 

Continued.

Other Forms of Mutations

Of the 23 samples that included small insertions, deletions, or indels, only one showed a statistically significant effect on site or local accessibility. This mutation (1032in150del8), a large insertion of 150 bp, reduced the site accessibility by 51.3% (t[28] = −6.25557; P < 0.05).³⁵ 1032ins150del8 is also the only insertion, deletion, or indel to reduce the local integration significantly, by 26.7% (t[248] = −3.99756; P < 0.05). Other than this very large insertion, eight others had a change site in accessibility of 10% or more from WT. Of these, seven are insertions or deletions of 3 or more nts, whereas only one was a single nucleotide insertion. Table 2 lists the relative accessibilities of each insertion, indel, or deletion. These results contradict the thought that mutations involving more than 1 nt would have larger influence on accessibility surveyed in the area proximal to the 725GUC↓ cleavage site (four total in the area of nt 725 ± 200). 

Table 2

 

Relative Accessibilities of Each Insertion, Indel, or Deletion

 

Discussion

PTGS therapy has great potential in the treatment of a wide variety of diseases, including adRP. Part of understanding the complexity of developing these treatments includes understanding how the molecular genetic basis of the disease might affect therapeutic efficacy. This bioinformatics study shows that using a predictive model of target RNA structure and ribozyme target accessibility provides expedited insight on how mutational variations that cause adRP could affect ribozyme annealing and subsequent target mRNA cleavage. Although not tested in this study, shRNA- or miRNA-type therapies also target local regions, and one would expect similar outcomes. shRNAs typically use a span of 19 nt for antisense binding, 4 bp larger than sampled for the standard 7 nt/7 nt hhRz tested here in the context of our lead agent with therapeutic potential. We found that most of the comprehensive hRHO mutations surveyed had no statistical effect on predicted accessibility at or around the Rz725GUC↓ cleavage motif. However, a small subset of mutations was shown to either increase or decrease accessibility around the 725GUC↓ site, and for most mutations the variability is not extreme. Overall, these outcomes are highly favorable for a mutation-independent PTGS strategy. 

As expected, the effect on accessibility at a single target annealing site increases with proximity of the mutation to the region of annealing. However, it was also shown that single nucleotide mutations in relatively distant positions can have an appreciable effect, as shown with the P347 mutation, the two mutations at nt 498 that reduced accessibility, and the global impact of mutation U355A. The effect of single nucleotide mutations on the efficacy of ribozymes has only been explored in the context of a mutation being located within the antisense binding region, where it directly affects initial target annealing and product leaving rates, or directly affects the NUH site, which has general inhibitory effects on the rate of ribozyme cleavage.²⁸ Our study surveyed a total of 199 mutants, 176 of which were single nucleotide variations. The results of this study strengthen the potential utility of the mutation-independent PTGS strategy but also raise caution that some mutations may affect accessibility in a target mRNA and alter the expected impact of cleavage-active lead therapeutic hhRzs (or other PTGS agents). Prior to using our Rz725GUC↓ therapeutic agent on mutations that induce significant predicted accessibility differences, it would be prudent to experimentally evaluate accessibility to annealing and cleavage on mutant hRHO mRNAs both in vitro and in cellulo. Only five of the known hRHO mutations directly affect nucleotides within the 15-nt Rz725GUC↓ annealing region in the hRHO transcript, and only one directly modifies the hhRz NUH↓ triplet recognition motif (G723U). Knowledge that only a small subset of known hRHO mutations could affect accessibility to ribozyme annealing, a factor critical to therapeutic efficacy, shows that a mutation-independent hhRz therapy targeting the 725GUC↓ site has utility for gene therapy in a large pool (>90%) of individuals with hRHO mutations. For those mutations that have a significant predicted effect on accessibility at the Rz725GUC↓ site/region, an experimental approach to evaluate accessibility and efficacy of the given therapy against model mRNAs in vitro and in cellulo is a reasonable approach. 

We also found that the ability for a hRHO mutation to affect target site accessibility is dependent on the domain the mutation occupies within the transcript. This is evidenced by areas of insulation where mutations had little to no effect on accessibility at the 725 site. The dominant area of insulation (nts 500 to 600) is an area of strong negative folding energy. This likely affects the impact of local single nucleotide mutations on the nearby 725 region as they are unable to propagate structural destabilization out of the insulating region. There were also areas in the target (nts 400 to 500 and 600 to 800) where accessibility around 725 was altered by mutation, albeit not always in a statistically significant manner. Although the impact of any mutation is most likely to have local structural impact on at least secondary RNA folding, as we showed in this study, there is also precedence for longer-range impacts of mutations. Several studies have shown that single nucleotide polymorphisms are most likely to have local effect on RNA structure but can also have longer range effects.^38–42 Having comprehensively surveyed the known hRHO mutations, it is reasonable to conclude that our results are representative of the impact of single nucleotide polymorphisms on mRNAs as found in other studies. We expect that newly identified mutations in hRHO will follow the same paradigm identified for the current set of variations as to how they could impact target mRNA accessibility. 

This study is limited to work with a single predictive model (mppRNA) to identify accessible regions in hRHO that are amenable to attack by hhRzs, both in vitro and in cellulo, and focuses on one known accessible region in hRHO at which successful PTGS agents have already been realized. mppRNA is a proven useful method of screening for target accessibility to design, screen, and develop ribozymes for cleavage of hRHO mRNA.^5,10 In its current state, the mppRNA model identified our lead candidate hhRz for hRHO: Rz725GUC↓. Although the model is proven useful in its current state, it is possible that predictive output could be improved by future bioinformatics studies. At this time, we know that mppRNA has predicted accessibility in hRHO that has resulted in meaningful therapeutic mRNA cleavage and knockdown by the lead 725 hhRz.¹⁰ The 725 region is therefore a valid reference site to evaluate the impact of hRHO mutations on lead agent target accessibility. There are other regions in hRHO to which we and other investigators have made successful PTGS agents. A broader bioinformatics study to evaluate the now established database of accessibility in the set of known hRHO mutants at these PTGS target sites is feasible (Froebel and Sullivan, unpublished results, 2017). A similar approach could be applied to autosomal dominant disease genes (e.g., PRPH2 and BEST1) for which there are many mutations and are likely to be candidates for PTGS therapeutics in the future. 

This study is a first-of-kind evaluation of the potential impact mutations have on the accessibility of a single mutation-independent PTGS therapeutic for autosomal dominant retinal degeneration. Although given mutations can increase or decrease accessibility at the target site, overwhelmingly this study has shown that most mutations do not significantly affect accessibility at a single PTGS cleavage site. For some mutations, there is a significant change in predicted accessibility, and for these mutations, experiments should be performed with mutant transcripts to compare in vitro cleavage and in cellulo target knockdown with other mutants and the WT mRNA. To test all known and predicted mutations would be time consuming and cost prohibitive, so a bioinformatics approach could be used to test effects of diverse mutations on lead target site accessibility. However, it may be valuable to test a subset of mutations to obtain experimental data to compare with predicted accessibility. 

Conclusions

A total of 199 mutations/variations of the hRHO gene known or suspected to cause adRP were surveyed for changes in predicted accessibility of the target mRNA, 176 of which were single nucleotide missense mutations. The effect of mutation on mRNA accessibility at the hhRz cleavage site is related to its position within the mRNA sequence and the proximity to the target cleavage site for the lead hhRz. Although many single nucleotide mutations exert little effect on predicted accessibility, some mutations can significantly increase or decrease accessibility when both proximate to and at large distances from the cleavage site. Mutations encompassing more than one nucleotide failed to show a greater chance of having a significant effect on accessibility, perhaps due to an insufficient number of sampled mutations of this sort in proximity to the 725GUC↓ cleavage site. An hRHO adRP mutation is more likely to have an impact on utility of the mutation-independent strategy of hhRz gene therapy if it occurs proximate to the target cleavage site or exerts gross mRNA folding perturbation. The outcomes in this study are likely to follow for other hhRz cleavage sites in hRHO and other targets. Although we found that some mutations can alter the predicted accessibility at a ribozyme cleaving region, overwhelmingly the mutation-independent strategy appears viable to treat diverse mutations in a given disease gene. 

Acknowledgments

The authors thank John (Pei-Wen) Chiang, PhD, FACMG, then Director of the Molecular Diagnostics testing service at the Casey Eye Institute of Oregon Health and Science University for his assistance with the HGMD database. The authors thank Jennifer B. Breen, BS (Sullivan Lab) for critically reading the manuscript. 

Supported, in part, by the Department of Veterans Affairs (VA), Veterans Health Administration, Office of Research and Development (Biomedical Laboratory Research and Development) (VA Merit Grant 1I01BX000669; JMS). JMS is employed, in part, as a Staff Physician-Scientist, Ophthalmology, by VA Western New York. This work was supported, in major part, by National Eye Institute, National Institutes of Health, R01 Grant EY013433 (JMS). This work was supported, in part, by an unrestricted grant from Research to Prevent Blindness to the Department of Ophthalmology/University at Buffalo. This work was supported, in part, by a collaborative SUNY Health Now award. BRF was supported by the State University of New York at Buffalo School of Medicine and the Biomedical Sciences 2014 Summer Research Fellowship Program. AJT was supported, partially, through the Interdisciplinary Science and Engineering Partnership (ISEP) program at the University at Buffalo-SUNY. 

The study was conducted at the Veterans Administration Western New York Healthcare System (Buffalo, NY, USA). Contents do not represent the views of the Department of Veterans Affairs or the US Government. 

Disclosure: B.R. Froebel, None; A.J. Trujillo, None; J.M. Sullivan, P