Temperature-dependent Allee effects in a stage-structured model for Bythotrephes establishment

Abstract

Whether the invasive freshwater cladoceran Bythotrephes longimanus can establish after introduction into a water body depends on several biotic and abiotic factors. Among these, water temperature is important because both development rates and mode of reproduction (parthenogenetic or sexual) in Bythotrephes are influenced by temperature. We built a stage-structured model for the population dynamics of Bythotrephes based on the temperature-dependency of events in its life cycle and used the density of resting eggs at the end of each year to track changes in population density. The model was parameterized using data from published laboratory experiments and data on the Bythotrephes population in Harp Lake, Canada, from 1994 to 2005. The parameterized model was then used to simulate the outcome of invasions with different initial resting egg densities under different temperature regimes. A strong Allee effect emerged from the model, i.e. there is a critical threshold density above which the population can establish and below which it goes extinct. We showed analytically that the existence of an Allee effect arises from the model structure and is therefore robust to the parameter values. An increase in temperature reduces the establishment threshold for introductions in the same year as well as for introductions in the previous years. We therefore hypothesize that climate warming might facilitate Bythotrephes invasions. Finally, we study how the establishment threshold is influenced by the timing of the introduction event and thus identify time periods during the year when lakes may be particularly susceptible to Bythotrephes invasions.

Appendix: Model parameterization details

Data on the Bythotrephes population of Harp Lake

Harp Lake, Ontario, (45°23′N, 79°07′W) is an inland lake on the Canadian Shield and has an area of 71.4 ha and a maximum depth of 37.5 m (Hall and Yan 1997). Bythotrephes was first detected in 1993 during routine monitoring by the Ontario Ministry of the Environment (Yan and Pawson 1997), but could have been present at low densities in the previous 1 or 2 years. The data we used for parameter estimation span the years 1994 to 2005. In the years 1994, 1995, and 1997 to 2004, following a stratified random sampling regime, on each sampling date, one randomly chosen location in each of three depth zones in each of three lake regions was sampled. Additionally, a sample was always taken at the deepest position of the lake, resulting in a total of 10 sample locations per sampling date (Hall and Yan 1997). Samples were taken as vertical hauls of a closing tow net with a diameter of 0.75 m, a length of 2.5 m and a 285-μm mesh (Yan and Pawson 1997). For each Bythotrephes individual found, instar number, gender, shape of the tail spine, state of the brood and clutch size were determined and stored in a database, with some exceptions: In 1994 and 1997 gender was not determined, so that the variables F ₁ and F ₂ also include first and second instar males, respectively. For these years, we only use the brood size data, not the densities of each stage. In general, it is difficult to distinguish first instar males from females, which might be a reason for the low observed densities of first instar males.

These data were used to create time series of the densities of the model variables F ₁, F ₂, F ₃, M ₁, and M ₂. We calculated densities for each station separately as number of individuals within a stage divided by the sample volume (corrected for filtration efficiency) in m³. The model variables F ₁ and F ₂ consist of sexually or parthenogenetically produced females of instars one and two, respectively. F ₃ is a combination of third instar sexually produced females, which need to molt once more to reproduce, mature fourth instar sexually produced females, and mature third instar parthenogenetically produced females. The density of resting eggs cannot be assessed by sampling the water column. There have been attempts to monitor the number of resting eggs in the sediment (Yurista 1997), but this was not done for Harp Lake during the time period of interest. We averaged the densities of the other stages at all 10 sampling sites to get the average density for each sampling date. The number of sampling dates throughout the growing season differed from year to year, ranging from only 4 in 1999 and 2000 to 25 in 2003. We also extracted information on clutch sizes and mode of reproduction from the Harp Lake Bythotrephes database. In 2003, there were additional horizontal tows at different depths that can be used to characterize the vertical distribution of the Bythotrephes population (Young and Yan 2008).

Temperature profiles from the lake surface to the bottom were measured at the deepest station in Harp Lake by the Ontario Ministry of the Environment on each visit. To obtain the temperature values for depth ϕ on a certain day of the year t, T(t, ϕ), we linearly interpolated between the measured time points and depths.

Accounting for the vertical distribution of Bythotrephes in the water column

At any moment in time, the individuals of a Bythotrephes population are distributed over a range of water depths and thus experience different temperatures, which leads to differences in the rates of temperature-dependent life-history processes. Therefore, we weighted all functions in our model that vary with temperature (e.g., development rates) by the probability of experiencing a certain temperature, which can be constructed from the vertical distribution and the temperature profile. The vertical distribution differs between the life-history stages and varies both seasonally and from day to night (Young and Yan 2008). For any time point and life-history stage, the cumulative Bythotrephes proportion Pr(ϕ ≤ x), i.e. the probability of the depth ϕ of an individual to be smaller or equal to x, can be well described by a function of the form:

$$ Pr(\phi \leq x)=\frac{k}{1+\exp(-a(x-x_c))}, $$

(17)

where k, a, and x _c are fitted parameters. By differentiating and normalizing the cumulative density function, we obtained the corresponding density function that gives the probability density of an individual to be located at depth x:

$$ \psi(x) = \frac{{ka\cdot\exp ( - a(x - x_{c} ))}}{{1 + \exp ( -a(x - x_{c} )) \cdot K}}, $$

(18)

where $K= \int_{0}^{{\phi _{{\max }} }} {\frac{{ka\cdot\exp ( - a(x - x_{c} ))}}{{1 + \exp ( - a(x - x_{c} ))}}dx}$ is a normalization constant that ensures that the probability of an individual to be located somewhere between the lake surface and its maximum depth is 1.

Young and Yan (2008) obtained cumulative density functions by fitting Eq. 17 to the distributions measured in Harp Lake in 2003. From their fitted parameters, we constructed probability distributions for day and night for first instar females (Pr ⁽¹⁾ _day/night ), second instar females (Pr ⁽²⁾ _day/night ), third instar females (Pr ⁽³⁾ _day/night ), and males (Pr ^(m) _day/night ). Young and Yan (2008) also modeled the distribution of females that are carrying brood, but since their vertical distribution did not differ from other third instar females, we will use Pr ⁽³⁾ _day/night for all processes involving third instar females. To evaluate the time-dependent population mean, f(t), for a given function of temperature, g(T) (e.g. the development rates), we weighted the value of the temperature-dependent function at a given depth and its specific temperature with the probability density of a Bythotrephes individual of the stage of interest to be at this depth and integrated this expression over all relevant depths. Additionally, we weighted the values of the integral for day and night by the proportions of day and night, p _day and 1 − p _day , respectively, which are functions of the time of the year (analogous to equation B2 in Young et al. 2009):

$$ \begin{aligned} \underbrace {{f(t)}}_{{{\text{population}}\,{\text{function}}\,{\text{value}}}} = \mathop {\overbrace {{\underbrace {{p_{{day}} (t)}}_{{{\text{proportion}}\,{\text{day}}}}\cdot\mathop {\mathop \int \limits^{{\phi _{{max}} }} }\limits_{0} \underbrace {{g(T(t,x))}}_{{{\text{function}}\,{\text{value}}\,{\text{at}}\,{\text{depth}}\,{\text{x}}}}\cdot\underbrace {{\psi_{{day}} (x)\,dx}}_{{{\text{daytime}}\,Bythotrephes\,{\text{probability}}\, {\text{density}}\,{\text{at}}\,{\text{depth}}\,{\text{x}}}}}}^{{{\text{value}}\,{\text{of}}\,{\text{the}}\,{\text{depth - weighted}}\,{\text{function}}\,{\text{during}}\,{\text{the}}\,{\text{day}}}}}\limits_{{}} + \mathop {\overbrace {{\underbrace {{(1 - p_{{day}} (t))}}_{{{\text{proportion}}\,{\text{night}}}}\cdot\mathop {\mathop \int \limits^{{\phi _{{max}} }} }\limits_{0} \underbrace {{g(T(t,x))}}_{{{\text{function}}\,{\text{value}}\,{\text{at}}\,{\text{depth}}\,{\text{x}}}}\cdot\underbrace {{\psi_{{night}} (x)\,dx}}_{{{\text{nighttime}}\,Bythotrephes\,{\text{probability}}\, {\text{density}}\,{\text{at}}\,{\text{depth}}\,{\text{x}}}}}}^{{{\text{value}}\,{\text{of}}\,{\text{the}}\,{\text{depth - weighted}}\,{\text{function}}\,{\text{during}}\,{\text{the}}\,{\text{night}}}}}\limits_{{}} \end{aligned} $$

(19)

Since the probability distributions have been parameterized only for a number of sampling dates, we linearly interpolated f(t) between these sampling dates. The vertical distribution of Bythotrephes has been studied in Harp Lake only for 2003. We assumed that the vertical distribution as a function of the day of the year will be similar in other years.

Development rates, d _e (t), d _f1(t), d _f2(t), d _m1(t)

Lehman et al. (1997) found that the relationship between temperature and the embryonic development rate can be described by three models that perform equally well over an ecologically relevant range of temperatures. Therefore, we chose their first and most simple model

$$ \frac{1}{D}=a\cdot T+b $$

(20)

to parameterize the development rates in our dynamical system. It describes the inverse of the development time D, i.e. the development rate, as a linear function of the temperature, T. For the development time from the embryonic red-eye stage, when the eye pigment first becomes visible, to the release as neonate, a and b are 0.0024 and 0.01151, respectively (Lehman et al. 1997). We assumed with Lehman et al. (1997) that all other development times scale linearly with the embryonic development time and used their scaling factors. Restricting the development rates to positive values, we obtained for the development rate between male and female instars 1 and 2:

$$ \delta_{1}(T)=\max(0.022879\left(\rm{^\circ \hbox{C\,day}}\right)^{-1}\cdot T-0.12908\rm{day}^{-1},0). $$

(21)

The development rate between female instar 2 and 3 is:

$$ \delta_{2}(T)=\max(0.02464\left(\rm{^\circ\hbox{C\,day}}\right)^{-1}\cdot T-0.13902\rm{day}^{-1},0). $$

(22)

The time from beginning of instar 3 to primiparity that Lehman et al. (1997) measured served as approximation for the total embryonic development time in our model, from the deposition of eggs in the brood pouch to their release. Using the appropriate scaling factor, we obtain:

$$ \delta_{e}(T)=\max(0.018991\left(\rm{^\circ\hbox{C\,day}}\right)^{-1}\cdot T-0.10715\rm{day}^{-1},0). $$

(23)

All three development rates are for an individual Bythotrephes that is experiencing temperature T; so we used Eq. 19 with g(T) = δ _i (T) and the corresponding stage-specific vertical distributions Pr ⁽³⁾ _day/night , Pr ⁽¹⁾ _day/night , Pr ⁽²⁾ _day/night , and Pr ^(m) _day/night to obtain d _e (t), d _f1(t), d _f2(t), and d _m1(t), respectively.

Population hatching rate of resting eggs, h(t)

We developed an individual resting egg hatching rate function dependent on the temperature, T, by taking Yurista’s (1992) hatching rate of 0.014 per day at 10.6°C, a critical temperature of 4°C for hatching to initiate, and the temperature dependency of Eq. 20:

$$ H(T)=\max(0.00214\left(^\circ \hbox{C\,day}\right)^{-1}\cdot T-0.00858\rm{day}^{-1},0). $$

(24)

To calculate the population hatching rate at time t, h(t), we incorporated the temperature the resting eggs were exposed to at the lake bottom, assuming that they were evenly distributed across the lake. With q(ϕ) as the density of Harp Lake’s surface area that has depth ϕ, we can write:

$$ \underbrace{h(t)}_{\substack{{\rm population\,hatching\,rate}}}=\underset{0}{\overset{\phi_{max}}{\int}}\underbrace{H\left(T(t,\phi)\right)}_{\substack{{\rm individual\,hatching\,rate\,at\,depth\, \phi}}}\cdot {q(\phi)}d\phi. $$

(25)

We approximated the integral by summing over a range of depth classes (Fig. 9), with temperatures corresponding to the temperature at the midpoint of each depth class.

Fig. 9

Bathymetric information for Harp Lake, Ontario, used to calculate the hatching rate of resting eggs h(t). The proportion of each depth class is plotted against the midpoint depth of the class

Clutch size, c(t)

Parthenogenetic clutch size, c(t), the number of asexual embryos per brood, has been observed to change over the course of a year (Pothoven et al. 2001; Yan et al. 2001; Young 2008). Where data were available for the early season, a rapid decline has been observed from clutch sizes as high as 9 to averages between 2 and 4 later in summer (Young 2008). Since we do not know of any functional form that has been suggested to describe the clutch size as a function of time or temperature, we fit a generalized linear model to relate clutch size data to time and effective temperature experienced by third instar females, T ⁽³⁾ _e (t). To calculate T ⁽³⁾ _e (t), we evaluated Eq. 19 with g(T(t, x)) = T(t, x) and the vertical distribution for third instar females (Pr ⁽¹⁾ _day/night ).

We pooled parthenogenetically reproducing females from all years whose clutch size had been scored (n = 1,047) and fit a zero-truncated Poisson generalized linear model with a log link (as described in Zuur et al. 2009) using the VGAM package (Yee 2008) in R version 2.10.1 (R Development Core Team 2009) with day of the year and effective temperature experienced by third instar females as factors. Based on Akaike’s information criterion (AIC), we chose the model that includes both temperature and day, as well as their interaction:

$$ c(t)= \exp(-7.83 + 0.581 \cdot T_e^{(3)}(t) + 0.0468 \cdot t - 0.00308 \cdot T_e^{(3)}(t) \cdot t). $$

(26)

Proportion of females that are reproducing sexually, σ(t)

To obtain a relationship between the proportion of females that are reproducing sexually and the predictors time and effective temperature, we fit a binomial generalized linear model with a logit link (Zuur et al. 2009) to data on the reproduction mode extracted from the Harp Lake database for the years 1994–2004 (n=5,064). The best model according to AIC contains temperature, time, and their interaction as factors and is given by:

$$ \sigma(t)=\frac{1}{1+\exp(-60.9 + 0.200 \cdot t + 4.43 \cdot T_e^{(3)}(t) - 0.0146 \cdot t \cdot T_e^{(3)}(t))}. $$

(27)