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).