Mammalian splicing divergence is shaped by drift, buffering in trans, and a scaling law

The lower triangle shows the scatter plots of percent spliced in of sample pairs, and the upper triangle shows Pearson’s correlation coefficients of sample pairs. Correlation coefficients between replicates in each tissue/cell line are highlighted in green. The lower correlation between two replicate ES samples is likely due to the fact that the two samples are from two distinct ES clones.

(A) Statistics of events expressed in different numbers of tissues. Numbers on the top of bars indicate the total number of events in corresponding categories. (B, C, D, E) Among the 11,729 events expressed in two or more tissues, 4,328 were classified as “Non-Differential” (switch score < 0.1), 3,229 as “Low” (0.1 ≤ switch score < 0.2), 1,443 as “Moderate-Low” (0.2 ≤ switch score < 0.3), 1,331 as “Moderate-High” (0.3 ≤ switch score < 0.5) and 1,398 as “Switch-Like” (switch score > 0.5). According to the error hypothesis, we would expect that genes with higher switch scores have lower expression levels and higher dN/dS ratios than those with smaller differences. However, we found patterns resembling U-shaped or hump-shaped distributions: “Non-Differential” and “Switch-Like” events both affect coding regions more often than events with “Moderate-Low” or “Moderate-High” switch scores (B; numbers in each bar denote the number of coding and non-coding SE events in the corresponding group). (C, D) Furthermore, the genes harboring events in either of these categories have higher expression levels (C) and those with “Switch-Like” events also have significantly lower dN/dS ratios than those with events in the “Low” and “Moderate-Low” categories (D). Importantly, the flanking regions of “Switch-Like” events also are more conserved (with higher PhastCons scores) than those of events in all three intermediate categories. (E) Events expressed in two or more tissues and with exon size ≥60 bp and intron size ≥200 bp were used. The mean PhastCons score of each position across events within the same tissue-regulatory group for exonic regions (including 30 bp at both 5′ and 3′ end of alternative exon) and intronic regions (200 bp flanking sequences on each side) around alternative splicing sites was compared between Switch-Like events and the other four tissue-regulatory groups. (F) Enrichment of micro-exons in Switch-Like events. Events expressed in two or more tissues are considered. Percentages of micro-exons (exons with length ≤ 30 bp) and regular exons are indicated for the five tissue-regulatory groups. Numbers within and on top of bars denote the number of regular exons and micro-exons in each group, respectively. (B, C) One-sided Wilcoxon rank-sum test was used in (B) and (C). Box plot elements: boxes span values from the first (Q1) to the third (Q3) quartile, with a horizontal line indicating the median (Q2). The whiskers extend from the edges of the box to the points representing the largest and smallest observed values within the 1.5 * IQR (interquartile range = Q3 − Q1) from the edges of the box. Outlier values are plotted as dots outside the whiskers. n.s. not significant; *P < 0.05; **P < 0.01; ***P < 0.001.

(A, B, C) Median, maximum, and minimum of absolute ΔPSI across divergent tissues are compared between All-Divergent (steel blue) and Some-Divergent (golden yellow) events (one-sided Wilcoxon rank-sum test). (D, E, F) Comparing median, maximum, and minimum of absolute ΔPSI between events in All-Divergent with (deep sky blue) and without (light sky blue) splice site change (one-sided Wilcoxon rank-sum test). (G) dN/dS ratios of genes associated with SE events belonging to different divergence groups (indicated by different colors). Events expressed in two to seven tissues whose associated genes have orthologs in rat are used (one-sided Wilcoxon rank-sum test). (H) Genes associated with “All-Divergent” events have lower expression levels than those associated with “Non-Divergent” events (one-sided Wilcoxon rank-sum test). Events expressed in two to seven tissues are used. (I) Comparison of dN/dS ratio between genes associated with divergent events (cyan) and those associated with non-divergent events (red) in each individual tissue (one-sided Wilcoxon rank-sum test). (J) Comparison of gene expression level between genes associated with divergent and those associated with non-divergent events in each tissue (one-sided Wilcoxon rank-sum test). Box plot elements: center line, median; box limits, lower and upper quartiles; whiskers, lowest and highest value within 1.5 IQR. n.s., not significant; *P < 0.05; **P < 0.01; ***P < 0.001.

(A, B, C, D, E, F, G) We down-sampled junction reads of each event in either allele across all samples to the same level of 20 reads, and determined divergent and non-divergent splicing based on percent spliced in (PSI) difference (divergent: ΔPSI ≥ 0.1). The x-axis and y-axis denote PSI difference before down-sampling and after down-sampling, respectively. The colors distinguish “Divergent” (red) and “Non-Divergent” (light blue) in the original dataset. The proportions of consistent classification of “Divergent” events and “Non-Divergent” events between down-sampled (use PSI difference of 0.1 as cutoff) and original datasets have been labeled near the corresponding points. (H) The mean and 95% confidence interval (calculated by down-sampling for 100 times) of proportions of consistent classification of “Divergent” events (red) and “Non-Divergent” events (blue) between down-sampled and original datasets. (I) The mean and 95% confidence interval (calculated by down-sampling for 100 times) of Pearson’s correlation coefficient (R) of PSI divergence between down-sampled and original datasets.

(A) Comparison of average gene expression level among the three divergence groups. Numbers in parentheses indicate the number of genes within corresponding groups. (One-sided Wilcoxon rank-sum test). (B) Comparison of average gene expression level of genes with “Some-Divergent” events in divergent tissues versus non-divergent tissues, the number in parentheses is the number of “Some-Divergent” events (One-sided Wilcoxon signed-rank test). (C) Distribution of Spearman’s correlation coefficients between |ΔPSI| and gene expression across tissues for events expressed in ≥3 tissues and divergent in ≥1 tissue(s) (one-sided Wilcoxon signed-rank test). Box plot elements: center line, median; box limits, lower and upper quartiles; whiskers, lowest and highest value within 1.5 IQR. n.s., not significant; *P < 0.05; **P < 0.01; ***P < 0.001.

The bar plot shows the proportions of the three divergence groups in each of the five tissue-regulatory groups. Allelic events expressed in two or more tissues which are also included in the set of expressed events based on total reads were used. The classification of tissue-regulatory groups is based on the switch score calculated from total reads. The percentages of “Non-Divergent” and “Some-Divergent” events are indicated in each bar, and the total number of events in each group is indicated in parentheses. When comparing the alternative splicing patterns across tissues, “Switch-Like” events also are more conserved between the C57BL/6J and the SPRET/EIJ allele (see also Table S8), indicating that for the majority of Switch-Like events, the tissue-switch between the two isoforms is functional in both alleles.

(A) Allele-specific splicing patterns of exon 5 in Pkp4 across tissues. Two single-nucleotide variants between the two species have been identified in the 5′ splicing site of this exon. (B) Tissue-dependent allele-specific splicing patterns of exon 7 in Apbb3. Three variants including two small indels (indicated with asterisks) and one single-nucleotide variant (ellipse) were predicted to contribute. (C) Allele-specific splicing pattern of exon 6 in Thyn1. This event is divergent between the two alleles in all expressing tissues. (D) Allele-specific splicing pattern of exon 19 in Limch1. This exon is partially included only in heart and cortex in the SPRET/EiJ allele, whereas it is completely skipped in all expressing tissues in the C57BL/6J allele. (E) Allele-specific splicing pattern of exon 7 in Mff. This exon is heart-specifically included in the SPRET/EiJ allele, while it is completely skipped in all expressing tissues in the C57BL/6J allele. The gene local structure around the alternative exon and sequence conservation are presented with the alternative exon highlighted in red.

The left panel shows the relationship between total reads count (all junction reads, log₂ transformed) and PSI difference between replicates based on the original dataset. The right panel shows the relationship between original read count (log₂ transformed) and PSI difference estimated based on the down-sampled dataset. The red line in each scatter plot was fitted by linear regression. Pearson’s correlation coefficient (R) and P-value are presented for each subfigure.

We changed the threshold of 0.1 as in Fig 5B to 0.2, and the number of events in “Scenario 3” (761) is still much higher than that in “Scenario 1” (1) and “Scenario 2” (121).

(A) Tissue-regulatory patterns of splicing divergence in four selected tissues. The x-axis presents the ΔPSI of the C57BL/6J (BL6) allele between the two tissues with maximal difference in allelic ΔPSI, and the y-axis presents the corresponding information for the SPRET/EiJ (SPR) allele. (B) Definition of the conserved allele by using two out-groups (Mus caroli and Mus pahari). (C) Micro-exon events in Scenario 2 have larger PSI differences between tissues than those falling in Scenario 3. The magnitude of the difference is even bigger than the difference between the two scenarios in all exons. The y-axis represents the maximum of {∆PSI between tissues of the BL6 allele, ∆PSI between tissues of the SPR allele}. The median value of the ∆PSI in each scenario is indicated in the box, and numbers in parentheses represent the numbers of events belonging to the corresponding scenarios. One-sided Wilcoxon rank-sum test was used to test for significance. Box plot elements: center line, median; box limits, lower and upper quartiles; whiskers, lowest and highest value within 1.5 IQR. n.s., not significant; *P < 0.05; **P < 0.01; ***P < 0.001.

We used the lower percent spliced in (PSI) value of the two alleles as starting PSI and the higher one as the final PSI to calculate the effects on splicing efficiency of cis-regulatory variants under treatment and control samples separately and compared “A” between the two conditions. Wilcoxon rank-sum test was used to test the difference.