Differential Expression of mRNA Splicing Isoforms Generates Transcript Length Variants in Atherosclerosis Progression

1. Background

The cell nucleus is pervasively transcribed into a plethora of different RNA molecules with different coding, structural or regulatory roles [1] for a recent review), among them microRNAs. These are small RNAs (over 20 nucleotides long) that exert a post-transcriptional control of gene expression by promoting mRNA degradation or by repressing translation through their binding to homologous sequences mainly at the 3’ Untranslated Regions (3’UTRs) of mRNAs [2]. 3’UTRs not only protect mRNAs from exonucleolytic digestión but also provide functional platforms for the miRNA regulation of mRNA function [3,4], which furthermore adapt to new regulatory requirements by changing the lengths of their sequences and the arrangement of regulatory motifs they harbour. In this way, sequence changes in 3’UTRs results in a deregulation of gene expression, linked to the loss (or acquisition) of specific miRNA binding sites, that eventually could lead to diseases ([5] for a review).

3’UTRs are highly polymorphic in length and sequence, and variant 3’UTRs provide mRNAs with different binding sites for miRNAs or RNA binding proteins (RBPs) which would contribute to the establishment of new regulatory environments [6]. Regulation of 3’UTR length is thus critical for the control of gene expression by its potential to regulate accessibility of miRNAs or RBPs to their cognate sequences in mRNAs. In this way, mRNA isoforms with alternative 3’UTRs has been seen to have differential mRNA stability [7], microRNA binding potential [8], or interactions with ceRNAs [9]. Alternative polyadenylation is one important mechanism for regulating 3’UTR length [10]. Pre-messenger RNAs (pre-mRNAs) harbour specific polyadenylation signals (AAUAAA in its canonical form) over 30 nts upstream the end of their 3’UTRs, to precisely position the cleavage/polyadenylation complex (CPSF) next the cleavage site where PABPN will extend the poly-A tail [11]. Over 70% of human genes have more than one polyadenylation site in their 3’UTRs and 50% have three or more [12], thus making the generation of alternative 3’UTRs by alternative polyadenylation (APA) a widespread mechanism for the generation of transcript variants that are heterogeneous in length and have different regulatory potentials [13,14]. On the other hand, alternative splicing has been shown to be also involved in the regulation of 3’UTR lengths, and their regulatory potentials, through different mechanisms such as intron retention, exon skipping, incorporation of one of two mutually-exclusive terminal exons of different length, use of 5’/3’ alternative sites or activation of cryptic splice sites [3,15].

The emerging picture on alternative 3’UTRs highlights a complex collection of isoforms whose individual expression switches in a controlled way during cell proliferation and differentiation, but that could be de-regulated under stress conditions. This would impact on the RNA interactome, causing a miscontrol of mRNA function, and facilitating the loss of physiologic homeostasis. In this way, 3’UTR shortening would increase mRNA stability by relaxing protein or miRNA-based mechanisms of mRNA degradation, while 3'UTR lengthening would strengthen accessibility to miRNAs [16].

On the other hand, 5’UTRs are also heterogeneous in length and play a key role in regulating translation efficiency and mRNA stability [17] since they contain a number of regulatory elements and miRNA binding sites [18], piRNA binding sites [19], translation-regulating secondary structures [20] and protein binding sites [21]. See [18] for a recent review on the impact of 5’UTR length variations in cell physiology.

Our group is interested in studying the disease-dependent alterations in the regulatory circuitry controlling mRNA stability and function, more specifically in the interactions among miRNAs/sponging ceRNAs and mRNAs. In a previous work, we described the association of adverse cardiovascular events in hemodialysis patients to SNPs in the gene of the lncRNA ANRIL [22]. Furthermore, we took advantage of the atherosclerosis-prone ApoE^-/- mouse model to investigate the role of CD40 signaling on atherosclerosis (ATS) progression. Mice treated with an anti-CD40 siRNA for 16 weeks showed a significantly reduced extension and severity of atherosclerotic lesions, as well as a diminished infiltration of macrophages (F4/80⁺ cells; galectin-3⁺ cells) in the intima of atherosclerotic plaques [23]. Next, a genome-wide miRNA/mRNA microarray profiling showed that miR-125b [together with Taf3, Xpr1 and Ikkβ] was significantly upregulated during ATS progression, effect that was reversed upon CD40 silencing [23]. Further work showed that inclusion of alternative 3’UTRs in the murine Cd34 or in the scavenger receptor SCARB1 transcripts had the potential to alter the pattern of interacting microRNAs, including miR-125 [24,25]. These results prompted us to perform a global study of length variation (including 5’UTR, CDS and 3’UTR domains) in transcripts downregulated during ATS progression and regression in the above ApoE^-/- model of ATS, and here we present the first analysis of the results obtained. These suggest that de-regulation of alternative splicing is an important mechanism to generate length variability in the 5’UTR, CDS and 3’UTR of transcripts in ATS progression. The consequences of this length variability on the regulatory circuitry controlling mRNA stability and function are discussed.

2. Materials and Methods

Mice, Tissue and RNA Extraction

In this work we used female ApoE^-/- mice (homocigous Apoe^tm1Unc in the C57BL/6 background, Jackson Laboratory, Charles River, Wilmington, MA, USA). Mice were treated twice weekly (Monday and Thursday) for 16 weeks with an intraperitoneal administration of 50 μg of an anti-CD40 siRNA (group T, for Treatment) or 50 μg of a scrambled siRNA as control (group C, for Control), and were euthanized the Monday following to the last Thursday injection [23]. Immediately after death, aortic tissue was extracted from mice of the basal group at week 8, for the Control group (ATS progression) at weeks 10 (C10) and 24 (C24), and for the treatment group (anti-CD40 siRNA) also at weeks 10 (T10) and 24 (T24) (see Figure 1). Aortic tissue was then split into two pieces, one of which was used for total RNA extraction as described [23]. All animal studies were in accordance with EU guidelines on animal care and the protocols here used were approved by the ethics committee for animal research of the University of Barcelona-HUB [CEEA, Comité Ètic d’Experimentació Animal], under the procedure name “Immunosuppressió de l’expressió de CD40 mitjançant un siRNA en un model d’arteriosclerosis en ratolins ApoE KO i evaluació dels efectes sobre la malaltía i la seva evolució”, with the approval code: Bell 457/13 and the approval date of 27/06/2013. Details on this experiment have been already published [23].

Microarray Expression Analysis

Total RNA extracted from aortas at the different time-points (10 and 24 weeks, groups C and T) were used for microarray expression analysis. All expression data were normalized to initial levels at week-8 [23].

Sequence Identifiers, Definitions and Nomenclature

In this analysis we used the “seqname” and “GeneSymbol” of the different mRNA species analysed. The “seqname” is analogous to the “Refseq” of the NCBI, and was used as unique identifier for each individual mRNA transcribed isoform, while the “GeneSymbol” identified a gene with all its different variants. Expression was measured by the “FoldChange [FC(Abs)]” or absolute ratio (no log scale) of normalized intensities between two conditions (experimental vs. initial w8), while its statistical significance (as PValue) was calculated from t-test.

Data Analysis

We first identified downregulated transcripts that appeared in the two timepoints (weeks 10 and 24) for each one of the two experimental conditions (groups C and T) by using the sorting function of Excel spreadsheets. Briefly, the two lists of downregulated transcripts (weeks 10 and 24) were sequentially loaded onto the same Excel column and alphabetically sorted acording to their “GeneSymbol”. The mixed list was examined looking for transcript isoforms that shared the same “GeneSymbol” but had different “seqnames”.

Structural Analysis

Length analysis of total transcripts, their CDS, 5’UTR and 3’UTRs was performed at the "ShinyGO v0.76.2” browser [26] at http://bioinformatics.sdstate.edu/go/. This program used Chi-squared and Student's t-tests to identify significant length differences among a list of query transcripts and the the list of exons from all protein-coding genes in the murine genome. Distribution of transcript/subdomain lengths was measured as a “density” function to represent the distribution of quantitative data against a continuous variable, in this case sequence-length in bp.

3. Results

We are interested in the generation of 3’UTR variants in coding transcripts as a mechanism to modulate accession of miRNAs to their target sequences. For this work we used the expression data previously generated on ATS progression/regression in the siCD40/ApoE^-/- model of mice submitted to a high fat diet [23]. In this disease model, 8-weeks old ApoE-deficient mice were treated with a scrambled-sequence siRNA as control (group control, C) or with an anti-CD40 siRNA as anti-ATS treatment (group Treatment, T) and RNA samples were extracted from aortas at weeks 10 or 24 of age (Figure 1). Upon microarray analysis of expression, we detected 489 transcripts downregulated in the C10 group, 1285 in the C24, 300 in the T10 and other 690 in the T24. Here, we will only describe the results obtained with the downregulated transcripts because the group of upregulated mRNAs from the same experiment includes a substantially much higher number of transcripts what makes their analysis more cumbersome. Work is currently in progress to address this point that will be published elsewhere.

In a parallel work, we identified miR-125b as a mediator highly upregulated in ATS progression and whose expression was normalized upon anti-CD40 siRNA treatment [23]. Since a first analysis revealed length variations in the 3’UTRs of miR-125b predicted targets we underwent a more systematic study of the ATS-3’UTRome. Transcripts downregulated at weeks 10 and 24 in the control and treatment groups were submitted to a structural study that was followed by the analysis of the dynamics of length variation in ATS progression (C24 vs. C10) and treatment-related ATS regression (T24 vs. T10/T24 vs. C24). Figure 1 draws the pipeline of the experimental approach.

3’UTR Lengthening in ATS Progression Is Partially Reversed by the Anti-CD40 siRNA Treatment

We studied length variations of the 5’UTR, CDS and 3’UTRs at the ShinyGO server, which compares transcript lengths from a query list (“List”) with those from the reference murine exome from protein-coding genes. As shown in Figure 2, results are represented as superposed plots of “density vs. length in bp”, in which the query is ploted in blue while the control is drawn in red, with T-tests identifying pairs of plots (List vs. control) that are statistically different among them. In our analysis, we used as query the lists of downregulated mRNAs at 10 and 24 weeks from the control (C, ATS progression) and from the treatment (T, ATS regression) groups of the siCD40/ApoE^-/- experiment. Furthermore, we used the same data to characterize changes in the patterns of mRNA isoform expression by taking advantage of the double identification (”GeneSymbol”/”seqname”) of transcripts in the microarray outputs, in which the “seqname” was used as an unique identifier of each individual mRNA isoform while the “GeneSymbol” was the common identifier to all mRNA isoforms transcribed from a single gene. This allowed the identification of expression switches in which one mRNA isoform was replaced by another one from the same gene and facilitated the study of the mechanisms involved.

We first compared total mRNA lengths of the reference murine exome with those of the downregulated transcripts identified at C10 (489 transcripts), C24 (1285 transcripts), T10 (300 transcripts) and T24 (690 transcripts) but no significant length differences were found among the reference and C24, T10 or T24 and only residual differences for C10 (p=0.033; data not shown). We next studied the dynamics of the 3’UTRs in ATS progression, as well as the effects of the anti-CD40 siRNA treatment on this. As can be seen in Figure 2, length distribution of 3’UTRs in the global murine exome followed a bimodal distribution with two peaks at aprox. 100 bp and 1000bp. ATS progression (C10 to C24) was characterized by a 3’UTR lengthening, from a “short” distribution in C10 with a reduced density of longer 3’UTRs in the second peak (Panel C10, Figure 2A, arrow #1) and a general shift of the query curve to the left (Panel C10, Figure 2A, arrow #2), to a “long” distribution in C24 characterised by a clear reduction in the density of 3’UTRs around 100 bp (Panel C24, Figure 2A, arrow #3), an increase in the density of 3’UTRs around 1000 bp (Panel C24, Figure 2A, arrow #1), and a general shift to the right of the query curve (Panel C24, Figure 2A, arrow #2).

On the other hand, the anti-CD40 siRNA treatment partially reversed the above changes detected in 3’UTR length distribution, with Figure 2B showing that 3’UTR length distribution in T10 and T24 almost recovered the normal shape of the control murine exome (p=n.s. among the query and the murine exome).

Length Changes in 5’UTRs and Coding Regions (CDS)

Next, we compared 5’UTR lengths in ATS progression (C24 vs. C10 and reference). Length distribution of 5’UTRs in the reference exome followed a normal curve peaking at 100 bp (red plots at Figure 3A), essentially the same distribution of C10 (p=ns, Figure 3A, blue plot, upper panel). On the contrary, query transcripts at C24 showed a slight but significant shift to longer 5’UTRs as seen by the displacement of sections of the peak to the right (arrows in Figure 3A, bottom panel, blue plot). When lengths of the coding sequences were analysed in ATS progression, a clear, significant, CDS shortening was detected for C10 and especially C24 when compared with the reference exome as seen by the shift of the query plots to the left (Figure 3B, arrows).

Lastly, lengths of 5’UTR and CDS were also examined in the treatment-induced ATS regression samples (T24 vs. C24 and reference). 5’UTRs at T24 were found to be significantly displaced towards longer sequences (p=0.0019, Figure 4A, arrows in bottom panel), suggesting that the siRNA treatment did not have any effect on the normalization of the length of 5’UTRs. A similar result was obtained for the CDS length in T24 (p=0.013, Figure 4B, bottom panel) also evidencing a lack of effect of the treatment on the normalization of CDS length.

Variations in the Number of Exons per Coding Gene during ATS Progression and Treatment-Induced Regression

To deepen on the mechanisms behind length variation during ATS progression and regression we first analysed exon density in the transcripts downregulated at C10, C24 and T24. Figure 5 displays the plots showing the actual number of expressed exons in coding genes (up to 10 exons, red bars), compared with the expected values from the mouse transcriptome (grey bars) and documents variations in the number of exons per coding gene from a distribution similar to that of the control in C10 (p=0.16) to a significant difference (p=0.0021) in C24 when compared with the control, and with an increase in the number of genes with few exons (1, 2 or 3 exons) and a diminution in the number of genes with more exons (8, 9 or 10 exons) when compared with C10. The treatment with the anti-CD40 siRNA partially reversed this effect (T24, p=0.036). Taken together, these results indicated the presence of changes in the exonic composition of transcripts during ATS progression and treatment-induced regression, pointing out to pre-mRNA splicing as the likely mechanism generating transcript variability.

Differential Expression of Transcript Isoforms during ATS Progression or Treatment-Induced Regression

The hypothetical impact of alternative splicing on the generation of mRNA variants was tested by studying isoform expression swifts from C10 to C24, or from T10 to T24. For ATS progression, we loaded onto an Excel spreadsheet the full list of transcripts expressed in C10 (489 entries) and C24 (1285 entries), then we mixed both lists and sorted the resulting list (1774 transcripts) by their “GeneSymbols”. Upon careful examination looking for “GeneSymbol” entries common to the two lists we were able to detect 114 genes simultaneously downregulated in C10 and C24. Although most of these corresponded to a single mRNA isoform (Seqname) expressed in both conditions, we detected 7 cases of isoform switching. In 6 of them a new transcript isoform was expressed “de novo” in C24 (Dusp22, Mta3, Rab11fip5, Rbmx, Skor and Srsf5), while in Ly6e two of the three isoforms expressed in in C10 were repressed in C24 so that only one isoform isoform was expressed in this last. Table 1 summarizes the length variations of the new transcript variants and their functional regions at C24, compared with those expressed at C10. Variations are encoded as longer (L), shorter (S) or equal length (=), and as can be seen, the isoforms differentially expressed in C24 vs. C10 showed length variability mostly in their 5’UTR and CDS, with cryptic or alternative splicing accounting for most of the variability found. Supplementary Table S1 shows the real length values of all these regions, their expression levels and identifies transcript variants.

The same analysis was performed for T10 (300 entries) and T24 (690 entries) and yielded 153 genes simultaneously downregulated in the two conditions, most of them corresponding to isoforms expressed in the two conditions although we also detected 3 cases in which a new isoform was expressed “de novo” in T24 (Ilf3, Pwwp2a and Rbmx, Table 2). Length analysis for these three new isoforms also highlighted variability at their 5’UTR, CDS and 3’UTR, with alternative or cryptic splicing being the most frequent mechanism of action. Supplementary Table S2 shows the real length values of all these regions, their expression levels and identifies transcript variants.

Dysregulation of 3’UTR Splicing Variants in ATS Progression

Lastly, we aimed the mechanisms causing length variations among the transcript variants described above. We first looked for the presence of alternative polyadenylation events [27] but none were detected. On the contrary, a careful analysis highlighted a major occurrence of alternative splicing events, including transcript variants generated by the activation of cryptic or latent splice donor/acceptor sites (Tables 1 and 2). Figure 6A documents the complex splicing of the 3’UTR of Mta3, in which activation of an alternative 5’ splice donor site at the middle of exon 13 and a far downstream acceptor at exon 14 led to the generation of a new transcript variant in C24. This new isoform had a longer 3’UTR of 820 bp vs. 93 bp in the transcript variant 3 expressed in C10. Furthermore, the length of the CDS was also affected by the activation of the new splice site, from the 1542 bp in the variant 3 from C10 to the 1757 for the variant 1 from C24. Figure 6B shows the alternative splice donor/aceptor signals activated in exons 13 and 14, that difer with the canonical donor/acceptor sites and the donor/acceptor signals at the constitutive splicing among exons 12 and 13.

4. Discussion

The sequencing revolution has led to the discovery of new families of regulatory RNAs (lncRNAs, miRNAs, piRNAs, etc.) that establish complex regulatory networks with mRNAs. Alterations in these networks contribute to the disease process and could be also considered as targets of therapeutic intervention [28]. Among these new RNAs, miRNAs have an important role by regulating mRNA stability and function through the targetting of complementary sequences at the 3’UTR of mRNAs. Nevertheless, recent reports are showing that the relationship miRNA/mRNA is much more complex than previously described, since mRNAs avoid miRNA targetting by regulating the length and sequence of their 3’UTRs by alternative polyadenylation or through the regulated use of alternative 3’ terminal exons [3,27]. In this context, we aimed to study the variability in the 3’UTRs of DEGs during ATS progression, or regression after treatment with an anti-CD40 specific siRNA, in aortas from ApoE-deficient mice. Our original microarray expression data identified mRNAs up- or downregulated in our experimental conditions (see Figure 1 and [23]), and here we report the results of the analysis of the downregulated transcripts, because the complexity of the data from the upregulated ones makes it advisable to make a sepparate publication. This work has benefited of the double identification of DEGs in the microarray output, that included the “GeneSymbol”, a gene identifier common to all isoforms, and the “Seqname” an individual isoform identifier.

Our results draw a complex picture of sequence length alterations in ATS progression, since these were not restricted to 3’UTRs but also affected 5’UTRs and the Coding Regions of mRNAs. In this sense, we found a clear lengthening of 3’UTRs during disease progression (C10 to C24, Figure 2A) that was partially reversed by the siRNA treatment (C10 vs. T10 and C24 vs. T24, Figure 2). Variations in 3’UTR length cause the acquisition of new regulatory potentials with regard of stability, translational or subcelular localization by controling accessibility of miRNAs or of RBPs to their target sequences [3,29] and have been mostly linked to Alternative PolyAdenylation (APA) events regulated by general or type-specific factors [30]. In this way, 3’UTR lengthening would increase the regulatory potential by including miRNA or RBP-responsive elements that in tumour suppressor genes have been associated to cáncer resistence by increasing density of mRNA stability regulating RBPs [31], while shortening of 3’UTRs has been linked to tumour progression [32,33] by releasing oncogene transcripts from miRNA control.

On the other hand, we also describe a small but significant lengthening of 5’UTRs from week 10 to week 24 of ATS progression (Figure 3A) not reversed by the antiCD40 treatment (Figure 4A). 5’UTR length variations have been linked to the activation of alternative or upstream start-codons [34], of alternative promoters [35] or to alternative splicing events [36]. Furthermore, 5’UTRs have been shown to be targeted by miRNAs [37,38] in a way depending on their secondary structure [39] (see Figure 7 for a summary of the findings decribed in this manuscript).

We next addressed the mechanisms promoting sequence lengthening by studying transcript isoform switching in ATS progression or regression as a way to describe evolution of gene expression. We were able to detect 9 different genes that were simultaneously expressed in two conditions (6 genes in C10 and C24, other 2 in T10 and T24 and one common to both experiments), which furthermore showed isoform switching with a total of 19 transcript variants expressed in the different conditions (Tables 1 and 2). In this group of transcript isoforms, we identified alternative splicing of 5’ and 3’UTRs as a frequent mechanism contributing to the sequence variability. Furthermore, among these 19 transcripts we found two, Mta3 and Ly6e, expressing alternative 3’ or 5’UTRs (respectively) conformed by the activation of alternative donor/acceptor splice sites. In the case of Mta3 this involved the activation of a cryptic/latent donor splice site [40] at the middle of one exon (Figure 6). This is interesting because alternative splicing has been considered as a minor contributor to the generation of transcript variability since analysis of the superfamily of odorant receptor (OR) genes showed that over 80% of OR mRNAs were submitted to alternative polyadenylation while only a few of these used alternative splicing to generate length variants of 3’UTRs [41].

Regulation of alternative splicing is very complex and involves multiple regulatory sites at the pre-mRNAs (splicing donor and acceptor sites, splicing enhancers and silencers, etc.) recognized by mRNA binding proteins and U-snRNPs ([42] for review). The first step in the maturation of a mRNA is the recognition of the 5’splice site through the functional integration of cis-acting splice signals and splicing regulatory elements (SREs) in the mRNA with the activity of U1-snRNPs and a number of trans-acting splicing factors, all of them working in the context of the secondary structure of the mRNA which regulates accessibility of the splicing machinery to the splice site [43]. Nevertheless, in the genome there are many potential exonic 5’ splicing sites that are not used under physiological conditions, raising the interesting question of the nature of the mechanisms restricting 5’SS selection in normal cells and how are these mechanisms altered in disease to allow the recognition of cryptic or latent 5’SS. This is especially important because improper activation of latent splice sites could result in the incorporation of intronic sequences with potential premature in-frame stop codons to mature mRNAs [44]. In this sense, a recent model of 5’SS selection by Brillen and cols. demonstrated that SREs and SRE-binding proteins were able to block “weak” cryptic 5’splice sites to facilitate recognition of the “strong” actual sites following a sequential, iterative, and position-dependent process [45] while Boehm and cols. proposed the Exon Junction Complexes to suppress cryptic 5’splice sites through the recruitment of the splicing regulator RNPS1 [46]. Furthermore, spliceosomes with noncanonical U1-snRNAs or changes in the stoichiometry of spliceosome components could also contribute to the recognition of variant 5’ss to generate cell/tissue specific patterns of alternative splicing [47,48], and Arafat and Sperling have proposed a quality control mechanism (SoS, suppression of splicing) that would distinguish among normal and latent 5’ splice sites while suppressing these last, through the recognition of a functional ORF (see [49] for a recent review). Lastly, an interesting mechanism has been proposed by Movassat and cols. that linked recognition of the 5’SS in the last exon with the cleavage and polyadenylation site, and identified PA factor CstF64 as a potential regulator of alternative splicing [50]. In this context it is clear that changes on the expression of splicing factors during ATS progression could have an impact on the regulation of the alternative splicing events here described, and work is currently on progress in our group to probe this hypothesis.

Our results should be interpretated with caution and taking into account that we used microarray and not RNA.seq to obtain our expression data [23]. Microarrays represent a quite static view of gene expression, being limited by probe design and by the selection of the arrayed mRNA isoforms. In this sense, we could not discard the presence of additional changes in the patterns of mRNA alternative polyadenylation in our mRNAs, not represented among the microarray probes, yet this apparent limitation has allowed us to highlight isoform switching and alternative splicing as contributors to the generation of transcript length variability in disease progression. In this sense, one interesting possibility is that splicing alterations were pervasive in ATS progression and that by using microarray data we had detected just a reduced part of these alternative/cryptic splicing events. Clarifying this interesting point will require the use of RNA.seq data to identify actual isoform variants. Lastly, splicing is becoming a subject of intense research to develop therapies aimed to normalize splice site usage, anomalous alternative splicing events and to correct splicing factor irregularities (see [51] for a recent review). Characterizing the splicing alterations in ATS will surely provide new targets of pharmacological intervention.

In conclussion, here we have characterized the length variability of transcript variants dowregulated during ATS progression. We have uncovererd general length variations, not only at the 3’UTRs of these transcripts as expected, but also at their 5’UTR and CDS, and have given evidences for the involvement of alternative splicing as a frequent mechanism for the generation of length variability. To the best of our knowledge, this is the first report linking ATS development with changes in the length of ATS-DEGs, and highlights splicing factors as possible targets of therapeutical intervention in ATS progression.