Results

Genomic signatures of local adaptation to round goby invasion and water calcium.

Using the core model in Baypass as our outlier analysis, we found 226,794 outliers SNPs with a p-value < 0.001 and either high or low XtX values, which represented ≈1.1% of the dataset (Fig. S1A). Outlier SNPs with high XtX values can be interpreted as putatively under positive selection, while low XtX values indicate balancing selection (Gautier, 2015). We also investigated the association of SNP allele frequencies with the selected environmental variables (invasion status and calcium concentration) using the STD model in Baypass. We found 778 outlier SNPs significantly associated with calcium concentration across the three replicate runs (BFis > 20, ≈ 0.004% of the dataset) and 88,277 SNPs associated with the goby presence/absence (q-value < 0.01, ≈ 0.4% of the dataset). Using the poolFreqDiff analysis with a FDR of 1%, we identified 54,285 outliers displaying consistent differences in allele frequency in the same direction between the low and high calcium habitats (0.3% of the dataset) and 23,651 outlier SNPs between the populations from the invaded and uninvaded habitats (0.1% of the dataset; Fig. S5). Of those calcium concentration outliers, 18 were in common between the Baypass STD model and the poolFreqDiff analysis, whereas 3,324 of the predator status outliers were in common between the Baypass STD model (using the C2-statistic) and the poolFreqDiff analysis (Fig. 2, Fig. S5). Most of the outliers were uniquely associated with a single covariable, with 2,668 SNPs associated with both calcium and invasion status across methods (Fig. S5). Overall, we found 1,050 SNPs in common between the Baypass core and STD models, as well as 7,009 SNPs in common between the Baypass STD model and the poolFreqDiff analyses including both the invasion status and calcium concentration (Fig. S6).

(Mal)adaptive responses to round goby invasion and water calcium levels

We found life history differences between the gastropod populations from the two environments, using fecundity and survival as fitness components (Fig. 3). For fecundity, the model with the random effect of population origin was not better than the model without (ΔAICc = 1.2), and only the origin water (Ottawa River: OR, LCGA vs St. Lawrence River: SL, HCGP) effect was significant (p = 0.015). Even though the interaction between origin water and treatment water was not significant (p = 0.076), fecundity was higher for SL and OR populations in home water (13.30 SD 19.70 and 2.42 SD 6.27 respectively) than in transplant water (3.00 SD 4.33 and 0.92 SD 1.83 respectively). SL populations produced ≈ 5 times more eggs than OR populations (4.90, 95% CI [1.54-15.58]). For survival, the model with a random effect of the population was better than without (ΔAIC = 250.7). The fixed effects of origin and treatment water were significant (p = 0.020 and p = \(1.496\times 10^{-9}\), respectively), but their interaction and the goby cue effect were not (p = 0.203 and p = 0.794, respectively). SL populations were 10 times more likely to survive compared to OR populations (9.55, 95% CI [1.70-53.73]). However, exposure to treatment water from the St. Lawrence River significantly lowered the odds of survival, with survival rates less than one-third that of populations exposed to Ottawa River water (0.31, 95% CI [0.21-0.46]). There was considerable variation in survival rates between populations, as demonstrated by the significant effect of population on survival. Variation in survival among populations (random effect of the population of origin) did not depend on geographical location or habitat of origin (Fig. S7): OKA-LCGA, PG-HCGP, and PON-HCGP had significantly higher survival rates, while BEA-HCGP, PDC-HCGA, and PST-HCGP had significantly lower survival rates.

Demographic and genetic effects of the invasion

Genome-wide nucleotide diversity π was relatively high overall with 0.011 (SD 0.006) on average. Diversity was similar between the populations from the invaded St. Lawrence River habitats (0.011, SD 0.007) and the populations of the uninvaded Ottawa River (0.011, SD 0.006), with a negligible effect size of habitat type (Fig. 4A; Hedges’g = 0.09, 95% CI [0.16, 0.02]). Estimates ofθ Watterson were identical between the two habitat types (Fig. 4B; SL populations: 0.015 SD 0.008; OR populations: 0.015 SD 0.008), and the effect size of habitat was therefore negligible (Hedges’g = 0.06, 95% CI [0.13, 0.01]). This resulted in a slightly negative overall Tajima’s D (-0.38 SD 0.43), which was lower for SL populations (-0.41 SD 0.44) than for OR populations (-0.36 SD 0.43), but the effect size of the difference was negligible (Hedges’ g =0.11, 95% CI [0.10, 0.13], Fig. 4C). Observed heterozygosity was not significantly different (p-value = 0.289, t = -1.153, df = 6.5; Fig. 4D) between the populations from uninvaded (0.168 SD 0.003) and invaded sites (0.171 SD 0.003).\(\frac{}{}\)
Consistent with the known population declines in invaded habitats, we found a significant association between invasion status and the effective population size of our study populations based on our demographic modeling (Table 1, Fig. S9). For the pair PB-LCGA (refuge) and PG-HCGP (invaded), the best model included bottlenecks in both populations, but the inferred parameters indicated that only the invaded population showed Ne recovery (Table S2), suggesting an alternative model that we did not test (bottlenecks in both population but Ne recovery only in the invaded population). In the case of IPE-LCGA and BEA-HCGP, we did not detect a significant bottleneck, but the invaded population had an effective population size 40 times lower compared to the refuge (Table S2). Finally, the best model for PDC-HCGA and GOY-HCGP included a bottleneck without recovery only for the invaded population (Table 1).

Isolation by environment and variable gene flow between habitat types

We found that populations clustered by environment type (Fig. 5), particularly by the presence/absence of the round goby, with RAF-LCGP and PDC-HCGA both clustering with the invaded populations and showing lower pairwise FST values within those clusters. Our results from the scaled covariance Ω matrix of population allele frequencies inferred with Baypass indicate that there is also positive covariance in allele frequencies within clusters and negative covariances between clusters. These population structure results were concordant between the Ω matrix (Fig. 5A) and the pairwise FST matrix (Fig. 5B). Population structure was also explained by isolation by distance and by the environment (Fig. S10), as the positive correlation between the genetic distance FST/(1-FST) and the log of the geographic distance or the Mahalanobis were significant (Mantel test, 9999 permutations: p = 0.007, r2 = 0.209 and p = 0.044, r2 = 0.078 respectively). We found low but significant asymmetric gene flow in most cases with the scaled migration rates 0.05 < 2Nem < 5.5 (Table S2). However, gene flow was non-significant (2Nem < 0.05) from the low calcium refuge toward high calcium invaded populations (IPE-LCGA toward BEA-HCGP and PB-LCGA-HCGP toward PG). For these two population pairs, the lower migration rates were consistent with higher values of pairwise FST and time since the population split (Fig. 5 and Table S2). Finally, the population pair GOY-HCGP and PDC-HCGA (wetland refuge) had much higher migration rates in both directions compared to the other population pairs, consistent with lower pairwise FST and these two populations belonging to the same cluster in Baypass (Fig. 5 and Table S2).