Alternative splicing of RNA triplets is often regulated and accelerates proteome evolution

Abstract

Thousands of human genes contain introns ending in NAGNAG (N any nucleotide), where both NAGs can function as 3' splice sites, yielding isoforms that differ by inclusion/exclusion of three bases. However, few models exist for how such splicing might be regulated, and some studies have concluded that NAGNAG splicing is purely stochastic and nonfunctional. Here, we used deep RNA-Seq data from 16 human and eight mouse tissues to analyze the regulation and evolution of NAGNAG splicing. Using both biological and technical replicates to estimate false discovery rates, we estimate that at least 25% of alternatively spliced NAGNAGs undergo tissue-specific regulation in mammals, and alternative splicing of strongly tissue-specific NAGNAGs was 10 times as likely to be conserved between species as was splicing of non-tissue-specific events, implying selective maintenance. Preferential use of the distal NAG was associated with distinct sequence features, including a more distal location of the branch point and presence of a pyrimidine immediately before the first NAG, and alteration of these features in a splicing reporter shifted splicing away from the distal site. Strikingly, alignments of orthologous exons revealed a ∼15-fold increase in the frequency of three base pair gaps at 3' splice sites relative to nearby exon positions in both mammals and in Drosophila. Alternative splicing of NAGNAGs in human was associated with dramatically increased frequency of exon length changes at orthologous exon boundaries in rodents, and a model involving point mutations that create, destroy, or alter NAGNAGs can explain both the increased frequency and biased codon composition of gained/lost sequence observed at the beginnings of exons. This study shows that NAGNAG alternative splicing generates widespread differences between the proteomes of mammalian tissues, and suggests that the evolutionary trajectories of mammalian proteins are strongly biased by the locations and phases of the introns that interrupt coding sequences.

Figure 1. Alternative splicing of tissue-specific NAGNAGs is highly conserved.

(A) Short reads were aligned to the intron-proximal and intron-distal splice junctions of NAGNAG splicing events in order to estimate isoform ratios. (B) Estimated proximal isoform usage (ψ) for a NAGNAG which inserts/deletes a predicted phosphorylation site in far upstream element binding protein 1 (FUBP1). Phosphorylation site and corresponding kinase were predicted by Scansite (Scansite z-score −3.024) . Error bars indicate the 95% binomial confidence interval. (C) Number of reading frame-preserving alternative splicing events in protein-coding regions, with both isoforms expressed at ≥5% in at least one tissue (see also Table S1). (D) A NAGNAG which inserts/deletes an arginine in RNA recognition motif 4 (RRM4) of the splicing factor PTBP2 is deeply conserved. Alignment of orthologous 3′ splice site sequences shown below the NMR structure (PDB accession 2ADC, displayed with PyMOL) of the highly homologous PTBP1 protein (green) complexed with RNA (red) . Boxed is K489 of PTBP1, which is homologous to the arginine shown in PTBP2, and hydrogen bonds to the RNA backbone (dotted yellow line). Putative branch point based on location of the first upstream AG, the sequence motif identified in , and the pattern of sequence conservation. (E) Conservation of alternative splicing between orthologous human and mouse NAGNAGs increases with tissue specificity. NAGNAGs that were alternatively spliced in human (left) and mouse (right) were grouped by switch score—defined as the maximum ψ difference between tissues—as indicated by colors, and the fraction of orthologs which were alternatively spliced in the other species is shown. Error bars indicate 95% binomial confidence intervals.

Figure 2. Increased sequence conservation upstream of tissue- and developmentally-regulated NAGNAGs.

(A) NAGNAG ψ estimates are highly consistent in brain RNA-Seq data from the mouse strains DBA/2J and C57BL/6J. Only NAGNAGs with both isoforms expressed at ≥5% in either strain are shown. The 75^th percentile of the deviation from the line y = x is shown in gray. (B) NAGNAG ψ estimates are quantitatively conserved between human and mouse brain. Only NAGNAGs with both isoforms expressed at ≥5% in either species and satisfying |proximal 3′ splice site score – distal 3′ splice site score|≤0.5 bits are plotted (splice sites scored by MaxEnt model [36]). Deviation from y = x shown as in (A). (C) Sequence conservation of human NAGNAGs, where all NAGNAGs are aligned by their 3′ splice site junctions and grouped by switch score. Mean (solid line) and standard error of the mean (shaded area about solid line) of phastCons score shown by position (averaged over a 2 nt sliding window) for each switch score category. Analysis restricted to human NAGNAGs for which the two AGs were conserved at the sequence level in mouse. (D) As in (C), but grouped by switch score in mouse and restricted to mouse NAGNAGs for which the two AGs were conserved at the sequence level in human. (E) As in (D), but for NAGNAGs in Drosophila melanogaster, with switch score defined across developmental stages rather than between tissues. Analysis restricted to D. melanogaster NAGNAGs for which the two AGs were conserved at the sequence level in D. yakuba.

Figure 3. Variation in 3′ splice site features are associated with differences in NAGNAG splicing.

(A) A simple biophysical model of NAGNAG splicing accurately models mean isoform usage across tissues as a function of difference in 3′ splice site score. Each point represents a single human NAGNAG, and the solid and dashed black lines show the mean ψ (across values for individual NAGNAGs with similar splice site score difference, with sliding window of 3.25 bits) and the standard deviation about the mean. The solid red line shows the prediction based on the model for parameters Q = 0.55 and B = 0.58, and the dashed red line indicates the standard deviation about the model mean expected from measurement error. The horizontal and vertical dashed lines indicate the splice site score difference (approximately 1 bit) at ψ = 50%. (B) The −3 bases largely determine whether a NAGNAG is alternatively spliced. We grouped NAGNAGs in the human genome according to their −3 bases and computed the fraction of each group which expressed the proximal (black) or distal (blue) isoform at ≥5% in at least one tissue. (C) Constitutive 3′ splice sites (top, YAG), YAGYAGs which express the proximal isoform at ≥75% in all tissues (middle, YAGYAG proximal major), YAGYAGs which express the distal isoform at ≥75% in all tissues (middle, YAGYAG distal major), and strongly regulated YAGYAGs (bottom, YAGYAG strongly regulated) all exhibit distinct upstream sequence preferences. The x-axis shows the position relative to the 3′ splice site (YAG) or proximal 3′ splice site (YAGYAG), and arrows indicate the 3′ splice site that is predominantly used. Figure was created with WebLogo . Human and mouse YAGYAGs were grouped together to increase the statistical signal for (C–F). (D) Distal major YAGYAGs have shorter polypyrimidine tracts (p<0.001 relative to proximal major class, Kolmogorov-Smirnov test). Plot shows median length of the polypyrimidine tract, estimated as the first stretch of ≥5 consecutive pyrimidines upstream of the −3 position. Error bars indicate the standard deviation of the median, estimated by bootstrapping (the error bars for “CJ” were too small to be visible). (E) Distal major YAGYAGs have higher CT and TC dinucleotide content (p<0.005 relative to proximal major class, Kolmogorov-Smirnov test). Median CT and TC dinucleotide content of the polypyrimidine tract, computed as the fraction of the polypyrimidine tract composed of CT dinucleotides, with an optional T at the beginning or C at the end. Error bars indicate the standard deviation of the median, estimated by bootstrapping. (F) The AG exclusion zone is more distally located in distal major YAGYAGs (p<0.001 relative to proximal major class, Kolmogorov-Smirnov test). Position of the first AG dinucleotide upstream of the −15 position is shown. Thick bars indicate the median positions, and boxes extend from the first to third quartiles.

Figure 4. Specific intronic sequence features regulate NAGNAG splicing.

(A) Illustration of NAGNAG minigene constructs, designed to test the roles of the branch point to 3′ splice site distance and of the −4 base in NAGNAG splicing. A short segment of intronic sequence spanning the branch point to the 3′ splice site of the PTBP2 NAGNAG was cloned upstream of the IGF2BP1 exon. To confirm the importance of a pyrimidine at the −4 position for distal NAG use, the effects of all four nucleotides at the −4 position were tested. The branch point to 3′ splice site distance was varied by introducing nucleotides (underlined in orange) in constructs containing the PTBP2 branch point sequence, or by removing nucleotides (indicated by green dots) in constructs containing the IGF2BP1 branch point sequence. Locations of RT-PCR primers are indicated by arrows. (B) Proximal isoform expression increased dramatically after the introduction of a purine at the −4 position. Splicing was monitored after minigene transfection into HEK293T cells by RT-PCR. Mean and standard deviation of at least three independent transfections are shown. A representative gel is shown below (top and bottom bands represent proximal and distal isoforms, respectively). (C) As in (B), but varying the branch point to 3′ splice site distance in the context of the native nPTB branch point sequence. The distance was increased by insertion of four or seven nucleotides of sequence of varying purine/pyrimidine composition as shown in (A). (D) As in (C), but decreasing the branch point to 3′ splice site distance in the context of the exogenous IGF2BP1 BPS by deletion of three or six bases as shown in (A).

Figure 5. NAGNAGs are associated with accelerated protein evolution at exon-exon boundaries.

(A) Alignment of portions of exons 10 and 11 of TRIM28 gene from three mammals, illustrating a shift in the upstream boundary of exon 11 between human and rodents. Exonic sequence shown in capitals; intronic sequence in lower case. (B) Gain/loss of exonic sequence between human and mouse occurs preferentially at 3′ splice sites (p<10⁻⁶, permutation test). The fraction of aligned orthologous human and mouse exons with gaps at each position is shown; the background level (mean fraction across the indicated region excluding the 3′ splice site) is shown by the dotted yellow line; the right-hand axis shows enrichment relative to this background. (C) As in (B), but restricted to gaps of length three bases. Preferential occurrence at 3′ splice sites was highly significant (p<10⁻⁶, permutation test). (D) Similar to (C), but based on alignments of orthologous D. melanogaster and D. yakuba exons. Preferential occurrence at 3′ splice sites was highly significant (p<10⁻⁶, permutation test). (E) Similar to (C), but based on alignments of orthologous C. elegans and C. briggsae exons. Preferential occurrence at 3′ splice sites was highly significant (p<10⁻⁶, permutation test). (F) Residual NAG motif at exons whose boundaries changed in the rodent lineage. Orthologous mouse and rat exons were classified as unchanged (top), expanded by three bases (middle), or contracted by three bases(bottom) based on comparison to an outgroup (human, cow, chicken, or Xenopus laevis), aligned to the inferred location of the ancestral 3′ splice site (dotted line). Information content of each position is shown relative to a uniform background composition. (G) Exons whose 3′ splice site boundaries differ by three bases between rat and mouse are 7.5 times as likely to have a NAGNAG in the human ortholog as exons whose boundaries did not change (p-value for difference<10⁻²⁴ by Fisher's exact test). Error bars indicate the 95% binomial confidence interval. (H) Rodent exons orthologous to alternatively spliced human NAGNAG exons (left) are much more likely to exhibit exon boundary changes of three base pairs than those orthologous to constitutively spliced human NAGNAGs (right) (p-value for difference<10⁻¹⁰ by Fisher's exact test). Blue and gray bars in (H) represent subsets of blue and gray bars in (G), respectively. Error bars indicate the 95% binomial confidence interval. (I) Frequency of encoded amino acids that occur opposite gaps at the 3′ splice site in alignments of human and mouse exons is plotted above, overall (pink) and separately by the phase of the upstream intron (i.e., the number of bases, if any, in the last incomplete codon of the upstream exon); amino acid frequency at background positions (4 codons downstream of the 3′ splice site) is shown below. The Shannon entropy (a measure of randomness) of each amino acid frequency distribution is also shown.