Spatial patterns in the contribution of biotic and abiotic factors to the population dynamics of three freshwater fish species

Studied species

We considered three species presenting different life-history strategies and patterns of elevational range shifts: the barbel (Barbus barbus), the roach (Rutilus rutilus) and the chub (Squalius cephalus) (Table 1). Following the three demographic strategies proposed by Winemiller (1992), the barbel is an ‘equilibrium’ strategist characterized by a long lifespan, a low fecundity and a large body size (Kottelat & Freyhof, 2007; Froese & Pauly, 2021), the roach is an ‘opportunistic’ strategist characterized by a small size and a low fecundity whereas the chub has the opposite characteristics and can be considered a ‘periodic’ strategist. The opportunistic strategy should maximize the colonizing capability of species in stochastic environments with frequent changes at small temporal and spatial scales. Alternatively, a periodic strategy is favored in environments with large scale cyclic variations (e.g., seasonal environment), whereas an equilibrium strategy is favored in environments with low temporal variation in habitat quality and strong biotic interactions (Winemiller, 1992).

These differences are expected to entail variations in the direction of the effect and the relative contribution of biotic and temperature-related factors to the population dynamics of the three considered species. For instance, while the abundance of all three species is expected to be positively affected by water temperature (Mills & Mann, 1985; Grenouillet et al., 2001; Piffady et al., 2010), we expect the barbel to show a stronger regulation by density, particularly for populations located at the center of the range, while the two other species are expected to be rather regulated by temperature-related factors with an increasing negative contribution as the distance to the geographic center increases (Williams, Ives & Applegate, 2003). Similarly, we expect temperature variability to have a larger positive contribution on roach abundances (opportunistic strategy) than on the abundance of the two other species where average temperature could have a stronger (and positive) influence (Winemiller & Rose, 1992).

The three species also present different patterns of elevational range shifts in the last decades (Comte & Grenouillet, 2013), with different responses observed along the elevational gradient (Table 1). For instance, the leading edges (high-elevation populations) of barbel and roach have shifted downward, whereas an upward shift was recorded for chub (Table 1). Thus, a positive influence of temperature is expected for high-elevation populations of chub while the opposite is expected for the two other species. The trailing edges (low-elevation populations) of the three species have been observed to shift upward but at a different pace. The larger shift observed for barbel (Table 1) could be associated with a stronger and negative contribution of temperature for low-elevation populations.

Species data

Fish population abundances and individual sizes were extracted from the freshwater fish monitoring database of the French Biodiversity Office (OFB, http://www.image.eaufrance.fr). We selected 182 sites where data was collected between 1982 and 2011 (3,143 sampling operations) by electrofishing. Streams were sampled by wading, during periods of low flow (from May to October), after the reproduction time, using a point sampling strategy covering the different habitats (e.g., pools, riffles, submerged vegetation) observed over the study sites (Poulet, Beaulaton & Dembski, 2011). Fish were identified to species level, measured for total body length, counted, and released to the river. For the three species, we selected time series that were composed of at least 15 years of data during which the sampling protocol remained the same and contained at least 50% of non-null captures at the population level (i.e., 0⁺ and >0⁺ individuals confounded). This selection was made to reduce the number of zeroes while keeping times series of sufficient length to allow for an appropriate estimation of the temporal dynamic of populations. We further discarded time series with more than three consecutive years missing to ensure that the loss of information in population change during the missing years is minimized (Engen et al., 2005). Missing values were ignored during the modelling process. This selection process ensures reliable model inference and left us with 71 (mean length = 17.09 years), 175 (mean length = 17.24 years) and 152 (mean length = 17.26 years) time series for barbel, chub and roach, respectively. In total, 326,234 individuals were collected. Further details about abundance and size data are provided in Appendix S1, Figs. S1 and S2.

Temperature data

Daily air temperatures from 1982 to 2011 were provided by Météo France and extracted from the high resolution (8 km by 8 km grid) SAFRAN atmospheric analysis over France (Le Moigne, 2002). Daily water temperature data measured from 2009 to 2012 at 135 sites located throughout France were provided by the French Biodiversity Office (https://ofb.gouv.fr/). From these two datasets, we used a random forest algorithm where we modelled model water temperature as a function of three covariates: air temperature, month and elevation. The model showed a very good performance and was therefore used to predict daily water temperatures for all sampling sites over the course of the study period. For further details see Chevalier et al. (2018). From daily predictions, we calculated the annual mean and intra-annual variability of water temperature between consecutive sampling occasions at each site and used these temperature variables as covariates in the population dynamic models. For each species, a summary of both variables is provided in Appendix S1; Fig. S3. These variables were transformed to z-scores before model fitting to improve model convergence.

Geographic range data

The above-mentioned abundance data do not encompass the full geographic range of the species, potentially leading to niche truncation issues and biased location of geographic range centers (Knouft, 2018; Soberón, Townsend Peterson & Osorio-Olvera, 2018; Dallas, Pironon & Santini, 2020). To obtain an unbiased estimate of the location of range centers, we used IUCN range maps (https://www.iucnredlist.org/resources/spatial-data-download). Specifically, for each species, we computed its geographic range center as the center of IUCN polygons (based on geographic coordinates) using the gCentroid function of the package rgeos (Bivand & Rundel, 2018).

Discriminating 0⁺ and >0⁺ individuals

For each species, we used the length-frequency histograms of each sampling event to separate individuals into two size classes (Bergerot & Cattanéo, 2017) using Gaussian components (McLachlan & Peel, 2004). This algorithm assumes that length data can be described by a mixture of two normal distributions which correspond in our case to the length frequency distributions of 0⁺ and >0⁺ individuals. The parameters of the two distributions were estimated using an expectation–maximization algorithm and the limit between the two size classes was fixed at the length where the two distributions crossed. However, because the algorithm performed poorly when the separation between the two size classes was not evident, e.g., when there were few individuals in each group, this procedure cannot be routinely applied to discriminate 0⁺ and >0⁺ individuals for each sampling operation. Therefore a few additional steps had to be considered (see Fig. 1; left column). For each species, we first selected 20 length-frequency histograms for which the discrimination between the two size classes was visually clear and assigned each individual to the 0⁺ or >0⁺ group based on the estimated length limit (Fig. 1). To discriminate 0⁺ from >0⁺ individuals for the remaining sampling events, we used a random forest approach (Liaw & Wiener, 2002) where individual status (0⁺ or >0⁺) was modeled as a function of individual size, individual numbers (to account for potential effects of density dependence) and annual cumulative degree-days where the water temperature was above 12 °C (i.e., the temperature below which growth is assumed inhibited; Nunn et al., 2003). The model was calibrated using the 20 sampling events for which individual status was assumed unbiased. The predictive performance of our model was tested by running a split-sample cross-validation procedure 100 times. This procedure revealed a very good model performance in predicting individual status, as measured with the Cohen’s kappa coefficient (κ > 0.99 for the three species; Cohen, 1960). We therefore used this model to predict individual’s status for the remaining sampling events. For each species, individuals in each size class were summed for each sampling event to obtain abundance time series (Appendix S1; Fig. S2).

Population dynamics model

The abundance of individuals in each size class was modelled on the log-scale using two normal distributions (one for each size class):

${\displaystyle X_{i,t}^{0^{+}}\sim Normal\left({\lambda }_{i,t}^{0^{+}},{\sigma }^{0^{+}}\right)}$ ${\displaystyle X_{i,t}^{>0^{+}}\sim Normal\left({\lambda }_{i,t}^{>0^{+}},{\sigma }^{>0^{+}}\right)}$ where X_i,t is the log-abundance of individuals in each size class (0⁺ and >0⁺) at site i and time t, λ_i,t is the expected log-abundance and σ is the associated process error variance. The means of the two distributions (i.e., ${\lambda }_{i,t}^{0^{+}}$ and ${\lambda }_{i,t}^{>0^{+}}$) were modeled with different functional forms to account for variation in the underlying demographic process (see Grenouillet et al., 2001 for a similar approach). Specifically, the dynamic of 0⁺ individuals was modelled as: ${\displaystyle {\lambda }_{i,t}^{0^{+}}={\alpha }_i^{0^{+}}+{\beta }_i^{0^{+}}\times \frac{X_{i,t}^{>0^{+}}}{log\left(S_{i,t}\right)}+\sum _{j=1}^J{\gamma }_{i,j}^{0^{+}}\times U_{j,i,t}+log\left(S_{i,t}\right)}$ where ${\alpha }_i^{0^{+}}$ is a site varying intercept, ${\beta }_i^{0^{+}}$ is a density-dependent parameter reflecting the dependency of 0⁺ individuals to the abundance of >0⁺ individuals at time t and can be interpreted as the apparent recruitment rate (since sampling takes place after reproduction) and S_i,t is the sampling area (offset term). Note that the recruitment rate is only ‘apparent’ because (1) the >0⁺ size-class includes both spawners and non-spawners (reported age at maturity for females is 3–4 years for the chub, 2–3 years for the roach and 5 years and more for the barbel; Keith et al., 2011; Keith et al., 2020) and thus also accounts for the effect of competition with non-spawners and (2) inferences are based directly on the abundance of recruits (i.e., 0⁺ individuals that successfully hatched), hence not accounting for variation in the per capita reproductive investment (i.e., the number of eggs laid by a given individual). A more accurate estimation of the recruitment rate could be achieved using egg data together with data on spawner biomass/abundance. The parameters ${\gamma }_{i,j}^{0^{+}}$ are regression coefficients applied to the array U_j,i,t which contains the raw and the squared values of the mean and the variability of water temperatures at site i and time t. Thus, ${\gamma }_{i,j}^{0^{+}}$ is a vector of coefficients representing the linear and the quadratic effect of the mean and the variability of water temperature on 0⁺ abundance at each site.

The dynamic of >0⁺ individuals was represented using a modified version of the stochastic Gompertz model of population growth as: ${\displaystyle {\lambda }_{i,t}^{>0^{+}}={\alpha }_i^{>0^{+}}+X_{i,t-1}^{>0^{+}}+{\beta }_i^{>0+}\times \frac{X_{i,t-1}^{>0^{+}}}{log\left(S_{i,t-1}\right)}+{\delta }_i^{>0+}\times \frac{X_{i,t-1}^{0^{+}}}{log\left(S_{i,t-1}\right)}+\sum _{j=1}^J{\gamma }_{i,j}^{>0^{+}}\times U_{j,i,t}+log\left(\frac{S_{i,t}}{S_{i,t-1}}\right)}$ where ${\alpha }_i^{>0^{+}}$ is a site varying intercept representing the intra-class productivity rate (an analog to the population growth rate with values above one indicating positive productivity rates), ${\beta }_i^{>0^{+}}$ is a density-dependent parameter representing the competition between >0⁺ individuals for access to resources and ${\delta }_i^{>0^{+}}$ represents the transition probability between the two size classes (i.e., the apparent survival rate of 0⁺ individuals). Similar to the recruitment rate, we note here that the survival rate is only ‘apparent’ because the >0⁺ size-class includes individuals in different ages. This survival rate thus also accounts for the survival probability of all other size-classes. The parameters ${\gamma }_{i,j}^{>0^{+}}$ are regression coefficients representing the linear and the quadratic effects of the two temperature variables on >0⁺ abundance.

Quadratic effects were included in both dynamics to account for potential bell-shaped response curves along the temperature gradient (Austin, 1999). The model was fitted to each species separately, and included random site effects for all population dynamic parameters, ultimately making it possible to analyze spatial patterns in the contribution of biotic and temperature-related factors to the dynamic of each size-class.

Parameter estimation and model goodness of fit

The model was fitted to each species using Bayesian inference and weakly informative priors. Site-specific parameters (${\alpha }_i^{0^{+}},{\beta }_i^{0^{+}},{\gamma }_{i,j}^{0^{+}},{\alpha }_i^{>0^{+}}$, ${\beta }_i^{>0^{+}}$, ${\delta }_i^{>0^{+}}$, ${\gamma }_{i,j}^{>0^{+}}$) were assumed to follow normal distributions with a vector of means µ{μ_α⁰⁺, μ_β⁰⁺, ${\mu }_{{\gamma }_j^{0^{+}}}$, μ_α^>0+, μ_β^>0+, ${\mu }_{{\gamma }_j^{>0^{+}}}$, μ_δ^>0+} and of standard deviations σ{ σ_α⁰⁺, σ_β⁰⁺, ${\sigma }_{{\gamma }_j^{0^{+}}}$, σ_α^>0+, σ_β^>0+, ${\sigma }_{{\gamma }_j^{>0^{+}}}$, σ_δ^>0+}. The vector µrepresents the average value of the parameters across all sampling sites whereas the vector σ represents departures from the mean and therefore the spatial variability in parameter values. We used normal distributions with mean zero and standard deviations of 10 as priors for all µ. For σ, σ⁰⁺ and σ^>0+, we used half-Cauchy distributions (Gelman, 2006). For each species, we generated three chains of length 11,000 with initial values in different regions of parameter space and discarded the first 1,000 iterations as burn-in. Chains were sampled every 10 iterations. Convergence was visually assessed and confirmed using the Gelman and Rubin statistic with a threshold value of 1.1 (Gelman & Rubin, 1992). Highest Posterior Density (HPD) intervals were used as 95% credible intervals. For each parameter, differences between species were assessed by computing the proportional overlap between the two posterior distributions. A low overlap (threshold set to 5% meaning that only 5% of MCMC samples were common between the two distributions) was taken as evidence that estimated parameters were different between species.

We used posterior predictive checks (Gelman, Meng & Stern, 1996) to assess the goodness of fit of our model for the three species. Specifically, we used χ² discrepancy metrics to compute the posterior predictive p-value, which quantifies the extent to which the proportion of samples in which the distance of observed data to the model is greater than the distance of replicated data to the model. Values close to 0.5 indicate a good model fit, whereas values close to 0 or 1 indicate lack of fit. Bayesian p-values were calculated regarding the log-abundance of both 0⁺ and >0⁺ individuals. We fitted the models using JAGS 4.3.0 (Plummer, 2003), run through the R environment (R Core Team, 2019) using the packages R2jags (Su & Yajima, 2013) and rjags (Plummer, 2014). The JAGS code is available in Appendix S2.

Elasticity analyses

We applied an elasticity analysis to species model outputs in order to highlight the relative contribution of biotic (${\alpha }_i^{0^{+}},{\beta }_i^{0^{+}},{\beta }_i^{>0^{+}}$ and ${\delta }_i^{>0^{+}}$) and temperature-related (${\gamma }_{i,j}^{0^{+}}$ and ${\gamma }_{i,j}^{>0^{+}}$) factors to the population dynamics of the three species following the framework developed by Koons et al. (2015). Specifically, for each species, we used the median of the posterior distribution of parameters obtained from the fitted model to project the log-abundance of both size-classes at each site over the study period (Haridas, Tuljapurkar & Coulson, 2009). This was done by iteratively updating the parameter λ_i,t for both 0⁺ and >0⁺ individuals using observed predictor values for the parameter for which elasticity needs to be calculated but average predictor values for the other parameters. We measured the contribution of biotic and temperature-related factors separately for each size-class by comparing the log-abundance computed using parameter values predicted by the model (θ_ori) to the log-abundance computed by changing each parameter value, one at a time by 10% (θ_per; both the linear and the quadratic terms were changed for temperature-related factors). Specifically, elasticities were computed numerically as: ${\displaystyle {\mathrm{e}}_{\omega ,\mathrm{i},\mathrm{t}}=\frac{{\theta }_{\mathrm{per},\mathrm{i},\mathrm{t}}-{\theta }_{\mathrm{ori},\mathrm{i},\mathrm{t}}}{{\theta }_{\mathrm{ori},\mathrm{i},\mathrm{t}}}\times \frac{1}{\delta }}$ where θ is the response parameter (i.e., log-abundance of the considered size-class original and perturbed) at site i and time t, ω is the parameter of interest and δ is the proportional change in ω (i.e., 10%). The mean of all e_ω,i,t, therefore represents the estimated elasticity of parameter ω at the species level (e_ω,sp). For each parameter, differences between species were tested using Wilcoxon signed-rank tests with p-values adjusted for multiple comparisons using the Bonferroni correction.

Spatial patterns in the contribution of biotic and temperature-related factors

In order to establish how the contribution of biotic and temperature-related factors varied spatially, we built four different linear models: an intercept-only model (null model against which the other models are compared to), a model with elevation as covariate (elevation range-shift model), a model with latitude as covariate (latitude range-shift model) and a model with the Haversine distance (i.e., Euclidean distance accounting for the curvature of the Earth) to the geographic range center (abundant-center hypothesis model; Soberón, Townsend Peterson & Osorio-Olvera, 2018) as covariate. The three models with covariates included quadratic terms to account for non-linear effects. Models were run for the eight parameters for which elasticity was computed and compared using AIC (Burnham & Anderson, 2002). For each species, further details for the three covariates is provided in Appendix S1; Fig. S4.