Prolonged migration windows buffer migratory animal populations against uncertainty in resource availability. Understanding how intensifying droughts from climate change influence the migration window is critical for biodiversity conservation in a warming world. We explored how drought affects the seaward migration of endangered coho salmon (Oncorhynchus kisutch) near the southern extent of their range in California, USA. We tracked stream departures of juvenile coho, measuring streamflow and temperature in 7 streams over 13 years, spanning an historic drought with extreme dry and warm conditions. Linear mixed effects models indicate that, over the range of observations, a decrease in seasonal streamflow (from 4.5 to 0.5 mm/day seasonal runoff) contracted the migration window by 31% (from 11 to 7 weeks). An increase from 10.2 to 12.8 ℃ in mean seasonal water temperature hastened the migration window by three weeks. Pacific salmon have evolved to synchronize ocean arrival with productive ocean upwelling. However, earlier and shorter migration windows during drought could lead to mismatches, decreasing fitness and population stability. Our study demonstrates that drought-induced low flows and warming threaten coho salmon in California and suggests that environmental flow protections will be needed to support the seaward migration of Pacific salmon in a changing climate.

Runoff and water temperature data

We estimated mean annual precipitation, averaged across each drainage area, using Google Climate Engine, March 2011 - February 2021.

Where multiple temperature loggers were present in a study stream, we selected a single location based on the completeness of data in the study season and proximity to the PIT antenna. Hourly temperature measurements were converted into mean daily values. 

Analysis

For data analysis and modeling, we excluded streams that had less than 3 years of biological data, leaving 47 stream-years. We conducted all analyses in R (version 4.0.4, R Core Team, 2018). We tested outmigration timing data for normal distribution among streams, years, and stream-years, using the shapiro.test function of the broom package. The Shapiro-Wilk test showed that all distributions were unlikely to be normally distributed (i.e. among years, p = 5.5 × 10^-9 – 7.6 × 10^-39 and W = 0.88 – 0.98). However, the Shapiro-Wilk test can provide small p-values for large samples and consequently provide a false negative, regarding normal distribution (among years, sample size range was 485 – 3453).

Therefore, we could not rule out the possibility of assumptions being met to perform ANOVA (Analysis of Variance) tests. We did so, using the aov function of the AICcmodavg package: i) one-way, by stream, ii) a one-way, by year, iii) a two-way, by stream and year, and iv) a two-way with stream-year interaction. To isolate the effects of stream and year on variance, we performed the ANOVA tests on the maximum subset of data for which each stream had the same years of outmigration (four streams, each with the same six years of data, totaling 24 stream-years). The aictab function of the AICcmodavg package demonstrated that the two-way model with stream-year interaction was the highest performing (lowest AICc value), followed by: the two-way model, one-way by year model, and one-way by stream model. In both ANOVA tests, the year, stream, and year-stream interaction terms each had “Pr(>F)” values < 2 × 10^-16. The “2-way ANOVA with interaction” (year F-value 646.58, stream F-value 349.85, year-stream interaction F-value 29.31, residuals 4.11 × 10^-16) had higher F values and lower residuals than the 2-way ANOVA (year F-value 629.3, stream F-value 340.5, residuals 4.22 × 10^-16). We used the TukeyHSD function of the AICcmodavg package to conduct pairwise tests for significant differences in outmigration timing distributions. Among streams, five of six pairwise differences were highly significant (p < 0.0001). Among years, all 15 pairwise comparisons were highly significant (p < 0.001). Among stream-years, 216 of 277 pair-wise comparisons were significant (p < 0.05). We checked for homoscedasticity in the interaction model, using the leveneTest function of the car library, and we found evidence that the variance across groups is significantly different. Consequently, we cannot assume homogeneity of variances in the different groups, which is typically a required assumption for conducting ANOVA tests.

Since the normal distribution assumption of the one-way ANOVA was not met, we applied the Kruskal-Wallis test, as a non-parametric alternative to test for variance among streams and years, using the package rstatix. As with the ANOVA tests, we performed Kruskal-Wallis tests on the maximum subset of data for which each stream had the same years of outmigration (24 stream-years), using the functions kruskal_test, kruskal_effsize, dunn_test, and wilcox_test. Among streams, we found significant variance (p = 2.16 × 10^-143), with a “small” effect size (eta-squared measure = 0.04) (Tomczak and Tomczak 2014), and 5 of 6 pairwise differences were highly significant (Dunn’s test & Wilcoxon’s test: p < 0.0001). Among years, we found significant variance (p = 0), with a “large” effect size (eta-squared measure = 0.17) (Tomczak and Tomczak 2014), and 13 of 15 pairwise differences were highly significant (Dunn’s test & Wilcoxon’s test: p < 0.0001).

Modeling the effects of streamflow and water temperature on outmigration timing

Modeling was limited to the 42 stream-years for which water temperature and outmigration timing data were collected. For the outmigration start date model, the runoff date range was March-April and the degree-days date range was March-April. For the outmigration end date and duration models, the runoff date range was March-June and the degree-days date range was March-April. Coefficient units are “days per daily runoff (mm)” and “days per 100 degree-days”.

In identifying top model(s), we did not consider degree-days to influence outmigration duration because: i) the AIC value of the runoff-only model was 1.99 less than the additive model, ii) the degree-days in the additive model had a p-value > 0.05, and iii) Mar-Jun runoff had similar coefficient effect sizes in the additive model and run-off only model (Appendix S1: Table S3). We calculated conditional coefficients (including stream, as a random effect) and marginal coefficients (excluding stream, as a random effect) of determination (R²) (Nakagawa and Schielzeth 2013), using the r.squaredGLMM function of the MuMIn package (Barton` 2020). We also reported the model coefficients and 95% confidence intervals, as measures of effect size, and generated partial dependence plots for using the plot_model function of the sjPlot package (Lüdecke 2021).

Literature cited

Barton`, K. (2020). MuMIn: Multi-Model Inference. R package version 1.43.17.

Lüdecke, D. (2021). sjPlot: Data Visualization for Statistics in Social Science. R package version 2.8.9.

Nakagawa, S., and H. Schielzeth. 2013. A general and simple method for obtaining R 2 from generalized linear mixed-effects models. Methods in Ecology and Evolution 4:133–142.

Tomczak, M., and E. Tomczak. 2014. The need to report effect size estimates revisited. An overview of some recommended measures of effect size 1:7.

Please see DataS1/data/README_Metadata.pdf.