Genomic changes underlying repeated niche shifts in an adaptive radiation

Abstract

In adaptive radiations, single lineages rapidly diversify by adapting to many new niches. Little is known yet about the genomic mechanisms involved, that is, the source of genetic variation or genomic architecture facilitating or constraining adaptive radiation. Here, we investigate genomic changes associated with repeated invasion of many different freshwater niches by threespine stickleback in the Haida Gwaii archipelago, Canada, by resequencing single genomes from one marine and 28 freshwater populations. We find 89 likely targets of parallel selection in the genome that are enriched for old standing genetic variation. In contrast to theoretical expectations, their genomic architecture is highly dispersed with little clustering. Candidate genes and genotype-environment correlations match the three major environmental axes predation regime, light environment, and ecosystem size. In a niche space with these three dimensions, we find that the more divergent a new niche from the ancestral marine habitat, the more loci show signatures of parallel selection. Our findings suggest that the genomic architecture of parallel adaptation in adaptive radiation depends on the steepness of ecological gradients and the dimensionality of the niche space.

Keywords: Adaptive radiation; Haida Gwaii; genomics; niche shift; niche space; threespine stickleback.

Figure 1

Adaptive radiation of threespine stickleback on Haida Gwaii, Canada. (a) Geography of the 28 freshwater populations, with habitat, gigantism, and plate morph indicated by symbol shapes and border colors. (b) Niche occupation of 28 freshwater and one marine stickleback populations along three major ecological axes: predation regime, light spectrum, and ecosystem size (Reimchen et al. 2013). Similar niches in this multidimensional niche space were colonized independently in different watersheds. Predation regime: invertebrate (0: purple) or vertebrate dominated, with main predators either being cutthroat trout Oncorhynchus clarkii (1: orange) or rainbow trout O. mykiss (2: blue). Stickleback drawings depict typical morphologies of selected populations, representing extremes in morphospace (Reimchen et al. 2013). Axis units: light spectrum = percent light transmission at 400 nm wavelength, ecosystem size = log₁₀‐transformed lake area in hectares. (c) Whole mitochondria maximum‐likelihood phylogeny of Haida Gwaii stickleback reveals the presence of Japan Sea stickleback (Gasterosteus nipponicus) derived haplotypes in three freshwater populations. Outgroups are blackspotted stickleback (G. wheatlandi) and nine‐spined stickleback (Pungitius pungitius). Branch labels: bootstrap support (%) shown if >80% support. “Broken” branches are long branches shortened for better visualization of the Haida Gwaii haplotypes. (d) Genomic variation mirrors geography along a SW‐NE gradient in a principal component analysis of genome‐wide autosomal SNPs (gPC). Square brackets contain percentages variance explained by each principal component. See Table S1 for full population names.

Figure 2

Dispersed genomic architecture of parallel adaptation and enrichment for old genetic variation. (a) Signatures of selective sweeps across the stickleback adaptive radiation on Haida Gwaii are distributed throughout the genome. Vertical gray bars indicate outlier regions enriched for both iHS (blue) and H12 (green) outliers (black points) at Bonferroni‐corrected alpha <0.001, with iHS and H12 outliers exceeding 99.9% quantiles of neutral demographic expectations. Several blocks of old standing variation, indicated by a large range in absolute divergence (d _XY) between population pairs (purple area), overlap such outlier regions. Color codes below each chromosome indicate d _XY range quantiles shown in panel c. Roman numerals are chromosome names; letters correspond to outlier regions (see Table S2). (b) Genotype distribution of Haida Gwaii stickleback for known large inversions: “marine” haplotypes are found in several freshwater populations, in particular in fully plated freshwater populations (DA, DW, SY). First principal components of SNP genotypes in the indicated genomic interval are shown, with the percentage of variation explained in brackets. (c) Distribution of d _XY ranges in 100‐kb windows based on all pairwise d _XY values between the 29 individuals, with quantiles indicated by colors. (d) Outlier regions are enriched for old standing genetic variation: 16% of the outlier regions fall into genomic regions containing old genetic variation.

Figure 3

Linkage disequilibrium (LD) between outlier regions indicate correlated evolution between physically unlinked outlier regions. (a) LD between outlier regions sorted by chromosome and physical position (above diagonal) and into 47 clusters in a hierarchical cluster analysis (below diagonal, see cluster tree on the right), with lines connecting outlier regions and colored lines corresponding to the 47 groups. Colored r ² values exceed the top 1% LD values among random SNPs on different chromosomes (gray distribution). (b) Observed mean interchromosomal LD between 89 outlier regions exceeds interchromosomal LD between random SNPs estimated from 10,000 permutations.

Figure 4

Niche divergence from the ancestral habitat and geography best predict the genomic architecture of parallel adaptation. (a) Number of sweep haplotypes in a population shown by color code in the three‐dimensional niche space of Haida Gwaii archipelago. (b) Niche divergence from the ancestral, marine habitat has the strongest effect (Table 1) on the number of sweep haplotypes in a population, followed by geography (gPC1, gPC2), likely due to spatial autocorrelation of the niche space in the Haida Gwaii archipelago (Reimchen et al. 2013). Population size N _e is negatively associated with haplotype number, as the marine and large freshwater populations with marine‐like phenotypes contain few sweep haplotypes, contrary to expectations from population size limiting the efficiency of selection. (c) Number of shared sweep haplotypes between populations shown by color code in the niche space of the Haida Gwaii archipelago. (d) Mean niche divergence from the marine habitat has also the strongest effect on the number of shared sweep haplotypes between population, followed by genetic relatedness of populations and differences in population size ΔN _e and pairwise niche divergence (Table 1).

Figure 5

Ecology‐associated outlier regions and their distribution across the genome. Associations between outlier regions (vertical gray bars) and the three major ecological axes are widely distributed across the genome. Shown are putative candidate gene targets of selection (upward pointing triangle), GWAS associations of peak H12‐SNPs with 50‐kb sliding window averaged P‐values <0.01 (circle), and random forest associations for SNPs with variable importance mean decrease in accuracy of >20 (downward pointing triangle). See Methods for grouping of ecological properties and phenotypic traits into associations with light spectrum, predation regime, and ecosystem size.