The persistence of populations facing climate shifts and harvest

Many species are expected to shift their geographic distribution as climates change, and yet climate change is only one of a suite of stressors that species face. Species that might, in theory, be able to shift rapidly enough to keep up with climate velocity (the rate and direction that isotherms move across the landscape) may not in actuality be able to do so when facing the cumulative impacts of multiple stressors. Despite empirical reports of substantial interactions between climate change and other stressors, we often lack a mechanistic understanding of these interactions. Here, we developed and analyzed a spatial population dynamics model to explore the cumulative impacts of climate with another dominant stressor in the ocean and on land: harvest. We found that critical rates of climate velocity and harvest depend on the growth rate and dispersal kernel of the population, as well as the magnitude of the other stressor. This allowed us to identify conditions under which harvesting and climate velocity can together drive populations extinct even when neither stressor would do so in isolation. Except in these extreme cases, we also found that the interaction between the declines in biomass caused by climate velocity and harvest is approximately additive. Finally, we have shown that threshold harvest rules can be effective management tools to mitigate the interaction between the two stressors, while protected areas can either help or hinder, depending on how harvesters reallocate their effort. We also have parameterized the model for black rockfish (Sebastes melanops) to demonstrate the model’s broad applicability.


INTRODUCTION
There are many stressors that can disturb an ecosystem. Ecologists have been working for decades to quantify the consequences of individual perturbations (Wilcove et al. 1998) and to measure the effects of multiple stressors and the interactions between them (Travis 2003, Crain et al. 2008, Darling and Côté 2008. If disturbances interact synergistically, a perturbation that has little effect when occurring alone may amplify the disturbance caused by a coincident perturbation (Gurevitch et al. 2000, Crain et al. 2008, Darling and Côté 2008, Nye et al. 2013). In the most worrying cases, interactions among multiple stressors could drive a population extinct, even though assessments of the individual impacts would not predict extinction (e.g., Travis 2003, Pelletier et al. 2006. Because disturbances rarely occur in isolation, measuring the effects of multiple disturbances provides a better understanding of likely impacts to an ecosystem (Folt et al. 1999, Doak and Morris 2010, Fordham et al. 2013. Climate change and harvesting, two of the largest anthropogenic impacts for both marine and terrestrial species (Milner-Gulland and Bennet 2003, Sekercioglu et al. 2008), provide an important example of two concurrent ecological disturbances. One effect of climate change is that isotherms-contour lines connecting places with the same temperaturemove across a landscape with a rate and direction referred to as climate velocity (Loarie et al. 2009, Burrows et al. 2011. Marine and terrestrial population distributions shift in response to climate change (Perry et al. 2005, Chen et al. 2011, and there is evidence that climate velocities can successfully explain these shifts (Pinsky et al. 2013).
Many of these shifting species are also subject to harvesting or fishing (Wilcove et al. 1998, Sala et al. 2000, Worm et al. 2009), so interactions between the two stressors are possible. For example, empirical data suggest that Atlantic croaker populations move poleward with warming temperatures, but do so less when heavily fished (Hare et al. 2010); several terrestrial species follow warming temperatures more effectively in protected areas than in unprotected areas (Thomas et al. 2012); and a number of studies concluded that harvest increases the sensitivity of populations to climate variability (Anderson et al. 2008, Planque et al. 2010, Shelton and Mangel 2011. While not specifically addressing range shifts and harvest together, there have been experimental indications of synergistic interactions between warming temperatures and harvesting (Mora et al. 2007). Taken together, this work underscores the importance of understanding in greater mechanistic detail how climate velocity and harvesting interact. Models provide a useful tool for building our intuition about this interaction.
A common approach to modeling climate impacts has been to use bioclimatic-envelope models (also known as species distribution models). These statistical models typically correlate presence-absence or abundance data with biophysical characteristics to predict how spe-cies' ranges will differ under climate change (Guisan and Zimmermann 2000, Guisan and Thuiller 2005, Elith et al. 2006. Despite these models' widespread adoption, many authors have criticized bioclimatic-envelope models as oversimplified because they lack dispersal, reproduction, species interactions, and other processes important for population dynamics (Kearney and Porter 2009, Robinson et al. 2011, Zarnetske et al. 2012. Recent work on range shifts has addressed some of these gaps by explicitly including dispersal and reproduction in models of species distributions under climate change (Berestycki et al. 2009, Zhou andKot 2011). In these models, the region in which a population can survive (e.g., the region of suitable temperatures) is shifting in space, and a population can only survive if it disperses to and grows in newly suitable habitat at a sufficient rate. Related models have been applied to study population persistence in advective environments (Byers and Pringle 2006). However, even these more mechanistic models only address one disturbance: climatedriven range shifts.
Here, we focus on a relatively simple ecological model that captures the dominant processes (reproduction, dispersal, and population growth) underlying climate-driven range shifts and population responses to harvesting pressure. We built this model originally for marine species, but because of its mathematical generality, it could also apply to any species with distinct growth and dispersal stages (e.g., plants and many insects). We identify the harvesting rate and climate velocity that drive populations extinct, investigate how the critical rate of one stressor depends on the other, and analyze the declines in biomass caused by each stressor. We also examine two different types of management strategies-threshold harvesting rules and protected areas-to determine how these management strategies affect population persistence and biomass. We chose to model protected areas because they are often recommended for conservation of biodiversity and improved yield from harvest (Pimm et al. 2001, Gaines et al. 2010b, Watson et al. 2011, and previous work has suggested protected areas can be a key form of climate insurance that provides stepping stones to help species keep up with a changing v www.esajournals.org environment (Hannah et al. 2007, Thomas et al. 2012. Finally, we demonstrate how to apply this model by using parameters describing black rockfish (Sebastes melanops) in California (Gaines et al. 2010a, White et al. 2010).

The model
We model the dynamics of populations along a one-dimensional line of longitude, similar to Zhou and Kot (2011). Individuals in the population can only reproduce within a defined segment of this one-dimensional coastline (hereafter simply ''patch''), which represents the range of thermally suitable conditions for the population. The patch shifts at a fixed rate towards the poles, and offspring disperse away from their parents according to a dispersal kernel. In its basic form, harvest removes a constant fraction of the local population density from each point along the coastline.
The above verbal description is represented well by integrodifference models, which have been used extensively for spatial population dynamics problems with discrete time (e.g., discrete growth and dispersal stages) and continuous space (Kot and Schaffer 1986, Van Kirk and Lewis 1997, Zhou and Kot 2011. More specifically, if n t (x) is the number of individuals settling after dispersal at position x and time t, then the number of individuals in the next generation is given by where f(n) is a recruitment function describing the number of juveniles that settle and survive to adulthood given that the juvenile population is of size n, g(n) is a function describing the number of adults that remain after harvesting given local density n. R 0 is the intrinsic growth rate of the population (i.e., number of offspring per adult), and k(x À y) is a dispersal kernel giving the probability of an offspring traveling from position y to position x. Reproduction only occurs within the suitable patch of length L, which shifts across space at a climate velocity c in units of distance per generation. In other words, if t is the number of generations that has passed, the center of the patch will be at location ct, and the upper and lower bounds of the patch will be found at ct þ L/2 and ct -L/2, respectively. Initially, we use g(n) ¼ n -hn as our function for those surviving harvesting, where h is the proportion of the population harvested. This assumes that harvest removes a constant fraction from each location x, as might be expected from an even distribution of harvesters across space. We used a Beverton-Holt stock-recruitment function to describe the settlement and survival of offspring f(n) while accounting for density dependent competition and mortality: As before, R 0 is the intrinsic growth rate, while K is the carrying capacity at a given point in space, which we assume to be constant (see Table  1 for a full description of parameters and functions). Since f(K ) ¼ K/R 0 , if n ¼ K, there will be K/R 0 surviving offspring, and when they reproduce at rate R 0 the population will remain at carrying capacity. As shown in the Appendix, the precise forms of g(n) and f(n) are not important to the persistence of the population. Persistence depends only on g 0 (0) and f 0 (0). The full functional forms do matter, however, for equilibrium biomass.
Analyzing this kind of model becomes easier if the dispersal kernel is separable into its dependence on sources and destinations of larvae, that is, if there are functions a i and b i such that k(x À y) ¼ P ' i¼1 a i ðxÞb i ðy) (see the Appendix for further details). In the analyses presented below, we used a separable Gaussian kernel (Latore et al. 1998) given by To derive analytical expressions for the critical rates of harvesting and climate velocity, we approximate the kernel to its first-order terms, as described in the Appendix. To examine the sensitivity of the model to the shape of the kernel, we also analyze a sinusoidal kernel (see the Appendix).
At demographic equilibrium, the population will move in a traveling wave, where the population density at a given point in space will v www.esajournals.org change, but the density at a location relative to the shifting patch will not (Zhou and Kot 2011). The traveling wave n* must satisfy wherex;ȳ 2 ½À L 2 ; L 2 describes the position within the patch. For a separable kernel, the equilibrium traveling pulse n*(x) must satisfy where the m i satisfy the equations m j a j ðyÞ dy ð6Þ (Latore et al. 1998). We show the derivation of these equations in the Appendix. While there are certainly interesting transient dynamics as the population reaches its equilibrium traveling wave, we focus on equilibrium biomass to make results from different dispersal kernels, parameters, and methods of analysis directly comparable, without the confounding effects of initial conditions and rates of approach to equilibrium.

Calculating persistence
At low harvesting rates h and low climate velocities c, populations will persist. However, above certain critical values, populations will be driven extinct. When the population is extinct, the system is in its trivial equilibrium: n*(x) ¼ 0 for all x 2 ½À L 2 ; L 2 , which satisfies Eq. 4. If a population is to persist, it must be able to avoid extinction and grow even when small (Zhou and Kot 2011). Population persistence is therefore equivalent to the trivial traveling pulse being an unstable equilibrium, where the introduction of a small population will grow rather than return to extinction. The critical parameters h* and c* are defined as the parameters that make the trivial pulse unstable. See the Appendix for further details of this analytical calculation.
Regardless of the functional form of the recruitment function f, the only property that determines whether or not a population can persist is f 0 (0), i.e., how quickly recruitment increases when the population size is near (but above) 0. For us, this number is 1, and any recruitment function with the same value will give the same results with respect to persistence. In addition to this property, the population's ability to persist depends on properties of the population itself (the shape of the dispersal kernel, and the expected distance a larva disperses hdi), properties of the environment (the length of the viable patch L and how quickly the environment shifts c), and the harvesting rate h. For a Gaussian kernel, the critical rates c* and h* are those values of c and h such that We derive a similar expression for a sinusoidal kernel in the Appendix. We realize that this Table 1. Parameters and functions used in the text.
Variable Definition n t (x) density of individuals at position x at time t n*(x) density of individuals at equilibrium at positionx relative to the patch k(x -y) dispersal kernel, the probability of offspring traveling from position y to position x hdi expected distance traveled by an offspring f(n) recruitment function, the number of offspring produced by a population of size n R 0 intrinsic growth rate of the population at low abundance K carrying capacity g(n) harvest function, the number of adults remaining after a population of size n has been harvested h proportion of adults harvested, when g(n) ¼ (1 -h)n L patch length c climate velocity in units of distance per generation v www.esajournals.org formula is not straightforward to interpret. For both Gaussian and sinusoidal kernels, however, we can approximate the critical harvesting proportion by a function that looks like where p is a decreasing function of the length of the viable patch and the intrinsic growth rate, and q describes how h* increases with patch length (L) and varies with expected dispersal distance and climate velocity (see the Appendix for details).

Calculating the interaction of climate velocity and harvest
We identify interactions between climate velocity and harvest in two ways. The first and simplest way is to see if there is an interaction between the critical rate of one stressor and the magnitude of the other stressor. We identify such an interaction if h* depends on c, or if c* depends on h. If this type of interaction exists, determining the critical level of one stressor requires knowing the severity of the second. Before the stressors are extreme enough to drive the population extinct, however, they will cause it to decrease in size. The second way of identifying interactions is to compare how the two stressors affect population biomass individually and jointly. In order to do measure these effects, we find the total biomass of the population when it reaches an equilibrium traveling pulse and compare this equilibrium biomass in the presence and absence of climate shift, harvesting, or both. Eqs. 5 and 6 allow us to numerically find the total biomass in the equilibrium traveling pulse under each of these conditions.
We use B 0 to denote the equilibrium biomass without either stressor, B h the equilibrium biomass with harvesting but with climate velocity equal to zero, B c the equilibrium biomass with climate velocity greater than zero but no harvesting, and B hc the equilibrium biomass with both stressors. For each stressor or combination of stressors, we calculate the decline in biomass caused by stressor s as Based upon this definition, there are three kinds of interaction types that can be defined. If the interaction is additive, then the cumulative response to both stressors together would be E hc ¼ E h þ E c . If the stressors instead interact synergistically, then E hc . E h þ E c . In contrast, if the stressors interact antagonistically, then E hc , E h þ E c . We can therefore quantify the interaction as where positive S indicates synergy, negative S indicates antagonism, and S of zero indicates an additive interaction. This is a common way to measure the interaction among stressors, though an alternative approaches would be to use the ratio of affected to unaffected biomass as a measure of effect size (multiplicative model) or to consider the effect of the single worst stressor (simple comparative effects model; Folt et al. 1999, Crain et al. 2008). The additive model we use here is the most conservative when quantifying negative effects, meaning that it is less likely to identify synergistic interactions (Folt et al. 1999, Crain et al. 2008).

Management strategies
We use simulations to implement two management strategies (threshold harvesting rules and protected areas) that make our basic integrodifference model analytically intractable. We also take advantage of the increased flexibility of simulations over mathematical analysis to use the Laplace dispersal kernel, kðx À yÞ ¼ 1 2 be ÀbjxÀyj ; a commonly used model of marine larval dispersal (Botsford et al. 2001) that is not amenable to the analytical methods we use above. This allows us to show that our results are not qualitatively dependent on our choice of dispersal kernel.
Under threshold harvesting rules, harvesting pressure is no longer implemented as a proportional removal from the population. Instead, we evaluate the abundance at each point in space to determine how much harvesting should occur. If the population abundance is below the designated threshold, no harvesting occurs. If the population exceeds the threshold, then all the 'surplus' individuals are available to be harvested. This approach is an extreme version of the harvest control rules proposed for many existing fisheries (Froese et al. 2011).
In addition, we introduce networks of protected areas into our simulations by designating v www.esajournals.org segments of space where the harvesting rate is equal to zero. Protected areas, particularly in the ocean, are typically designed to meet either harvest management or conservation goals (Agardy 1994, Holland and Brazee 1996, Gaines et al. 2010a, and their spacing and size differ according to which goal is being pursued. Harvest-oriented protected areas are often designed such that they maximize adult spillover into harvestable areas by creating many small, closely spaced reserves (Hastings and Botsford 2003, Gaylord et al. 2005, Gaines et al. 2010a). To mimic this management scheme, we implemented protected areas with a length one third of the average dispersal distance and an inter-reserve spacing two thirds of the average dispersal distance. Conservation-oriented protected areas seek to protect entire ecosystems and reduce adult spillover by creating fewer, larger protected areas (Toonen et al. 2013). To mimic this scheme, we implement protected areas with a length four times the average dispersal distance and an interreserve spacing eight times the average dispersal distance between them . In both harvest-oriented and conservation-oriented protected area networks, one third of the coastline is protected. With protected areas present, we test two ways in which harvesting pressure could respond to reserves: either total harvesting is reduced to two thirds of what it would be without reserves (i.e., harvest effort in reserves is eliminated), or harvesting is shifted to available, unprotected habitat such that total harvesting pressure remains constant (i.e., harvest effort is displaced).
For every simulation, we seed the model with 50 individuals at a single location and iterate for 2000 generations to reach equilibrium without harvesting or climate shift (more than sufficient based on initial tests). We then add harvesting pressure, allow the population to again reach equilibrium (2000 generations), and finally add a changing climate by moving the viable patch with a certain velocity. After 6000 generations we calculate equilibrium biomass as the mean biomass of 2000 additional generations. If population abundance declines below 0.001, the population is considered extinct. Implementing protected areas makes the population abundance cycle, but averaging over 2000 generations is sufficient to erase the effects of periodicity in our results. For most systems, these long timespans are not biologically realistic. However, they ensure that the population reaches its equilibrium traveling wave and that initial conditions do not affect our results. We find qualitatively similar results with shorter simulation times.

Parameters
For our general model investigation, we used the following parameters: R 0 between three and ten, d between 0.1 and 2, K ¼ 100, and L ¼ 1. In this parameterization, hdi is expressed in fractions of the habitable patch width, while c is expressed in fractions of the patch width per generation. In addition, we used life history parameters for black rockfish as an example of how our model can be applied (White et al. 2010). We chose black rockfish because the fish is of both conservation and commercial interest. The parameters for black rockfish in the California Current were as follows: R 0 ¼ 2.86, hdi ¼ 73 km, K ¼ 1, and L ¼ 1000 km (White et al. 2010). We used protected areas with length and spacing representative of the marine reserves put in place by California's Marine Life Protection Act (20 km wide, spaced 76 km apart; Gaines et al. 2010a). For the black rockfish example, we tested climate velocities from zero to 200 km/decade, which was the upper limit observed globally (Burrows et al. 2011). See the Appendix: Table A1 for additional parameter details. While our results depend quantitatively on the parameters of the model, our results are qualitatively robust and we chose a representative set of parameters to analyze.

Persistence with harvesting and climate velocity
We begin by examining the critical rates of harvesting and climate velocity, i.e., those rates sufficient to drive the population extinct. As one might expect, we identify an interaction between the critical rate of one stressor and the magnitude of the other. Specifically, the critical rate of each stressor is lower if a population faces higher intensities of the other stressor (note the negative slope of the lines in Fig. 1). This means that a harvesting rate that is sustainable in the absence of environmental shift (c near zero) may no longer be sustainable if the environment begins v www.esajournals.org to change rapidly (c ) zero). We also found this negative relationship when we parameterized the model for black rockfish (Appendix: Fig. A1).
We also examine the sensitivity of critical rates to growth and dispersal. In our model, it is always the case that increasing the intrinsic growth rate (R 0 ), all else being equal, will increase the critical climate velocity c* and the critical harvesting rate h*, since a population that grows more quickly can recover more effectively from losses caused by these stressors (compare lines with different shading in Fig. 1). However, whether or not dispersing farther is better depends on how quickly the environment is shifting (compare solid and dashed lines in Fig.  1). When the environment is shifting slowly, populations with wider dispersal kernels have a lower critical harvesting rate because dispersing farther results in too many larvae dispersing off the viable patch. When the environment is shifting quickly, on the other hand, populations with wider dispersal kernels can better withstand harvesting because larvae dispersing long distances more effectively colonize the habitat patch that will be viable in the next generation.

Interactions between stressors
We next consider how a population responds to moderate cumulative impacts that are insufficient to drive it extinct. Whenever climate velocity or harvesting pressure exceeds its critical rate, the biomass of the population at equilibrium will be equal to zero (by the definition of the critical rate). Before the stressors reach those thresholds, however, the equilibrium biomass of the population decreases smoothly as either the harvesting pressure or the rate of environmental shift increases ( Fig. 2A). We found the same results when we parameterized the model for black rockfish (Appendix: Fig. A2). The similarity between the shape of the equilibrium biomass surface from our mathematical analysis of an approximation of a Gaussian dispersal kernel ( Fig. 2A) and from our simulations of a Laplace dispersal kernel (Fig. 3A) shows that this result and the following results are robust both to changing our method of analysis and to changing the dispersal kernel.
When we compare the cumulative impacts of the stressors to the sum of each stressor individually, we find low levels of positive synergy between the two stressors (Fig. 2B). The stressors display a synergistic interaction most strongly at high harvest and climate velocity rates, close to where they would drive the population extinct. However, the degree of synergy is low and concentrated in a limited part of parameter space. Throughout much of the range of harvest rates and climate velocities, the interaction between the effects of the stressors is essentially additive. We note that results are robust to changes from a Gaussian to a sinusoidal dispersal kernel.

Alternative management strategies
With harvest thresholds in place, there is a threshold population density below which harvesting is not allowed. Therefore, the population can only be driven extinct by harvesting alone if the threshold is zero, i.e., the whole population is harvested, and otherwise a small population can always escape harvesting. In addition to making it impossible for harvesting to drive a population extinct, the harvest thresholds remove the interaction between the critical climate velocity c* and the harvesting rate h (notice the vertical line dividing positive and zero biomass in Fig. 3B). In this case, the effect of the stressors follows a simple comparative model: the cumulative impacts of the two stressors are equal to the individual effect of the worst stressor. Fig. 1. Lines indicate the critical harvesting rate as a function of climate velocity on the x-axis. The shade of grey corresponds to the growth rate, with darker lines corresponding to higher growth rates. Line style indicates the average dispersal distance. These results are from an approximated Gaussian dispersal kernel with parameters L ¼ 1, K ¼ 100.
v www.esajournals.org If the harvesting rate in unprotected areas is not increased with implementation of the protected areas (i.e., if harvest effort is eliminated instead of displaced), the population withstands combinations of higher climate velocities and higher harvesting rates than without the protected areas. This result applies to either strategy for implementing protected areas (many small versus few large; compare Fig. 3C and D to Fig. 3A). Despite these similarities, there are differences between the strategies of having many small and few large protected areas. At lower climate velocities, small protected areas spaced less than one average dispersal distance apart result in smaller fluctuations of population biomass relative to large spaced protected areas further apart (Appendix: Fig. A3).
If, on the other hand, harvesting effort is reallocated rather than eliminated by the protected areas, the existence of protected areas reduces the critical climate velocity and harvesting rate. In other words, implementation of protected areas in these cases causes extinction of the population at lower climate velocities and harvesting rates than with the case of no protected areas (compare Fig. 3E and F to Fig. 3A, C and D). We find the same qualitative results in our black rockfish parameterization: threshold har-vesting changes the interaction between range shifts and harvesting pressure to a comparative model, and displacing effort outside of protected areas can result in lower population biomass than without protected areas at all (see Appendix: Fig. A4 for details).

DISCUSSION
Climate change and harvest are two of the dominant human impacts on marine species and many terrestrial species, but our understanding of their interaction and joint effects remains limited. By analyzing a general model that incorporates dispersal and reproduction with a set of representative parameters and parameters describing black rockfish, we find an interaction between the critical rate of the each stressor and the magnitude of the other, such that the critical harvesting rate decreases as climate velocity increases and vice versa. In other words, the more quickly the environment shifts, the less harvesting it takes to drive the population extinct. We then find that climate velocity and harvesting interact essentially additively in their effects on biomass for most combinations of stressor levels, with weak synergy only appearing close to population extinction.  Fig. 2A). (B) Equilibrium biomass for simulations with threshold management. For threshold management, the maximum threshold below which no harvesting is allowed is set to be the largest population size observed at a given time step before harvesting. For a less severe threshold, we use a proportion of this maximum threshold, so that a lower proportion gives a lower threshold and allows for more harvesting. We show this proportion on the y-axis. (C) Equilibrium biomass for simulations with many small protected areas with harvesting pressure outside reserves unchanged. (D) Equilibrium biomass for simulations with few large protected areas with harvesting pressure outside reserves unchanged. (E) Equilibrium biomass for simulations with many small protected areas with harvesting pressure reallocated outside reserves. (F) Equilibrium biomass for simulations with few large protected areas with harvesting pressure reallocated outside reserves. These results are from a simulation with a Laplacian dispersal kernel with parameters L ¼ 1, R 0 ¼ 5, K ¼ 100, and hdi ¼ 2.
v www.esajournals.org Our results suggest that particular combinations of harvesting and climate velocity will affect certain species more than others. Species with a higher intrinsic population growth rate (i.e., growth rate at low abundance) and a longer average dispersal distance will better track rapid climate velocities, as compared to species with a low intrinsic population growth rate and short dispersal distances. This finding matches previous expectations: higher growth rates make a population more resistant to the removals from harvesting or the losses associated with tracking climate velocity. It is worth pointing out that a higher population growth rate can be generated by shorter generation times, higher fecundity, or higher survival. Empirical work also suggests that marine fish and invertebrates with faster life histories, as well as terrestrial birds and plants with greater dispersal abilities, shift their distributions more quickly in response to warming (Perry et al. 2005, Angert et al. 2011, Pinsky et al. 2013. While higher reproductive rates improve a population's ability to persist in our model, higher dispersal distances do not necessarily do so. In agreement with related results from Zhou and Kot (2011), we found that at low speeds, a short dispersal distance improved the maximum harvesting rate a population could sustain, while at higher speeds a longer dispersal distance improved the maximum climate velocity under which the population could persist. It appears that climate velocity could selectively favor species with dispersal distances best matched to the rate of shift.
Our finding that the interaction between harvest and climate velocity on biomass is effectively additive would appear to contrast with demonstrations of synergy between harvest and climate in the literature. For example, a number of modeling and empirical studies have found that fishing increases the sensitivity of populations to climate variability (including Anderson et al. 2008, Shelton and Mangel 2011, and a recent review reaches the same conclusion (Planque et al. 2010). Positive feedback loops involving the loss of predators due to fishing have also been identified that amplify climate impacts on prey species (Kirby et al. 2009, Ling et al. 2009, Planque et al. 2010. Similarly, synergy between harvesting and temperature was detected in experimental populations of rotifers (Mora et al. 2007).
A partial explanation for the differences between our model results and the previous evidence for synergy may be that we analyze the ability of populations to keep pace with climate velocity, while many previous studies examined other aspects of changing climate. In the rotifer experiment, for example, populations were subjected to warming temperatures, but organisms were unable to relocate to thermal optima (Mora et al. 2007). In many other fishing and climate studies, the impacts of climate variability on stationary populations have been the focus, rather than cumulative climate change or shifting distributions (Walters and Parma 1996, Anderson et al. 2008, Planque et al. 2010, Shelton and Mangel 2011. Work that does incorporate shifting species distributions typically examines regional or global scenarios for climate change, making it difficult to isolate the effect that different species interactions, climate and harvesting each play (Cheung et al. 2010).
Another explanation for the discrepancy may be that the only effect of harvesting in our model is a reduction in the amount of the adult biomass. In reality, populations often contain a diversity of subpopulations, ages, and genotypes that can buffer them against climate variability and climate change (Schindler et al. 2010). Harvest tends to simplify this diversity within populations, making them more sensitive to climate variability (Mora et al. 2007, Planque et al. 2010). In addition, some synergistic interactions between climate and harvesting identified in previous studies involved the loss of predators and the release of prey (Kirby et al. 2009, Ling et al. 2009), but our model did not include food web dynamics or species interactions and thus was unable to capture these dynamics. Our simple, single-species, non-age-structured model suggests that additive interactions between climate velocity and harvesting constitute a reasonable baseline or ''null'' expectation in the absence of more complicated mechanisms. Future work considering food web processes and genetic, spatial, and age diversity will be important to examine other possible sources of synergistic (or antagonistic) interactions between harvesting and climate velocity.
We also examine whether two frequently v www.esajournals.org recommended management approaches, protected areas and harvest control rules, could help ensure species persistence in the face of multiple stressors, again both for a general set of parameters and for parameters describing black rockfish. Threshold harvesting rules in particular appear to fundamentally alter how the two stressors interact. In particular, the interaction between the critical rates is fundamentally altered: the critical climate velocity no longer depends on harvesting and as long as the climate velocity is below this critical rate, the population size is determined only by the magnitude of harvesting. In our model, thresholds appear to have this effect because they effectively prevent harvesting of the leading edge and allow colonization to occur as if these individuals were moving into un-harvested areas.
While we framed our model as one that describes a population following a shifting climate gradient, it shares many features with a population that is invading new territory. Our results match well with invasion theory, which has shown that populations move into new territory at a rate approximately equal to 2 ffiffiffiffiffiffi ffi R 0 l p , where l is the mean squared displacement of individuals per unit time (Fisher 1937). With a constant harvest rate applied everywhere, the invasion rate would drop to 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ð1 À hÞR 0 l p , whereas the invasion rate would be unaffected if harvesting avoided the leading edge, in accordance with our finding that protecting the low-abundance leading edge from harvesting can mitigate the effect climate shift. Since this elegant early result, theoretical and empirical work in invasion biology has shown that a low growth rate at the leading edge of a moving population, which could, for instance result from an Allee effect caused by the low population density there, can slow down or prevent an invasion (Lewis and Kareiva 1993, Kot et al. 1996, Veit and Lewis 1996, Hastings et al. 2005. It is interesting to note that newly colonized populations, which initially appear at low abundance, are commonly unregulated in fisheries systems (Beddington et al. 2007, Dowling et al. 2008. Whether fisheries and other harvesting activities rapidly exploit newly colonizing species depends in part on the interaction of social, economic, and regulatory factors (Pinsky and Fogarty 2012). Our work highlights the fact that a low (or zero) harvest rate on species that have recently colonized new habitats can be important for helping them keep up with rapid climate velocities.
Previous work has advanced protected areas as a way to help organisms keep pace with shifting climates, as well as to ameliorate anthropogenic disturbances like harvesting and habitat fragmentation (Botsford et al. 2001, Hastings and Botsford 2003, Gaylord et al. 2005, Hannah et al. 2007, Lawler et al. 2010, Watson et al. 2011, Thomas et al. 2012. We find that protected areas can actually make the population more vulnerable to climate change and harvesting pressures than a scenario in which no reserves are present if harvesting pressure is reallocated to unprotected areas. If, on the other hand, harvesting pressure within reserves is removed from the system, our results show that protected areas increase the critical climate velocity and harvest rate of harvested populations. Since reallocation of harvesting effort has the effect of increasing the harvest rate in unprotected areas, this result matches our earlier finding that high harvest pressures at the leading edge of a population can make it more vulnerable to climate velocity. In a theoretical model of an initially small population invading a patchy environment, decreasing the growth rate in the unfavorable patches made it harder for the population to invade (Shigesada et al. 1986, Kinezaki et al. 2003). This agrees with our finding that reallocating harvesting pressure to unprotected areas increases sensitivity to stressors.
We also find that the details of protected-area design affect our results. Few, large protected areas increase population fluctuations at low climate velocities as the population moves through protected and unprotected areas. Many smaller protected areas, on the other hand, maintain a population whose minimum biomass is higher, which could potentially provide a buffer against extinction caused by stochastic events. This occurs because harvest drives the population to lower levels while between protected areas. The larger those gaps are, the more diminished the population will be during its transit.
Whether many small or fewer large protected areas is better depends on many factors and is often species-or system-specific (Gaines et al. 2010b, McCarthy et al. 2011. Halpern (2003) found in a meta-analysis of empirical studies of MPAs that the benefits from implementing an MPA did not depend strongly on its size, though Claudet et al. (2008) found that fish density increased with reserve size. Using a theoretical model, Neubert (2003) found that the optimal MPA spacing to maximize harvesting yield depended on the length of the region in which the population could survive: as the length increased, more and smaller MPAs became preferable. Increasing the length of the viable region is equivalent to increasing the size of the habitable patch in our model, so that his results are similar to our findings. On the other hand, McLeod et al. (2008) argued that having fewer larger MPAs should increase an ecosystem's resilience to climate change by protecting selfpersistent populations. Similarly, Moffitt et al. (2011) used a theoretical model to compare MPAs that were 10 km long spaced 50 km apart to MPAs that were 20 km long spaced 100 km apart and found that the larger more widely spaced MPAs would support the persistence of a greater number of species types. However, none of these studies considered a population moving across a network of MPAs. By considering how a population will track a moving isotherm, we contribute to this body of work by showing that small gaps between protected areas may help species keep up with climate velocities in the face of harvest and that considering a shifting climate is important for making recommendations about MPA spacing.
The advantage of a simple model like ours is that it is potentially general enough to apply to a wide range of species. Our discrete-time, continuous-space model captures the processes important to species with distinct growth and dispersal stages, which includes most marine organisms, plants, and many insects. Our approach does not capture all the complexities of real populations or of harvesting dynamics, however. For example, we do not include the potential for negative per capita growth at low densities, often called Allee or depensation effects. Allee effects can make it more difficult for a population to invade a new environment (Lewis and Kareiva 1993, Kot et al. 1996, Veit and Lewis 1996, Hastings et al. 2005. We would also expect that populations with Allee effects would be more sensitive to the combined effects of harvest and climate velocity than our model initially suggests. We also did not include age structure or other aspects of subpopulation diversity (e.g., spatial or genetic) in our model. As described above, these forms of diversity have been important for studying the joint effects of harvesting and climate variability (Planque et al. 2010, and they will likely be important for understanding climate velocity impacts as well. Besides these species-specific extensions, our modeling framework could be extended to consider species interactions, such as between predator and prey (Gilman et al. 2010). There are some rules of thumb to predict how multiple stressors will affect multispecies systems. For example, ecosystems that contain at least some species tolerant to a wide range of stressors (positive species co-tolerance) can more effectively maintain functioning in the face of climate change (Vinebrooke et al. 2004). End-to-end simulation models, which incorporate physical environmental drivers and describe the dynamics of species at multiple trophic levels (e.g., Travers-Trolet et al. 2014), are increasingly popular as a framework for modeling multispecies systems (Fulton 2010). Because our model is not specific to a particular region or set of species, it can be used as a complement to these larger simulation studies. Hollowed et al. (2000) recommend caution in building overly detailed models because determining model sensitivity and understanding (sometimes hidden) assumptions becomes difficult.
A final important extension would be to represent harvesting dynamics more realistically. Our results show that the success of protected areas is diminished if harvest is reallocated to unprotected areas. Previous studies have also found that the details of how effort is reallocated can change the predicted effects on population dynamics (Kellner et al. 2007). Whether or not harvesting pressure is reallocated, fishermen often focus their efforts at the boundaries between protected and unprotected areas, where the spillover from the MPAs is likely to be highest. There are circumstances under which fishing the line can lead to comparable biomass and overall catch relative to uniform harvesting pressure in unprotected areas (Kellner et al. 2007). However, in our model, fishing the line would reduce the low-abundance leading edge as it moves into an unprotected area and we therefore expect that it would make it more difficult for a population to persist. To the extent that harvester behavior has been considered in fisheries, there is considerable uncertainty in how vessels allocate effort over space and respond to changes in environmental and regulatory conditions (Wilen et al. 2002, Fulton et al. 2011, Pinsky and Fogarty 2012. Harvest behaviors are rarely integrated into modeling efforts, and an important next step will be integrated assessments of social-ecological systems.
Using a simple, mechanistic model like the one we present here helps to build intuition about the conditions under which species can survive the cumulative impacts of climate and harvesting. This work highlights the importance of considering stressors in combination, as outcomes deviate from what we would predict in isolation. It also shows the importance of management choices, as the location of harvest greatly affects the interaction between harvesting and climate. While fisheries management strategies only change harvesting practices and do not directly address climate change, understanding how regulations can affect interactions between harvesting and range shifts can help to improve harvesting rules and the development of protected areas. Our results offer encouraging evidence that management practices can help protect marine populations from the cumulative impacts of harvesting and climate change, particularly if the location of harvesting can be controlled.