Predominant contribution of cis-regulatory divergence in the evolution of mouse alternative splicing

Abstract

Divergence of alternative splicing represents one of the major driving forces to shape phenotypic diversity during evolution. However, the extent to which these divergences could be explained by the evolving cis-regulatory versus trans-acting factors remains unresolved. To globally investigate the relative contributions of the two factors for the first time in mammals, we measured splicing difference between C57BL/6J and SPRET/EiJ mouse strains and allele-specific splicing pattern in their F1 hybrid. Out of 11,818 alternative splicing events expressed in the cultured fibroblast cells, we identified 796 with significant difference between the parental strains. After integrating allele-specific data from F1 hybrid, we demonstrated that these events could be predominately attributed to cis-regulatory variants, including those residing at and beyond canonical splicing sites. Contrary to previous observations in Drosophila, such predominant contribution was consistently observed across different types of alternative splicing. Further analysis of liver tissues from the same mouse strains and reanalysis of published datasets on other strains showed similar trends, implying in general the predominant contribution of cis-regulatory changes in the evolution of mouse alternative splicing.

Keywords: alternative splicing; cis‐regulation; evolution.

Figure 1

Study design Fibroblast cells were isolated from adult C57BL/6J, SPRET/EiJ, and the F1 hybrid mice and cultured. PolyA RNAs prepared from each cell line were sequenced on an Illumina HiSeq 2000/2500 platform.

Figure EV1

FDR estimation for each |ΔPSI| cutoff FDR for parental (triangle) and allelic (circle) splicing comparison (y-axis) was plotted against different |ΔPSI| cutoffs (x-axis). For each value of x from 0.01 to 0.20 increasing by 0.01, we performed independent 100 bootstrapped label permutations of replicate 2 and replicate 3, respectively. For each of the 100 shuffled sets, we calculated the number of events passing the threshold (false positives), that is BF > 5 in all the replicates and average |ΔPSI| > x. Then, for each of the 100 permutations of each value x, the FDR was estimated as false positives divided by the number of real events passing the threshold, including both false positives and true positives.

Figure EV2

Illustration of data filtering based on mock F1 hybrid

1. A MA plot comparing the PSI values in parental strains and their downsampling datasets. The local standard deviation for each comparison was also indicated (see Materials and Methods).

2. B MA plot comparing the PSI values in parental strains and those estimated based on mock F1 dataset. The red dots represented the outliers with inconsistent PSI values between parental strain and mock F1 dataset.

3. C–E After filtering, the PSI values for C57BL/6J (C), SPRET/EiJ (D), and their difference (E) correlated well between parental strains and mock F1 hybrid (R² = 0.99, 0.99 and 0.94, respectively).

Figure EV3

Illustration of PacBio sequencing of splicing event-spanning cDNA products For each candidate event, RT–PCR primers were designed in the conserved regions of the constitutive exons to amplify both isoforms from the two alleles/strains. The PCR products were then sequenced at full length using PacBio RS system.

Figure 2

Dissection of cis- and trans-regulatory contributions in alternative splicing

1. Scatterplot comparing parental splicing differences (dots, denoted as F0 hereafter) or allelic splicing differences (triangles) estimated based on Illumina RNA-seq results (y-axis) to those based on PacBio sequencing of splicing event-spanning cDNA products (x-axis) (R² = 0.91 and 0.92 for comparison of parental and allelic difference, respectively).

2. Scatterplot comparing splicing difference in parental strains (y-axis) versus the allelic difference in F1 hybrid (x-axis). After filtering using mock F1 hybrid, 5,802 AS events were expressed in F1 hybrid (gray dots). Among these, 417 AS events were divergent between parental strains (black dots), of which 255 (indicated as “+”) and 62 (indicated as “×”) exhibited significant cis- and trans-regulatory divergence, respectively.

3. Examples of cis (upper panel)- and trans (lower panel)-regulatory divergence in alternative splicing. The RNA-seq read densities supporting the inclusion and exclusion of exons were shown in the left plot. The estimated PSI values and 95% confidence intervals were shown in the right plot.

4. Percentage of cis- and trans-divergent events for the five AS types separately (numbers of events for each type were indicated above bars).

Figure EV4

Dissection of cis- and trans-regulatory contributions in alternative splicing at different |ΔPSI| cutoffs

1. A–C Scatterplot comparing splicing differences in parental strains (y-axis) versus the allelic differences in F1 hybrid (x-axis) at different |ΔPSI| cutoffs [|ΔPSI| > 0 (A), 0.05 (B) and 0.15 (C)]. After filtering using mock F1 hybrid, 5,802 AS events were expressed in F1 hybrid (gray dots). Among these, 615 (A)/530 (B)/336 (C) AS events were divergent between parental strains (black dots), of which 376 (A)/320 (B)/209 (C) (indicated as “+”) and 115 (A)/86 (B)/43 (C) (indicated as “×”) exhibited significant cis- and trans-regulatory divergence, respectively.

2. D–F Percentage of cis- and trans-divergent events for the five AS types separately at different |ΔPSI| cutoffs [|ΔPSI| > 0 (D), 0.05 (E), and 0.15 (F)].

3. G Contributions of cis (indicated as triangle)-/trans (indicated as circle)-regulatory divergence (y-axis) to parental divergent AS events with different effect sizes (|ΔPSI|, x-axis). A total of 417 divergent events between parental strains (see Fig2B) were grouped into 7 categories according to the |ΔPSI| values: (0.1, 0.2), (0.2, 0.3), (0.3, 0.4), (0.4, 0.5), (0.5, 0.6), (0.6, 0.7), and (0.7, 1.0). The number of events in each category was marked. While cis-regulatory divergence always played the predominant role in determining parental AS divergence with different effect sizes, its relative contribution slightly decreased with the decreasing effect size.

Figure EV5

Divergent AS events identified using Fisher’s exact test

1. A, B Venn diagram showing the overlap of the divergent events identified by Fisher’s exact test and MISO in parental strains (A) and in F1 hybrid (B).

2. C Scatterplot comparing splicing difference in parental strains versus the allelic difference in F1 hybrid identified by Fisher’s exact test. After filtering using mock F1 hybrid, 5,802 AS events were expressed in F1 hybrid (gray dots). Among these, 626 AS events were divergent between parental strains (black dots), of which 357 (indicated as “+”) and 72 (indicated as “×”) exhibited significant cis- and trans-regulatory divergence, respectively.

3. D Percentage of cis- and trans-divergent events for the five AS types separately using Fisher’s exact test.

Figure EV6

Dissection of cis- and trans-regulation in alternative splicing between C57BL/6J and SPRET/EiJ liver samples

1. Scatterplot comparing splicing difference between C57BL/6J and SPRET/EiJ liver samples versus their allelic difference in F1 hybrid liver sample. After filtering using mock F1 hybrid, 4,124 AS events were expressed in F1 hybrid (gray dots). Among these, 336 AS events were divergent between parental strains (black dots), of which 196 (indicated as “+”) and 38 (indicated as “×”) exhibited significant cis- and trans-regulatory divergence, respectively.

2. Percentage of cis- and trans-divergent events for the five AS types separately.

Figure EV7

Dissection of cis- and trans-regulation in alternative splicing between C57BL/6J and CAST/EiJ

1. Scatterplot comparing splicing difference between C57BL/6J and CAST/EiJ versus their allelic difference in F1 hybrid. After filtering using mock F1 hybrid, 2,042 AS events were expressed in F1 hybrid (gray dots). Among these, 79 AS events were divergent between parental strains (black dots), of which 44 (indicated as “+”) and 6 (indicated as “×”) exhibited significant cis- and trans-regulatory divergence, respectively.

2. Percentage of cis- and trans-divergent events for the five AS types separately.

Figure EV8

Illustration of the regions flanking the AS events For SE, the alternative exons and their flanking 100 nt intron sequences were considered; for RI, the retained introns and their flanking 100 nt exon sequences were considered. For A3SS or A5SS, the alternative exon regions and their flanking 100 nt exon/intron sequences were considered. For MXE, both alternative exons and their flanking 100 nt intron sequences were considered.

Figure 3

Genomic features that correlate with cis-regulatory alternative splicing divergence

1. The cumulative distribution function (CDF) of frequencies of nucleotide variants in the AS flanking regions for the events with cis-regulatory divergence (black) and controls (grey). Compared with controls, the events with significant cis-regulatory impact had higher sequence divergence in the flanking regions. The P-values were calculated by the Mann–Whitney U-test.

2. 36.2 and 11.5% of the events with significant cis-regulatory divergence (black) and control events (gray) had sequence divergence at their exact splice sites, respectively (***P = 9.21e-14, Fisher’s exact test).

3. CDF of allelic differences in splicing site strengths due to sequence variants at the exact splicing sites plotted for cis-regulatory divergent events (black) and control events (grey), separately. The splicing site strengths changed more in the events with cis-regulatory events than in those without. The P-values were calculated by the Mann–Whitney U-test.

4. An example showing that a SNV at the canonical GU/AG sites (indicated as an arrow) resulted in complete functional abortion of the corresponding splice sites. The substitution of the AG to GG in SPRET/EiJ disrupted the splicing site and thereby facilitated the use of a downstream splicing acceptor.

Figure EV9

Genomic features that correlate with cis-regulatory alternative splicing divergence for each AS type separately

1. CDF of frequencies of nucleotide variants in the AS flanking regions for the events with cis-regulatory divergence (black) and controls (grey) for A3SS, A5SS, MXE, RI, and SE, respectively.

2. Percentages of the events with significant cis-regulatory divergence (black) and controls events (gray) that had sequence divergence at the exact splice sites for A3SS, A5SS, MXE, RI, and SE, respectively.

Figure EV10

Genomic features that correlate with cis-regulatory alternative splicing divergence identified in the liver sample

1. CDF of frequencies of nucleotide variants in the AS flanking regions for the events with cis-regulatory divergence (black) and controls (grey) identified in liver sample. Compared with controls, the events with significant cis-regulatory impact also had higher sequence divergence in the flanking regions.

2. In liver sample, 37.0 and 10.8% of the events with significant cis-regulatory divergence (black) and controls (gray) had sequence divergence at the exact splice sites, respectively.

3. CDF of allelic differences in splicing site strengths due to sequence variants at the exact splicing sites plotted for cis-regulatory divergent events (black) and controls (grey) identified in liver sample. The splicing site strengths changed more in the events with cis-regulatory events than in those without.

Figure EV11

Sashimi plot for the splicing patterns of the SE event in Trim26 gene from fibroblast cell line as well as brain tissues of five mouse strains The top four rows represented splicing patterns for C57BL/6J and SPRET/EiJ strains and their alleles in F1 hybrid. The bottom five rows represented splicing patterns for brains tissues of the five mouse strains. PWK/PhJ and CAST/EiJ had a similar splicing pattern as C57BL/6J, but different from SPRET/EiJ. Four variants located in the flanking regions, two of which correlated with the species-specific splicing pattern and were analyzed using minigene assays (see Fig4).

Figure 4

Minigene analysis for the cis-divergent SE event in Trim26 gene

1. Schematic diagrams of minigene constructs for validating the cis-divergent SE event identified in Trim26 gene. Two candidate variants, one SNV and one insertion (INS), were indicated. Four constructs were prepared in C57BL/6J background with no variant, only insertion, only SNV, and both insertion and SNV, respectively (see Materials and Methods).

2. Minigene assays of the four constructs transfected into HEK293T cells suggested only the insertion contributed to this divergent SE event. The gel image illustrated RT–PCR products from these constructs. The barplot below the gel image represented the PSI values calculated from triplicates of RT–PCR products using Agilent Bioanalyzer 2000 system (see Materials and Methods, for minigene assays in NIH3T3 cells, see FigEV12).

Figure EV12

Minigene analysis for the cis-divergent SE event in Trim26 gene in NIH3T3 cells Label is the same as in Fig4B.