Absence of Significant Off-Target Splicing Variation with a U7snRNA Vector Targeting DMD Exon 2 Duplications

Abstract

Exon skipping therapies for Duchenne muscular dystrophy that restore an open reading frame can be induced by the use of noncoding U7 small nuclear RNA (U7snRNA) modified by an antisense exon-targeting sequence delivered by an adeno-associated virus (AAV) vector. We have developed an AAV vector (AAV9.U7-ACCA) containing four U7snRNAs targeting the splice donor and acceptor sites of dystrophin exon 2, resulting in highly efficient exclusion of DMD exon 2. We assessed the specificity of splice variation induced by AAV9.U7-ACCA delivery in the Dmd exon 2 duplication (Dup2) mouse model through an unbiased RNA-seq approach. Treatment-related effects on pre-mRNA splicing were quantified using local splicing variation (LSV) analysis. Filtering the transcriptome for differences in treatment-related splicing resulted in only 16 candidate off-target LSVs. Only a single candidate off-target LSV was found in both skeletal and cardiac muscle tissue and occurred at a known variable cassette exon. In contrast, four LSVs represented significant on-target correction of Dmd exon 2 splicing and transcriptome analysis showed correction of known dystrophin-deficient gene dysregulation. We conclude that the absence of off-target splicing induced by treatment with the U7-ACCA vector supports the continued clinical development of this approach.

Background

Exon 2 of the DMD gene, expressed in all 427 kilodalton (kDa) isoforms of the dystrophin protein, is the most commonly duplicated exon associated with Duchenne muscular dystrophy (DMD), accounting for ∼10% of all DMD duplications.^1,2 Single exon duplications are promising targets for molecular therapies, as skipping of a single copy of a duplicated exon would restore a wild-type (WT) reading frame and expression of full-length dystrophin, which may be expected to lead to better outcomes than existing therapies for DMD,^3,4 as well as other approaches currently in clinical trials.⁵

As a platform for testing single exon skipping, we created a novel mouse model that contains a duplicated DMD exon 2 (the Dup2 mouse), which largely recapitulates the findings in the standard mdx mouse.⁶ For any single exon duplication resulting in an out-of-frame transcript, highly efficient skipping of both copies of the exon would be predicted to result in an out-of-frame exon deletion transcript, with no improvement in phenotype. Exon 2, however, is a special case that presents a wide therapeutic window, as complete exclusion of the exon results in the expression of a highly functional dystrophin isoform that results from ribosome assembly on an internal ribosome entry site (IRES) within exon 5 and translational initiation from a methionine within exon 6 of DMD. ^7,8

The extent of this therapeutic window is demonstrated by the exceedingly mild phenotype associated with the first North American founder allele demonstrated in DMD, a nonsense mutation within exon 1 (c.9G>A; p.Trp3X).⁹ Despite the fact that a nonsense mutation would be predicted to result in expression of no dystrophin, most males who carry the mutation only have myalgias with exertion; the most severely affected individual stopped walking at age 62, whereas others had no reported symptoms into their eighth decade.⁹ We demonstrated that this mild phenotype was due to the expression of significant amounts of an isoform of dystrophin resulting from utilization the glucocorticoid-sensitive exon 5 IRES.⁸ The resulting isoform is ∼413 kDa in size and lacks the calponin homology domain 1 (CH1) that makes up the first half of the actin-binding domain 1 (ABD1) at the N-terminus of dystrophin.^7,8 Despite the canonical view that an intact ABD1 is required for dystrophin functionality, the human “experiment”—the presence of the founder allele in humans without significant symptoms, and clearly no effect on reproductive fitness—argues that the actin-binding domain 2 (ABD2), located within the dystrophin rod domain, is sufficient for significant preservation of muscle function.⁸

We noted that duplications of exon 2 were nearly always associated with DMD, even though the duplication resulted in an altered reading frame and hence a premature termination codon in the second copy. At the same time, we noted that deletions of exon 2, which are similarly out of frame, had never been reported in either our large cohort or in the exhaustive dystrophinopathy mutational databases.^1,2,10 Assuming that such a mutation had never been ascertained because of IRES activity, we hypothesized that the IRES was active in the absence of exon 2 but not in the presence of a duplicated exon 2, a result confirmed by in vitro studies, and by the identification of an essentially asymptomatic patient with a deletion of exon 2.⁸ The therapeutic implication of these findings is that there is a wide margin of safety for skipping of exon 2, as boys with exon 2 duplications cannot be “overtreated.” Efficient skipping will result in either the WT transcript containing one copy of exon 2, or in the complete exclusion of exon 2, resulting in IRES activation and the expression of the highly protective isoform lacking CH1.

To induce highly efficient exon 2 skipping, we have developed the scAAV9.U7.ACCA vector, packaging untranslated U7 small nuclear RNAs (U7snRNAs) containing antisense sequences targeting both the splice acceptor and donor sites within a self-complementary adeno-associated virus (AAV) genome, and shown its ability to induce skipping in the Dup2 mouse model of DMD.^6,8 Canonically involved in histone pre-mRNA 3′-end processing, U7snRNA is a component of the U7 small nuclear ribonucleoprotein (snRNP) and contains a Sm-binding sequence (SmWT) responsible for binding the Sm core proteins involved in histone maturation. When this Sm-binding sequence is modified to a consensus sequence (SmOpt) derived from the spliceosomal snRNAs, the snRNPs accumulate more efficiently in the nucleus and become nonfunctional in histone RNA processing. When combined with antisense sequences, the U7snRNA SmOpt is able to compete with the spliceosome and prevent recognition of splice signals and thus prevent inclusion of a specific exon into the mature mRNA.^11
–13 This approach has been used in the mdx mouse, which has a nonsense mutation in exon 23¹⁴; in the mdx/Utr ^−/− double knockout mouse, which also lacks utrophin¹⁵; and in the golden retriever muscular dystrophy dog, which has a splice site mutation in intron 6.¹⁶ In each case, restoration of an open reading frame—although one in which there is an internal deletion in the mRNA—resulted in improved skeletal muscle and/or cardiac function.

The purpose of this study is to identify and quantify potential off-target splicing effects induced by systemic delivery of the scAAV9.U7-ACCA vector in the Dup2 mouse. We delivered clinically relevant doses to P0 mice by a single systemic injection through the facial vein, and 3 months later examined RNA-seq data generated from skeletal muscle, heart, and liver from treated and control Dup2 mice, analyzing gene expression and splicing using an empiric and unbiased approach. We measured on and off-target splicing effects by quantifying “local splice variation” (LSV) at splice junctions defining the start or end of each exon. A single LSV measures the relative level of splice events detected by RNA-seq reads mapping across exon–exon, intron–exon, or exon–intron boundaries. An LSV quantifies the relative amount of splice junction reads on a “Percent Spliced In” (PSI) scale from 0 to 1. For example, an idealized exon skipping event at 20% skipping efficiency would be measured by counting splice junction reads at both the “source” and “target” LSVs. The source LSV is the splice donor site upstream of the skipped exon while the target LSV is the splice acceptor site of the exon downstream of the skipped exon. They will share a common set of RNA-seq junctions reads that skip the intervening exon but will differ in their junction reads to and from the skipped exon. Splicing viewed from the source LSV donor site will count 20% of its RNA-seq junction reads ( = 0.2 PSI value) mapping to the target LSV acceptor site, while 80% of its junction reads (0.8 PSI) will map to the skipped acceptor site. The same exon skipping event will be measured from the target LSV by counting the junction reads in common with the source LSV (0.2 PSI) and the 80% (0.8 PSI) junction reads that map from the skipped donor site to the target LSV acceptor site. Previous RNA-seq studies have shown that quantifying LSVs in this fashion identifies differentially spliced events with high accuracy.¹⁷ LSV analysis shows no evidence for meaningful off-target splicing alteration, whereas significant on-target exon 2 splicing variation induced gene expression changes consistent with a reversal of known dystrophic gene expression signatures, results that provide support for the continued clinical development of this approach.

Methods

Mouse injection

As part of ongoing clinical development studies, male Dup2 mice⁶, were injected through facial vein at P0 with 1.8e14 vg/kg of scAAV9.U7.ACCA, and sacrificed at 3 months for evaluation of exon skipping and dystrophin expression. The full report of that study is under review elsewhere, but an outline of the study is included in Supplemental Data.

RNA isolation and sequencing

Frozen heart, gastrocnemius muscle, and liver tissue from 3 C57/Bl6 WT control, 3 untreated Dup2, and 3 treated Dup2 mice were sectioned on a cryostat and homogenized with a 26-gauge needle in 400 μL of 20 mM Tris*Cl pH 7.4, 150 mM NaCl, 5 mM MgCl₂, 1 mM DTT, and 1% Triton X-100. RNA was extracted using TRIzol^® and precipitated with isopropanol according to the manufacturer's instruction (Life Technologies). Total RNA (∼1 μg) was depleted of mouse cytoplasmic and mitochondrial rRNA using tiling oligodeoxynucleotides and the RNase H digestion method.¹⁸ The Illumina TruSeq Stranded Total RNA library Kits were used to prepare indexed libraries and paired-end, 2 × 125 bp reads were generated on an Illumina HiSeq instrument using v4 chemistry.

RNA-seq read alignment and differential gene expression

Paired-end FASTQ reads were preprocessed for quality trimming and adaptor removal using Trimmomatic¹⁹ with the MAXINFO parameter set to 50:0.2. The effective read lengths of the libraries were extended using the FLASH²⁰ tool that merged overlapping paired-end reads. All reads (extended and nonoverlapping paired-end) were mapped by the STAR RNA-seq aligner²¹ (version 2.6.0b) to the GRCm38/mm10 mouse reference genome with splice junction annotations built using the 1-pass mapping strategy with known junctions from the GENCODE comprehensive gene annotation GTF file (Release M17, https://www.gencodegenes.org). STAR mapping options included a minimum overhang for annotated splice junctions of 1 nt., a minimum overhang for unannotated splice junctions of 8 nts., the –outFilterType = BySJout option to reduce the number of spurious junctions, and the –outFilterMultimapNmax = 10 option as the maximum number of multiple alignments allowed for a read.

Gene-level read counts were extracted from the STAR output Binary Alignment/Map (BAM) files using featureCounts²² with parameters set to strand-specific read counting, filtering secondary alignments, retaining only uniquely mapped reads and counting reads from multioverlapping gene features. Gene-level read count extraction was performed using featureCounts with Ensembl gene_id annotations from GENCODE release M17. Differential expression (DE) was performed using the edgeR Bioconductor package.²³ Read counts were normalized using the trimmed mean of M values method and filtered to remove low-count genes, requiring a count-per-million threshold value of 5 in at least 3 samples. Library size per sample averaged for heart 35.0 million read counts (range 30.2 to 40.3 million), for gastrocnemius muscle 29.1 million read counts (range 22.2 to 32.3 million), and for liver 8.9 million read counts (range 4.5 to 14.8 million). The design matrix and dispersion estimates were fit to a linear model using the edgeR glmFit function and DE was tested with the glmTreat function between treated and control samples relative to a fold change (FC) threshold of 2. Gene-level results were displayed as mean-difference plots, with log₂ FC for each gene plotted against the average abundance in log₂ counts per million mapped reads. Gene Ontology and pathway enrichment analysis was performed with the PANTHER (Protein ANalysis THrough Evolutionary Relationships) Classification System (Version 16.0, released 2020-12-01; www.pantherdb.org).

LSV detection

The Modeling Alternative Junction Inclusion Quantification (MAJIQ/Voila v1.0.5) software¹⁷ was used to detect and quantify splice variation from the filtered bam alignment files. MAJIQ Builder used the GENCODE comprehensive gene annotation file (gencode.vM17.annotation.gff3) downloaded from https://www.gencodegenes.org/mouse/ and the heart, gastrocnemius muscle, and liver samples configured as separate experiments to construct splice graphs and detect LSVs. The splicegraph was generated from the bam alignment files using the MAJIQ build command with the min-denovo parameter = 5, which is the threshold for the minimum total number of reads for a de novo junction to be included. Relative inclusion values (PSI) of splice junctions within LSVs were calculated using the MAJIQ psi command with parameters, minreads = 2, and minpos = 1. The Voila psi command generated tab-delimited output with sample level PSI values and genomic coordinates for exons and splice junctions. Pairwise delta PSI values (ΔPSI) were calculated for rAAV9.U7.ACCA-treated versus untreated groups for each tissue with MAJIQ/Voila deltapsi parameters set to minreads = 5 and minpos = 1. Assignment of mouse LSVs to known alternative splicing (AS) events was achieved by indexing LSV coordinates to an AS event in the EVENTS table (Vertebrate Alternative Splicing and Transcription Database [VastDB] release v1.1, http://vastdb.crg.eu/wiki/Downloads#Mouse) downloaded from the VastDB.²⁴ VastDB contains mouse mm10 annotated events that show a minimal alternative splicing level defined as 0.10 ≤ PSI ≤0.90 in at least 10% of mouse tissues/cells, or by a range of PSIs ≥0.25 across all mouse tissues/cells.

Thermodynamics of off-target interactions

To calculate RNA–RNA-binding energies, we used the ViennaRNA Package 2.0 server²⁵ running the RNAup algorithm (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAup.cgi). The modified U7.ACCA snRNA query sequences used were: Along = 5′-GUUUUCUUUUGAAGAUCUUCUCUUUCAUCUA-3′ and C = 5′-AUUCUUACCUUAGAAAAUUGUGC-3′. The Dmd exon 2 target sequence was from chrX:83136250-83137074 and the Ssbp4 MmuEX0045034 target sequence was from chr8:70599588-70600585, both downloaded from the UCSC genome browser using the GRCm38/mm10 mouse reference genome. For transcriptome-wide analysis of off-target binding energies, we downloaded from the UCSC Table Browser a multi-fasta file with 195,829 entries of RefSeq curated exon sequences with 30 nts. of flanking intronic sequence. The interaction-free energy (ΔGi) between the modified U7 sRNA query sequences and the pre-mRNA targets were calculated using the command line equivalent of “RNAup -b -d2 –noLP -c “S” < sequences.fa > RNAup.out, computed using RNAup 2.4.8.²⁶

Results

Splice junction quantification in treated versus untreated tissues using RNA-seq

To determine treatment-specific mRNA splicing and gene expression differences in the Dup2 mouse, total RNA was purified from gastrocnemius muscle, heart, and liver tissue, using three Dup2 mice from the treated condition (=rAAV9.U7.ACCA virally-mediated Dmd exon 2 skipping) and 3 Dup2 mice from the untreated condition. RNA-seq libraries were sequenced from 18 total tissue samples: 3 treated Dup2 mice × 3 tissues and 3 untreated Dup2 mice × 3 tissues. After quality filtering, the 125 nt. paired-end reads were mapped to the mm10 reference genome of Mus musculus strain C57BL/6J. The average number of uniquely mapped reads per individual sample was: 41.9 × 10⁶ reads (muscle), 49.6 × 10⁶ reads (heart), and 35.9 × 10⁶ reads (liver). An average of 8,401 and 5,547 reads per sample were mapped to the Dmd muscle-isoform transcript (RefSeq NM_007868.6, 79 exons) in skeletal muscle and heart, respectively; while no reads mapped to Dmd exons 1 and 2 in any of the liver samples.

The effect of rAAV9.U7.ACCA treatment on Dmd mRNA splicing is shown in Fig. 1A with RNA-seq read depths and counts from genomic mapping to the mm10 chromosome X reference sequence modified by insertion of the duplicated Dmd exon 2. Dmd exon 1 to exon 3 junction reads were seen only in treated mice, consistent with exon 2 skipping induced by rAAV9.U7.ACCA treatment. Exon 2 skipping results in an increased intron length due to the formation of extended introns spanning the 187 kb intron 1 and 210 kb intron 2 regions and the subsequent delayed degradation of the extended introns due to cotranscriptional splicing followed by rapid degradation of the excised lariat intron.²⁷ The observed increase in intronic depth of coverage surrounding the Dup exon 2 region (Fig. 1B) was consistent with treatment-induced exon 2 skipping and the formation of a fused 397 kb intron between exon 1 and exon 3.

Figure 1.

RNA-seq visualization of alternative spliced exons for Dmd transcripts in skeletal muscle and heart. (A) Aligned RNA-seq reads displayed in the Broad Integrated Genome Viewer with depth of coverage plotted on the y-axis and the Dmd genomic coordinates (mm10) on the x-axis. Merged RNA-seq data from 3 animals is shown for the untreated/treated conditions in skeletal muscle (blue) and heart tissue (brown), and the y-axis scale ranges from 0 to 1082 (untreated skeletal muscle), 612 (treated heart), 841 (treated skeletal muscle), and 523 reads (treated heart), respectively. The exon–intron locations are shown for RefSeq transcript NM_007868.6 and an arrow points to the location where the duplicated exon 2 sequence was inserted into the mm10 genomic reference sequence. Junction reads are plotted as arcs connecting spliced exons and the number of junction-spanning reads is indicated for each arc. (B) Depth of coverage plot showing the relative increase in the intronic depth of coverage for the region between exon 2 and Dup2. (C) Junction read counts from the best alignment to different Dmd exon 1:2:Dup2:3 transcript isoforms shown for each condition and tissue (right) and each transcript isoform (left).

Since the genomic mapping of RNA-seq reads requires splitting of reads by the mapping algorithm, we extracted reads that mapped to Dmd exons 1, 2, or 3 and then realigned them directly to transcript sequences representing exon 1-3, exon 1-2-3 and exon 1-2-2-3 isoforms. Figure 1C shows the number of reads that completely span exon–exon junctions for exons 1-3 (skipping both copies of exon 2), exons 1-2-3 (skipping of a single exon 2), and exons 1-2-2-3 (the expected spliced transcript from the Dup2 mutation), plus reads that cannot be unambiguously assigned to an isoform but bridge the exon 1-2 or 2-3 junction. In both tissues, the level of exon 2 skipping indicated that the virally expressed modified U7snRNAs were active at their primary Dmd exon 2 targets and confirmed the general validity of quantifying treatment-induced alternative splicing from the mapped RNA-seq reads.

Global analysis of differential mRNA splicing associated with rAAV9.U7.ACCA treatment

To quantify on- versus off-target alternative splicing events from global RNA-seq mapping data, we used the MAJIQ software to detect differential splicing between groups of replicates from the rAAV9.U7.ACCA-treated versus untreated conditions. MAJIQ used the RNA-seq genomic alignment (BAM) files and transcriptome annotations to quantify LSVs, which are the ratios of exon junction-spanning reads from a single “source” (splice donor) or “target” (splice acceptor) exon. An example of MAJIQ output for Dmd exon 2 LSV-centric splice junction analysis is shown in Fig. 2A. Skipping of exon 2 leads to multiple splice junction reads from the Dmd exon 1 splice donor site and the exon 3 splice acceptor site. Alternative splicing was quantified by the PSI values at the exon 1 “source” LSV for exons 1:2 and 1:3 junction reads, and the exon 3 “target” LSV for exons 2:3 and 1:3 junction reads.

Figure 2.

LSVs using MAJIQ. (A) The constitutive splice junctions for Dmd exon 1 and exon 3 are connected by blue lines, with source LSV splice junction reads spanning exon 1-2 (upper panel) and the target LSV splice junction reads spanning exon 2-3 (lower panel). The splice junction reads spanning exon 1-3 produced by exon 2 skipping are connected by green lines for the source LSV and red lines for the target LSV. (B) Violin plots indicate the posterior distribution of estimated PSI values for each junction in the treated and untreated condition. The constitutive exon 1:2 and 2:3 junctions are colored in blue, while the induced exon 1:3 junctions are colored in green (source LSV) or red (target LSV) and with their estimated PSI values. The ΔPSI violin plots correspond to the posterior distribution of differential splicing for each splice junction, with positive values indicating increased differential inclusion in treated tissues. ΔPSI, delta percent spliced in; LSV, local splicing variation; MAJIQ, Modeling Alternative Junction Inclusion Quantification.

In untreated skeletal muscle and heart tissue, PSI values >0.99 indicated that exon 1:2 and exon 2:3 junctions were constitutively spliced and Dmd exon 2 skipping (exon 1:3 junctions) was not detected (PSI value threshold <0.01) in the untreated condition from either source or target LSVs. The rAAV9.U7.ACCA treatment altered the exon 1:2, exon 2:3, and exon 1:3 junction PSI values in both tissues (Fig. 2B). Decreased constitutive splicing was supported by reduced PSI values for the exon 1:2 junctions from 1.00 to 0.84 (muscle) or 0.55 (heart) and exon 2:3 junctions from 1.00 to 0.83 (muscle) or 0.52 (heart). Treatment-induced exon 2 skipping was supported by increased PSI values for the exon 1:3 junctions from 0.00 to 0.16 (muscle) or 0.54 (heart), and the relative change in LSV differential inclusion (ΔPSI) between the two conditions (Fig. 2B) was quantified as 0.14 in muscle (both source and target LSVs), and 0.44 (source LSV) or 0.47 (target LSV) in heart.

Since the LSV methodology was capable of detecting and quantifying treatment-induced alternative splicing at Dmd exon 2, we cataloged transcriptome-wide LSVs in skeletal muscle, heart, and liver tissues. Differences in transcriptome complexity in each tissue led to variation in the total number of LSVs detected, where a candidate LSV was defined by two or more splice junctions emanating from a single splice donor site (source) or splice acceptor site (target). We used two sequential data filters to narrow the search for treatment-induced differences in constitutive splice junction usage. The first filter selected for constitutive splice sites in the untreated condition by requiring PSI values >0.99. This threshold is level of constitutive splicing observed for Dmd exon 2 in the untreated condition. In heart tissue, 35,880 constitutive LSVs were detected in 9,489 genes; in skeletal muscle, 29,606 constitutive LSVs were detected in 7,964 genes; and in liver, 8,084 constitutive LSVs were detected in 3,113 genes. The second filter required a minimum ΔPSI value of 0.10 between the treated versus untreated condition at these constitutive LSVs.

Only 20 candidate LSVs mapping to 12 genes exceeded a ΔPSI ≥0.10 threshold difference: 11 LSVs from heart, 8 LSVs from skeletal muscle, and 1 LSV from liver (Table 1). Four of these candidate LSVs were the on-target effects for Dmd exon 2 skipping. When ranked by highest ΔPSI values, the Dmd exon 2 source and target LSVs ranked 1 and 2 in heart tissue, and 5 and 6 in skeletal muscle. Five of these candidate LSVs were not splicing LSVs but instead represented alternative promoter LSVs associated with upregulation of promoter activity in the treated samples. Figure 3 shows two examples of upregulation of alternative promoter LSVs for the Tacc2 (transforming, acidic coiled-coil-containing protein 2) gene in heart tissue, and the Abacb (acetyl-Coenzyme A carboxylase beta) gene in liver tissue.

Figure 3.

Differences in alternate promoter usage in candidate LSVs. Merged RNA-seq data from 3 animals is shown for the untreated/treated conditions as described in Figure 1, with the scale for RNA-seq read depth indicated by the bracketed numbers. The yellow boxes highlight the alternate promoter junction reads that are unique to the treated condition. (A) RNA-seq coverage plot of splicing patterns in heart tissue from multiple alternate promoters in the 5′ end of the Tacc2 gene. The arrow indicates the location of the target LSV (chr7:130577438–130577516) with a ΔPSI value of 0.15 in treated versus untreated samples. (B) Splicing patterns in liver from two alternate promoters in the 5′ end of the Acacb gene. The arrow indicates the location of the target LSV (chr5:114165509–114166140) with a ΔPSI value of 0.14 in treated versus untreated samples.

Table 1.

Treatment-induced differences in splice junction usage from constitutive exons

Tissue Event	Gene; LSV Name	Event Type (VastDB ID)	ΔPSI Junction #1	ΔPSI Junction #2
Heart 1	Dmd; chrX:83346262–83346354;target	Targeted exon 2	0.47	−0.44
Heart 2	Dmd; chrX:82948870–82949147;source	Targeted exon 2	0.44	−0.44
Heart 3	Ssbp4;chr8:70600396–70600485;target	Known alt. splice: MmuEX0045034	0.25	−0.23
Heart 4	Rnf185;chr11:3432400–3432620;source	Alt. promoter	0.23	−0.18
Heart 5	Tdp1;chr12:99935005–99937362;target	Alt. promoter	0.21	−0.21
Heart 6	Ssbp4;chr8:70599688–70599757;source	Known alt. splice: MmuEX0045034	0.21	−0.19
Heart 7	Tbcel;chr9:42451613–42451762;target	Known alt. splice: MmuEX0046565	0.19	−0.19
Heart 8	Fhl3;chr4:124705613–124705789;target	Alt. promoter	0.18	−0.18
Heart 9	Tbcel;chr9:42444407–42444588;source	Known alt. splice: MmuEX0046565	0.17	−0.17
Heart 10	Ppp3cc;chr14:70266630–70266754;target	Known alt. splice: MmuEX0036802	0.16	−0.16
Heart 11	Tacc2;chr7:130577438–130577516;target	Alt. promoter	0.15	−0.1
SkMuscle 1	Ssbp4;chr8:70600396–70600485;target	Known alt. splice: MmuEX0045034	0.3	−0.28
SkMuscle 2	Ssbp4;chr8:70599688–70599757;source	Known alt. splice: MmuEX0045034	0.19	−0.19
SkMuscle 3	Fam204a;chr19:60221054–60221304;target	Known alt. splice: MmuEX0013516	0.18	−0.22
SkMuscle 4	Rian;chr12:109607414–109607496;source	Known alt. splice: MmuEX0002672	0.15	−0.15
SkMuscle 5	Dmd;chrX:82948870–82949147;source	Targeted exon 2	0.14	−0.14
SkMuscle 6	Dmd;chrX:83346262–83346354;target	Targeted exon 2	0.14	−0.14
SkMuscle 7	Fam204a;chr19:60220330–60220360;source	Known alt. splice: MmuEX0013516	0.14	−0.14
SkMuscle 8	Sobp;chr10:43020982–43022918;source	Known alt. splice: MmuEX0044325	0.13	−0.13
Liver 1	Acacb;chr5:114165509–114166140;target	Alt. promoter	0.14	−0.14

Bold marks targeted exon 2 sequence junctions. ΔPSI, delta percent spliced in; LSV, local splicing variation; Ssbp4, single-stranded DNA-binding protein 4; VastDB, Vertebrate Alternative Splicing and Transcription Database.

Analysis of off-target effects for the remaining candidate splicing LSVs

The 11 remaining candidate off-target LSVs mapped to 6 genes and were all known alternative splicing events annotated in the VastDB, a large annotated resource of alternative splicing events detected by RNA-seq from a diverse set of vertebrate cells and tissues.¹⁰ Since we observed on-target exon 2 skipping of Dmd transcripts in both skeletal muscle and heart, we assumed that ΔPSI differences due to off-target U7.ACCA binding would be observed in both tissues if the pre-mRNA targets were detectable by RNA-seq. For each of these six genes with candidate off-target LSVs, transcript levels were nearly equivalent in skeletal muscle and heart (Figs. 4 and 5). Two LSVs in Tbcel and one in Ppp3cc exceeded the ΔPSI ≥0.10 threshold difference in heart only, whereas two LSVs in Fam204a and one each in Rian and Sobp exceeded the ΔPSI threshold in skeletal muscle only. Failure to observe ΔPSI differences in both tissues for Tbcel, Ppp3cc, and Rian LSVs (Fig. 4A, B) suggests that off-target binding effects were not responsible for these splicing changes. For Fam204a and Sobp, the alternative splicing pattern seen in treated muscle is similar to the splicing pattern seen in both treated and untreated heart (Fig. 4C) and this splicing pattern is also not consistent with off-target binding stimulating the alternative splice.

Figure 4.

Tissue-specific AS candidate LSVs. Merged RNA-seq data from three animals are shown for the untreated/treated conditions as described in Figure 1, with the scale for RNA-seq read depth indicated by the bracketed numbers. The shaded green boxes indicate the location of candidate LSVs shown in Table 2 and the yellow boxes highlight the AS junction reads that are responsible for the detected LSVs. (A) Splicing patterns from treated heart-only multijunction LSVs, (B) treated skeletal muscle-only multijunction LSVs, and (C) untreated skeletal muscle-only single-junction LSVs. AS, alternative splicing.

Figure 5.

Off-target analysis for the Ssbp4 candidate LSVs. (A) Merged RNA-seq data from three animals are shown for the untreated/treated conditions as described in Figure 1, with the scale for RNA-seq read depth indicated by the bracketed numbers. The shaded green boxes indicate the location of candidate LSVs shown in Table 2, and the yellow boxes highlight the alternative splice junction reads responsible for the detected LSVs (B) The plot shows the ΔGi interaction-free energy for the interaction of U7.Along (blue) or C (red) snRNA with pre-mRNA from the target region (Sspb4 or Dmd exon 2) and the energy needed to open existing structures in the target sequence (grey). Ssbp4, single-stranded DNA-binding protein 4.

In addition to Dmd exon 2 skipping, the remaining 4 LSVs in Ssbp4 were the only example of ΔPSI differences seen in both skeletal muscle and heart (Fig. 5A). In VastDB, these Ssbp4 LSVs mapped to a known mouse cassette exon (VastDB ID MmuEX0045034) whose inclusion/exclusion leads to alternative protein isoforms of single-stranded DNA-binding protein 4 (Ssbp4), which is alternatively spliced in multiple tissues, including heart, brain, lung, and adipose tissue. To identify putative off-target binding sites in Ssbp4 pre-mRNA, we computed the predicted binding energy between U7.A or C snRNA and the Ssbp4 pre-mRNA using the RNAup algorithm,¹² which computes the interaction energy of small RNA—target RNA hybrids combined with the probability that potential binding sites are accessible on the target pre-mRNA (not paired in secondary structure). Figure 5B shows the binding energy profiles for the U7.A or C snRNA interactions for every position along the target pre-mRNA sequences from the Ssbp4 MmuEX0045034 region and Dmd exon 2 region. While the RNA–RNA interaction scan detected an optimal binding energy overlapping the Dmd exon 2 splice sites, the Ssbp4 profile gave no indication of a strong attractive interaction between the U7.A or C snRNA and Ssbp4 pre-mRNA. This suggests that the Ssbp4 ΔPSI differences are an indirect effect of treatment on alternative splicing instead of a direct off-target effect due to modified U7 snRNA—target pre-mRNA binding.

Since treatment-induced Dmd exon 2 skipping was greater in heart versus skeletal muscle, we also evaluated off-target binding sites predicted by the RNAup algorithm for the heart-specific LSV events in Ppp3cc and Tbcel. The skipped Ppp3cc exon did not show putative binding sites for Along or C duplex formation, but the 140 nt. Tbcel skipped exon had an RNAup-predicted 13 base pair Along duplex at positions +59 to +71 within the exon (Supplementary Fig. S1A). The interaction-free energy (ΔGi) of this putative off-target site was −13.8, compared with a predicted ΔGi of −26.2 for the on-target Along Dmd exon 2 splice acceptor site. To evaluate how frequent putative binding sites below this interaction-free energy may occur in the transcriptome, we examined the distribution of putative Along-binding sites in 195,829 annotated exons from the mm10 reference mouse genome. For each exon plus 30 nts. of flanking intron sequence, we predicted the lowest interaction energy for that exon and then rank ordered this list of sites. The on-target Along Dmd exon 2 splice acceptor site was at the top of the list, while the ΔGi of the Tbcel site was seen at rank 1,581 or within the top 1% of sites. The putative off-target site in Tbcel maps to exon 3 of the gene, and reverse transcription PCR (RT-PCR) analysis of Tbcel mRNA from heart shows low-level skipping (0.8–2.1% of transcript) in WT BL6 mice, which increased to 20.6–23.0% in treated mice (Supplementary Fig. S2), comparable to results from RNA-seq.

Differential gene expression associated with rAAV9.U7.ACCA treatment

We also analyzed gene expression to evaluate whether restoration of dystrophin in skeletal muscle and heart led to specific transcriptome-wide changes. In skeletal muscle, there were 16 DE genes observed at a 5% false discovery rate (Fig. 6A and Supplementary Table S1). Upregulation of insulin-like growth factor 2 (Igf2) and sarcomeric components for the embryonic myosin heavy chain genes (Myh3 and Myh8) and troponin T2 (Tnnt2) are consistently observed in both mdx mouse²⁸ and human²⁹ studies of dystrophic skeletal muscle tissue. Expression of these four genes were restored in the treated condition to levels observed in skeletal muscle from C56BL/6J control mice (Fig. 6B). In a comparative protein and mRNA-level omics study using the mdx mouse,³⁰ the prune homolog 2 (Prune2) was among the 10 most upregulated genes in dystrophic muscle and Prune2 is downregulated in the treated condition. A similar but opposite effect was seen for the protein phosphatase 1, regulatory subunit 1A (Ppp1r1a) and the dual specificity phosphatase 26 (Dusp26) genes, which are downregulated in the mdx mouse, but upregulated in the treated condition (Fig. 6B).

Figure 6.

Gene-level expression responses associated with treatment. (A) Volcano plot of skeletal muscle gene expression differences from the treated versus untreated conditions with the statistical significance (log₁₀-FDR, false discovery rate) shown on the y-axis and the magnitude of change (log₂-FC) shown on the x-axis. The labeled genes (red points) are differentially expressed in treated versus untreated skeletal muscle at an FDR cutoff of 0.05 and have a log₂-FC >2, while the green points are genes that have only a log₂-FC >2. The heat map columns show the relative expression levels of the labeled genes (columns) passing both the FDR and FC threshold, with rows for each animal in the treated versus untreated conditions, and in control BL/6 mice. The dendrogram indicates the hierarchical clustering of both rows and columns. (B) Volcano plot and heat map as in (A) for gene expression in treated versus untreated heart. FC, fold change; FDR, false detection rate.

There were 30 DE genes observed in treated versus untreated heart (Fig. 6C and Supplementary Table S2) at a 1% false discovery rate, not including U7snRNA (Rnu7) reads that mapped to the precursor RNA from virally delivered U7.ACCA snRNA gene. Unlike skeletal muscle, clustering of the heart samples revealed that the DE transcripts were specifically changed in the treated condition compared with untreated and control C57BL/6J tissue (Fig. 6D), with 26 of the 30 DE transcripts upregulated in treated samples. Beta-1,4-N-acetyl-galactosaminyl transferase 2 (B4galnt2, formerly Galgt2) is highly upregulated (Fig. 6C) and overexpression of B4galnt2 is known to have therapeutic effects mediated by glycosylation of α dystroglycan in dystrophic skeletal muscle.³¹ Gene Ontology enrichment analysis for biological process terms indicated overrepresentation of genes annotated for positive regulation of protein phosphorylation (Gdf15, Fgf21, Nrg1, Gdf10, and Crlf1; GO:0001934, false discovery rate (FDR) = 7.7E-03; fold enrichment = 16) and genes for carboxylic acid metabolic process (Phgdh, Psat1, Asns, Acsm5, and Mthfd2; GO:0019752, FDR = 3.9E-02; fold enrichment = 8). Pathway analysis suggested enrichment for genes involved in the serine glycine biosynthesis pathway (PANTHER pathway accession = P02776; FDR = 2.4E-02; fold enrichment >100) that converts 3-phosphoglycerate into serine with the first two steps catalyzed by phosphoglycerate dehydrogenase (Phgdh) and phosphoserine aminotransferase 1 (Psat1), both of which are upregulated in the treated condition. Genes in the serine (Phgdh and Psat1) and one-carbon metabolism (Mthfd2) pathway are known to be regulated by the transcription factor Atf4,³² which along with Atf5 is upregulated in treated heart (Fig. 6C, D and Supplementary Table S2).

Discussion

The development and deployment of novel therapies directed toward alteration of mRNA splicing understandably requires studies to address the specificity of exon splicing effects. To address this question in the case of the rAAV9.U7ACCA vector, we chose an unbiased approach using a global RNA-seq analysis for several reasons. First, in silico analysis failed to identify regions of significant homology to the exon 2 target sequences that would serve as likely candidates for altered splicing amenable to a targeted RT-PCR approach. The second reason, a corollary to first, is that global RNA-seq analysis allows an unbiased assessment of all transcripts to search for unanticipated splicing events, or unannotated regions of target homology. Third, RNA-seq allows us to directly quantify the target DMD mRNA exon 2 junction sequences. Finally, RNA-seq allows characterization of expression of the entire transcriptome, confirming alterations known to be associated with muscle from DMD patients or DMD models and quantifying the degree of correction following delivery of the vector.

Among these, the most critical for assessing safety is identification of off-target and unanticipated splicing events. Using conservative quantitative thresholds to quantify the degree of LSV among all exon/exon junctions, our data confirmed altered splicing of the intended DMD exon 2 target junctions, which may be interpreted as a control for the sensitivity of the approach. We note that alteration of exon 2 splicing was more pronounced in the heart than in skeletal muscle, a finding consistent with greater vector tropism in heart than in skeletal muscle in non-human primates treated with the same vector.³³ At the same time, these data demonstrate no convincing evidence for off-target splicing alterations. The only gene that showed significant alteration of splicing in both cardiac and skeletal muscle was Ssbp4, where LSVs were found at a cassette exon location that is already known to be widely alternatively spliced. An increase in endogenous alternative splicing in Tbcel could be confirmed in heart by RT-PCR, but is unlikely to be of biological significance. The resulting Tbcel transcript isoform ablates the tubulin-specific chaperone cofactor E-like protein reading frame leading to a short N-terminal peptide, but ∼80% of the transcripts remained WT. While this study detected splice junctions in the large number of genes, a limitation is the detection of transcripts expressed at low levels, which will require substantially greater read depths for a more complete survey of transcriptomes in target tissues and cell types.

In evaluating the transcriptome, the pattern of DE of genes in the tissues studied suggested treatment-related, specific transcriptional responses as an indirect effect of dystrophin restoration. In dystrophic skeletal muscle from humans as well as animal models, the upregulation of the embryonic/developmental myosin heavy chain genes are robust markers of regenerating muscle fibers. The reduction in Myh3 and Myh8 transcript levels observed in skeletal muscle after rAAV9.U7.ACCA treatment is similar to results seen after restoration of dystrophin by oligonucleotide-mediated exon 23 skipping in the mdx mouse.³⁴ Also, the increase in Ppp1r1a and Dusp26 transcript levels reverse their downregulation observed in dystrophic skeletal muscle, suggesting that dystrophin expression from virally mediated exon skipping restored the skeletal muscle transcriptome.

In contrast, the transcriptional signature in the heart was unique to the treated condition, differing from both untreated Dmd Dup2 and control C57BL/6J tissue. Many of the top 30 DE genes are direct targets for transcriptional regulation by the stress-induced activating transcription factors Atf4 and Atf5 in response to mitochondrial perturbations, including Fgf21, Asns, and Trib3. In a mouse model of mitochondrial myopathy,³⁵ deletion of mtDNA induces an integrated stress response and upregulation of Atf4 and Atf5, which in turn upregulate genes that have amino acid starvation response elements (AARE) in their upstream regulatory regions. We observed treatment-associated upregulation of both Atf5 and Atf4, and of genes known to be controlled by the AARE response, including: Fgf21, Mthfd2, Gdf15, Phgdh, and Psat1. Enrichment analysis of the treatment-induced upregulated genes in heart implicated the serine glycine biosynthesis (Phgdh and Psat1) and one-carbon metabolism pathways (Mthfd2, the mitochondrial folate-coupled dehydrogenase). Previous studies with mice overexpressing the CnAβ2 (Ppp3cb) catalytic subunit of the calcium-regulated phosphatase calcineurin have observed ATF4-dependent upregulation of Phgdh, Psat1, Mthfd2, and Aldhl2 specifically in cardiac tissue, and it has been proposed that this transcriptional response may have a protective effect on cardiac function mediated by increased antioxidant production.³⁶ The calcineurin gamma catalytic subunit encoded by Ppp3cc that we observed with a heart-specific exon skipping LSV event may be linked to this calcineurin-mediated upregulation of the serine glycine pathway genes. It is noteworthy that besides the four differentially spliced LSVs detected at the targeted Dmd exon 2, all the other 16 candidate LSVs occurred at known alternative splice junctions. It is possible that the observed differential splicing at these LSVs were also due to an indirect effect of dystrophin restoration.

Altogether, these results demonstrate that systemic delivery of the rAAV9.U7.ACCA vector does not result in nonspecific splicing alteration, findings that are entirely consistent with a similar study using a U7snRNA vector directed toward exon 53.³⁷ Rather than nonspecific splicing, delivery of rAAV9.U7.ACCA results in targeted exon 2 exclusion, along with altered gene expression from a disease pattern to that of WT muscle. We note although these results do not by themselves constitute the complete preclinical safety analysis—a pathway that included, for example, testing of vector safety in non-human primates³³—these data are strongly supportive of the rationale for use of systemic U7snRNA therapy, and served as preclinical data leading to initiation of an ongoing gene therapy trial using the vector we describe (NCT04240314).

Footnotes

Acknowledgments

The authors acknowledge the technical assistance of Hannah Huang.

Author Disclosure

Drs. K.M.F. and N.W. hold patents related to the scAAV9.U7 vector, and both they and Nationwide Children's Hospital receive royalties from the licensing of that patent to Audentes Therapeutics.

Funding Information

This work was supported by the Beauhawks Foundation (KMF).

Supplementary Material

Supplementary Data

Supplementary Table S1

Supplementary Table S2

Supplementary Figure S1

Supplementary Figure S2

References

White

, Aartsma-Rus

, Flanigan

, et al. Duplications in the DMD gene. Hum Mutat, 2006; 27:938–945.

Flanigan

, Dunn

, von Niederhausern

, et al. Mutational spectrum of DMD mutations in dystrophinopathy patients: application of modern diagnostic techniques to a large cohort. Hum Mutat, 2009; 30:1657–1666.

Charleston

, Schnell

, Dworzak

, et al. Eteplirsen treatment for Duchenne muscular dystrophy: exon skipping and dystrophin production. Neurology, 2018; 90:e2146–e2154.

Heo

. Golodirsen: first approval. Drugs, 2020; 80:329–333.

Duan

. Systemic AAV micro-dystrophin gene therapy for Duchenne muscular dystrophy. Mol Ther, 2018; 26:2337–2356.

Vulin

, Wein

, Simmons

, et al. The first exon duplication mouse model of Duchenne muscular dystrophy: a tool for therapeutic development. Neuromuscul Disord, 2015; 25:827–834.

Gurvich

, Maiti

, Weiss

, et al. DMD exon 1 truncating point mutations: amelioration of phenotype by alternative translation initiation in exon 6. Hum Mutat, 2009; 30:633–640.

Wein

, Vulin

, Falzarano

, et al. Translation from a DMD exon 5 IRES results in a functional dystrophin isoform that attenuates dystrophinopathy in humans and mice. Nat Med, 2014; 20:992–1000.

Flanigan

, Dunn

, von Niederhausern

, et al. DMD Trp3X nonsense mutation associated with a founder effect in North American families with mild Becker muscular dystrophy. Neuromuscul Disord, 2009; 19:743–748.

10.

Bladen

, Salgado

, Monges

, et al. The TREAT-NMD DMD Global Database: analysis of more than 7,000 Duchenne muscular dystrophy mutations. Hum Mutat, 2015; 36:395–402.

11.

De Angelis

, Sthandier

, Berarducci

, et al. Chimeric snRNA molecules carrying antisense sequences against the splice junctions of exon 51 of the dystrophin pre-mRNA induce exon skipping and restoration of a dystrophin synthesis in Delta 48–50 DMD cells. Proc Natl Acad Sci U S A, 2002; 99:9456–9461.

12.

Brun

, Suter

, Pauli

, et al. U7 snRNAs induce correction of mutated dystrophin pre-mRNA by exon skipping. Cell Mol Life Sci, 2003; 60:557–566.

13.

Goyenvalle

. Engineering U7snRNA gene to reframe transcripts. Methods Mol Biol, 2012; 867:259–271.

14.

Goyenvalle

, Vulin

, Fougerousse

, et al. Rescue of dystrophic muscle through U7 snRNA-mediated exon skipping. Science, 2004; 306:1796–1799.

15.

Goyenvalle

, Babbs

, Wright

, et al. Rescue of severely affected dystrophin/utrophin-deficient mice through scAAV-U7snRNA-mediated exon skipping. Hum Mol Genet, 2012; 21:2559–2571.

16.

Le Guiner

, Montus

, Servais

, et al. Forelimb treatment in a large cohort of dystrophic dogs supports delivery of a recombinant AAV for exon skipping in Duchenne patients. Mol Ther, 2014; 22:1923–1935.

17.

Vaquero-Garcia

, Barrera

, Gazzara

, et al. A new view of transcriptome complexity and regulation through the lens of local splicing variations. Elife, 2016; 5:e11752.

18.

Adiconis

, Borges-Rivera

, Satija

, et al. Comparative analysis of RNA sequencing methods for degraded or low-input samples. Nat Methods, 2013; 10:623–629.

19.

Bolger

, Lohse

, Usadel

. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics, 2014; 30:2114–2120.

20.

Magoc

, Salzberg

. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics, 2011; 27:2957–2963.

21.

Dobin

, Davis

, Schlesinger

, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 2013; 29:15–21.

22.

Liao

, Smyth

, Shi

. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics, 2014; 30:923–930.

23.

McCarthy

, Chen

, Smyth

. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res, 2012; 40:4288–4297.

24.

Tapial

, Ha

KCH

, Sterne-Weiler

, et al. An atlas of alternative splicing profiles and functional associations reveals new regulatory programs and genes that simultaneously express multiple major isoforms. Genome Res, 2017; 27:1759–1768.

25.

Lorenz

, Bernhart

, Honer Zu Siederdissen

, et al. ViennaRNA package 2.0. Algorithms Mol Biol, 2011; 6:26.

26.

Muckstein

, Tafer

, Hackermuller

, et al. Thermodynamics of RNA-RNA binding. Bioinformatics, 2006; 22:1177–1182.

27.

Ameur

, Zaghlool

, Halvardson

, et al. Total RNA sequencing reveals nascent transcription and widespread co-transcriptional splicing in the human brain. Nat Struct Mol Biol, 2011; 18:1435–1440.

28.

Porter

, Khanna

, Kaminski

, et al. A chronic inflammatory response dominates the skeletal muscle molecular signature in dystrophin-deficient mdx mice. Hum Mol Genet, 2002; 11:263–272.

29.

Haslett

, Sanoudou

, Kho

, et al. Gene expression comparison of biopsies from Duchenne muscular dystrophy (DMD) and normal skeletal muscle. Proc Natl Acad Sci U S A, 2002; 99:15000–15005.

30.

Roberts

, Johansson

, McClorey

, et al. Multi-level omics analysis in a murine model of dystrophin loss and therapeutic restoration. Hum Mol Genet, 2015; 24:6756–6768.

31.

Thomas

, Xu

, Martin

. B4GALNT2 (GALGT2) gene therapy reduces skeletal muscle pathology in the FKRP P448L mouse model of limb girdle muscular dystrophy 2I. Am J Pathol, 2016; 186:2429–2448.

32.

Ben-Sahra

, Hoxhaj

, Ricoult

SJH

, et al. mTORC1 induces purine synthesis through control of the mitochondrial tetrahydrofolate cycle. Science, 2016; 351:728–733.

33.

Gushchina

, Frair

, Rohan

, et al. Lack of toxicity in non-human primates receiving clinically relevant doses of an AAV9.U7snRNA vector designed to induce DMD exon 2 skipping. Hum Gene Ther 2021 [Epub ahead of print]; DOI: 10.1089/hum.2020.286.

34.

Guiraud

, Edwards

, Squire

, et al. Embryonic myosin is a regeneration marker to monitor utrophin-based therapies for DMD. Hum Mol Genet, 2019; 28:307–319.

35.

Forsstrom

, Jackson

, Carroll

, et al. Fibroblast growth factor 21 drives dynamics of local and systemic stress responses in mitochondrial myopathy with mtDNA deletions. Cell Metab, 2019; 30:1040–1054 e1047.

36.

Padron-Barthe

, Villalba-Orero

, Gomez-Salinero

, et al. Activation of serine one-carbon metabolism by calcineurin abeta1 reduces myocardial hypertrophy and improves ventricular function. J Am Coll Cardiol, 2018; 71:654–667.

37.

Domenger

, Allais

, Francois

, et al. RNA-seq analysis of an antisense sequence optimized for exon skipping in Duchenne patients reveals no off-target effect. Mol Ther Nucleic Acids, 2018; 10:277–291.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.02 MB

0.49 MB

2.27 MB

0.99 MB

1.18 MB