Alternative splicing and gene expression play contrasting roles in the parallel phenotypic evolution of a salmonid fish

Abstract

Understanding the contribution of different molecular processes to evolution and development is crucial for identifying the mechanisms of adaptation. Here, we used RNA-sequencing data to test the importance of alternative splicing and differential gene expression in a case of parallel adaptive evolution, the replicated postglacial divergence of the salmonid fish Arctic charr (Salvelinus alpinus) into sympatric benthic and pelagic ecotypes across multiple independent lakes. We found that genes differentially spliced between ecotypes were mostly not differentially expressed (<6% overlap) and were involved in different biological processes. Differentially spliced genes were primarily enriched for muscle development and functioning, while differentially expressed genes were involved in metabolism, immunity and growth. Furthermore, alternative splicing and gene expression were mostly controlled by independent cis-regulatory quantitative trait loci (<3.4% overlap). Cis-regulatory regions were associated with the parallel divergence in splicing (16.5% of intron clusters) and expression (6.7%-10.1% of differentially expressed genes), indicating shared regulatory variation across ecotype pairs. Contrary to theoretical expectation, we found that differentially spliced genes tended to be highly central in regulatory networks ("hub genes") and were annotated to significantly more gene ontology terms compared to nondifferentially spliced genes, consistent with a higher level of pleiotropy. Together, our results suggest that the concerted regulation of alternative splicing and differential gene expression through different regulatory regions leads to the divergence of complementary processes important for local adaptation. This provides novel insights into the importance of contrasting but putatively complementary molecular processes in rapid parallel adaptive evolution.

Keywords: adaptive divergence; alternative splicing; convergent evolution; ecological speciation; gene expression; parallel evolution; salmonid.

FIGURE 1

Sampling and differential gene expression. (a) Map showing the locations of the focal populations (Awe, Tay and Dughaill [Dug]). Representative pictures of the benthic and pelagic ecotypes per lake are shown. Coordinates and additional sample information are given in Table 1. (b) PCA based on expression data (GEx‐PCA) of all genes (n = 19,623) for all individuals (n = 24). Round large points show the respective centroids for each ecotype, and sympatric ecotypes are connected by vectors. Small unicolour points represent individuals. Populations are coded by shape and ecotypes by colour. (c) Venn diagram displaying the amount of overlap of differentially expressed genes between sympatric ecotypes (FDR < 0.05) across lakes. (d) Distribution of z‐scores derived from the RDA, depicting association of gene expression with ecotype across the genome. The red dashed lines highlight the significance threshold (|z| >2). Genes with negative z‐scores were associated with the benthic ecotype across lakes, and genes with positive thresholds were associated with pelagic ecotypes. Illustrations highlight differences in head shape between benthic and pelagic ecotypes. Core genes that were differentially expressed in at least two ecotype pairs and significantly associated with ecotype across lakes are highlighted in red, with gene IDs from the Arctic charr genome annotation. Chromosomes are highlighted by alternative colours, and unplaced scaffolds are located at the end

FIGURE 2

Shared patterns of differential splicing across lakes. (a) PCA based on intron excision ratios (n = 18,207 introns) for all individuals. See Figure 1 for explanation of the key. (b) Venn diagrams showing the amount of overlap of differentially spliced genes based on differential exon usage and differential intron excision between sympatric ecotypes across lakes. Two genes were detected using both approaches. (c, d) Gene models illustrating alternative splicing patterns for (c) U2AF2 and (d) SH3BGRL based on exon usage. The expression of each exon corrected for overall gene expression (exon usage) is shown for each ecotype. Exons that are differentially spliced in at least two of the three ecotypes are highlighted in blue. (e,f) Sashimi graphs highlighting patterns of differential intron excision (DIE) across all ecotypes and lakes for intron clusters associated with (e) U2AF2 and (f) SH3BGRL. The amount of DIE is measured as “change in the per cent spliced in (ΔPSI)”. The arcs represent splice junction connected exons, with the colour displaying if they are either up‐ or down‐regulated in the benthic compared to the pelagic ecotype

FIGURE 3

Sharing, connectivity and pleiotropy of differentially spliced and differentially expressed genes. (a) Venn diagrams showing the amount of overlap of differentially spliced (by analyses) and differentially expressed genes for each lake. (b) Boxplots (bar = median; notch = confidence interval around the median, box range = range between third and first quartile [interquartile range]; whiskers = extend to furthest point [highest or lowest] no further than 1.5 times the interquartile range; points = outliers) showing differences in intramodular connectivity between candidate and noncandidate genes for DEGs (n = 964 vs. n = 15,391), DEU (n = 926 vs. n = 16,355) and DIE (n = 103 vs. n = 16,355). (c) Differences in the number of associated gene ontology (GO) terms (biological processes) for candidate genes and noncandidate genes, based on DEGs and DSGs (DEU and DIE). DS genes were associated with more GO terms compared to the genomic background (noncandidate genes), suggesting higher pleiotropic effects

FIGURE 4

Functional pathways of divergence. (a) Gene ontology (GO) term interaction network for GO terms enriched for ecotype‐associated genes (RDA). Benthic ecotype associations are highlighted in orange, while pelagic ecotype associations are in blue. Central processes are highlighted by larger filled dots. (b) Multidimensional scaling (MDS) plots for shared GO terms (Biological processes; Molecular functions) enriched for differentially spliced genes (DIE and DEU). Clustering was performed based on semantic similarity of GO terms. Circles are coloured based on the number of populations they were enriched in and if they were detected based on DIE (leafcutter), DEU (dexseq) or both. The most representative and highly shared GO terms are named. The size of the circles corresponds to the number of genes associated with a GO term. (c) Venn diagrams showing the overlap between GO terms (Biological processes, Molecular functions) associated with differentially spliced (DS) and differentially expressed (DE) genes. The names of overlapping processes and functions are listed

FIGURE 5

Genetic variation underlying regulatory variation. (a) Manhattan plots showing the association of SNPs with variation in alternative splicing (intron excision ratios) across ecotypes and lakes (top) and with variation in gene expression (bottom). SNPs that are highlighted in red were detected in both analyses (n = 117). Chromosomes are highlighted by alternative colours, and unplaced scaffolds are located at the end. The blue dashed line indicates a false‐discovery rate (FDR) of 5% and the red dashed line an FDR of 1%. (b) cis‐sQTL associated with the intron excision of an intron cluster located in SRSF7. The y‐axis shows the intron excision ratio for the intron cluster by genotype (p = pelagic allele, b = benthic allele) across individuals (points) and ecotypes (colour; orange = benthic, blue = pelagic). The grey dot and ranges show the mean intron excision ratio by genotype and the standard deviation. (c) cis‐eQTL associated with the normalized gene expression of EGFR. The plot shows how the expression of EGFR differs with genotype (p = pelagic allele, b = benthic allele) across individuals (points) and ecotypes (colour; orange = benthic, blue = pelagic). The grey dot and ranges show the mean expression per genotype and the standard deviation. Both SRSF7 and EGFR are also differentially spliced or expressed in at least two ecotype pairs. (d) Predicted effects and locations of SNPs under selection. A larger proportion of SNPs under selection (in at least one lake) are located in 3′‐UTRs or are synonymous compared to proportions in the full SNP data sets (p < .001)