Evolution of movement rate increases the effectiveness of marine reserves for the conservation of pelagic fishes

Abstract Current debates about the efficacy of no‐take marine reserves (MR) in protecting large pelagic fish such as tuna and sharks have usually not considered the evolutionary dimension of this issue, which emerges because the propensity to swim away from a given place, like any other biological trait, will probably vary in a heritable fashion among individuals. Here, based on spatially explicit simulations, we investigated whether selection to remain in MRs to avoid higher fishing mortality can lead to the evolution of more philopatric fish. Our simulations, which covered a range of life histories among tuna species (skipjack tuna vs. Atlantic bluefin tuna) and shark species (great white sharks vs. spiny dogfish), suggested that MRs were most effective at maintaining viable population sizes when movement distances were lowest. Decreased movement rate evolved following the establishment of marine reserves, and this evolution occurred more rapidly with higher fishing pressure. Evolutionary reductions in movement rate led to increases in within‐reserve population sizes over the course of the 50 years following MR establishment, although this varied among life histories, with skipjack responding fastest and great white sharks slowest. Our results suggest the evolution of decreased movement can augment the efficacy of marine reserves, especially for species, such as skipjack tuna, with relatively short generation times. Even when movement rates did not evolve substantially over 50 years (e.g., given long generation times or little heritable variation), marine reserves were an effective tool for the conservation of fish populations when mean movement rates were low or MRs were large.


| 445
MM Met al in a single year (Bonfil et al., 2005;Sibert & Hampton, 2003), and reserves may afford protection only for a limited portion of a species' life cycle, especially for migratory species (McAllister, Barnett, Lyle, & Semmens, 2015). Yet, even within such species, variation exists, with many individuals displacing less than 50 nautical miles (NM) between tagging and recapture, and some displacing more than 1,000 NM (Sibert & Hampton, 2003). Whereas the effect of marine reserves on evolution of size at maturation has received considerable attention (Baskett, Levin, Gaines, & Dushoff, 2005;Dunlop, Baskett, Heino, & Dieckmann, 2009;Miethe, Pitchford, & Dytham, 2009), a previously neglected effect of marine reserves (but see Miethe et al., 2009) is that more mobile individuals are more likely to leave a reserve and be fished, while less mobile (or more philopatric) individuals can remain within (or return to) reserves and survive. As a result, over time, the number of fish exhibiting low movement rates should rise. An evolutionary change in movement rate may, in turn, result in an increased density of fish in the marine reserves. There is a growing recognition that fisheries management and conservation should account for evolutionary effects (Fraser, 2013;Palkovacs, 2011). By neglecting the potential evolution of reduced movement rate following reserve establishment, managers may underestimate the efficacy of marine reserves.
For a trait to evolve, it must be heritable. Estimates of heritability for movement rate for pelagic species such as tuna and sharks do not, however, currently exist. There is a large amount of variation among individuals in movement rate (Sibert & Hampton, 2003), and it is unlikely that none of the variation in movement rate is heritable. Furthermore, there is an extensive literature demonstrating the heritability (and genetic basis) of movement rate across taxa. For example, in one of the best-studied cases, a large amount of the variation in movement rate of larval and adult fruit flies depends on the genotype at the for ("foraging") locus (Edelsparre, Vesterberg, Lim, Anwari, & Fitzpatrick, 2014;Sokolowski, 1980). A heritable genetic basis for many behaviors related to movement rate has also been identified in a variety of fish species (Alos, Palmer, & Arlinghaus, 2012;Booth, Hairston, & Flecker, 2013;Breau, Weir, & Grant, 2007;Conrad, Weinersmith, Brodin, Saltz, & Sih, 2011;Coombs & Rodriguez, 2007;Diamond et al., 2007;Laine, Primmer, Herczeg, Merila, & Shikano, 2012;Morrissey & Ferguson, 2011;Roy, Roy, Grant, & Bergeron, 2013;Sneddon, 2003;Steingrimsson & Grant, 2003;Wilson & Godin, 2009;Wright, Nakamichi, Krause, & Butlin, 2006). In general, studies of mobile animals usually find that some individuals move longer distances than the majority of individuals and that variability in movement is the direct result of heritable differences in movement rate, which in turn has fitness consequences (Fraser, Gilliam, Daley, Le, & Skalski, 2001).
In this study, our objective was to determine the extent to which the evolution of movement rate improves the efficacy of large marine reserves in protecting pelagic fish populations. Our hypothesis was that heritable variation among individuals in mean movement distance, combined with higher mortality among individuals that move outside of marine reserves (i.e., natural selection), should promote the evolution of low movement rates, and thereby result in more individuals staying within the boundaries of marine reserves. We simulated individual-level processes (births, deaths, and movement) in a spatially explicit (i.e., grid-based) model combining population dynamics with evolutionary dynamics. We chose four distinct parameter sets, informed by empirical data, to represent the range of life histories (e.g., short life span and small size vs. longer life span and large size) among tuna species (skipjack tuna vs. Atlantic bluefin tuna) and shark species (great white sharks vs. spiny dogfish). We predicted that the evolution of decreased movement rate following the establishment of marine reserves would lead to higher within-reserve population size relative to situations with no heritable variation in movement rate. We also predicted that higher fisheries-induced mortality (e.g., overfishing) would lead to faster evolution of reduced movement.

| METHODS
Our model builds upon the individual-based model of marine reserves developed by Miethe et al. (2009), which considered the evolution of movement rate between two patches and the evolution of size at maturity for a haploid, asexual organism. We also considered the evolution of movement rate, but we did so in a spatially explicit framework.
Furthermore, we modeled fish more realistically as diploid, sexually reproducing organisms. We tracked individual-level processes across multiple patches arranged in a grid, wherein each patch (or grid cell) is analogous to a square area of ocean of 100 NM in length ( Figure 1).
The size and arrangement of grids were adjusted based on speciesspecific data. For example, because the "tuna belt" in the Pacific Ocean consists of about 8,000,000 NM 2 stretching about 2,000 NM from latitude 20°N to 20°S and about 4,000 NM from 150°E to 140°W, we used a 20 × 40 grid to represent the Pacific Ocean tuna belt. Typically, a single marine reserve encompassing the entire Exclusive Economic Zone (EEZ) of an oceanic island consists (roughly) of a circle of 400 NM in diameter, which we modeled as a square reserve of equal length (16 cells in a 4 × 4 configuration).
We tracked the dynamics at a single two-allele gene causing variation among fish in their movement distances. Each individual's movement was determined by its diploid genotype: AA, Aa, or aa, with a 10% initial frequency of allele a (the low-movement allele). Note that these genotypes did not cause a switch between migratory and nonmigratory behavior, but rather affected the net distance travelled, which might be due to changes in angles of movement, rather than total distance traveled. Hence, we do not consider that life history changes necessarily accompany these shifts in net distance traveled, although further work could examine pleiotropic effects or subsequent evolution of life history alongside movement rate evolution. Each genotype conferred a mean (or expected) net movement distance, μ i (Table 1), with the actual net movement distance of an individual in NM drawn randomly from a negative binomial distribution with variance equal to μ i + μ i 2 /2. The resulting distribution of movement distances (based on the negative binomial) is similar to a gamma distribution with rounding, which is appropriate for our purposes because (1) the range is appropriate (including fish that remain in the same cell and fish that can move a large distance), and (2) the distribution closely resembles empirically observed distributions of movement distances (Sibert & Hampton, 2003). Additive gene action was assumed (dominance coefficient: dc = 1/2). Movement direction, θ, was determined by a random draw from a uniform distribution ranging from 0 to 2π (0 representing east and π/2 representing south). We then calculated the net distance moved in the longitudinal direction (cosθ × distance) and latitudinal direction (sinθ × distance). If the distance in either direction brought an individual to the edge of the grid, any additional distance was reflected in the opposite direction. Each individual had a nonzero probability of remaining within the same cell (given by the probability of drawing a net movement distance smaller than the distance from the center of a cell to its edge), and otherwise it moved to other cells probabilistically, according to the direction and angle of movement.
Modeling movement probabilistically across a grid greatly increased the speed of the simulations; given the relatively large number of grid cells, the loss of accuracy is not expected to be substantial, relative to explicitly tracking movement of each individual in continuous space (as noted in Berec, 2002;Mooij & DeAngelis, 2003). Sex (male or female) was assigned randomly to each individual at birth. At the time of reproduction, mating occurred only if there was at least one mature male within the grid cell of a female, with all fertilized females producing the same mean number of eggs (determined by the fecundity parameter in Table 1). Within a cell, the number of eggs carrying the A allele was drawn randomly from a Poisson distribution with a mean of fecundity times the number of AA females + fecundity times the number of Aa females/2, and the number of eggs carrying the a allele was similarly based on the number of aa and Aa females. We assumed no sperm limitation within cells containing a male, and the allele frequency among sperm was determined by the allele frequency among males in the cell. Allowing more individual variance in fecundity (e.g., by using a negative binomial distribution) would increase genetic drift (i.e., decrease effective population size) but should not have any other effect on our results. The number of zygotes (and, subsequently, hatchlings and fry) of each genotype was determined by randomized binomial trials, given the number of eggs carrying each allele and the allele frequency among sperm. The low F I G U R E 1 Simulation grids and reserve designs (a-c) used for simulations involving tuna and shark species. Arrows indicate that migration occurs between feeding and spawning grounds for bluefin tuna or, for great white sharks, between a local population (South Africa-on the left of each panel) and a nonlocal population (Western Australia-on the right). Dark gray areas indicate reserve locations. Each grid cell represents a 100 × 100 NM patch of ocean fecundity values used in our simulations (Table 1) relative to the actual number of eggs produced by tuna (e.g., a large bluefin tuna female can produce hundreds of millions of eggs in a single spawning bout) were meant to reflect the high mortality between the egg and fry stage.
Ovoviviparous shark species, on the other hand, typically produce few offspring, and shark pups likely have a higher survival rate than tuna fry (Compagno, 1984(Compagno, , 2002. There were three age classes in our model: fry, juveniles, and adults. Fry either died or survived and were then recruited to the juvenile age class in the year of their birth (i.e., all fry were age 0). We assumed no larval dispersal, meaning that fry did not move outside of the cell in which they were spawned. Fry survival was negatively density dependent based on the number of fry in the grid cell: where s b is the proportion of fry that would survive if there were no density dependence, and dd determines the magnitude of density de- The mortality value for each species was selected from within the range of values reported on fishbase.org. b The effect of density on survival (i.e., the value of the dd parameter) was initially set to the value used in Miethe et al. (2009) and then adjusted (for the shark species), after running several pilot simulations, to a value that resulted in population persistence without any fishing pressure. c See Table 2 for the source of these age of maturity values.
of being caught was f i,j , and the number of trials was equal to the number of individuals in the grid cell after natural mortality. Fishing pressure was assumed to be sufficiently high that adults are effectively released from density-dependent competition, and, aside from fry survival, we did not include density-dependent processes in our model. A direction for future work might be to relax this assumption (e.g., Gårdmark, Jonzén, & Mangel, 2006), especially in cases of highly successful MPAs where local densities can rebound and shift the life history stage at which density dependence is strongest.
After reserve establishment, the grid cells that were designated as reserves experienced no fishing mortality ("no-catch" reserves). Total fisheries pressure was, however, kept constant (i.e., fishing effort was displaced) by considering only nonreserve patches in the calculation of mean.patch.pop, and by increasing mean fishing mortality outside of the reserves in proportion to the area of the reserves (i.e., if a fraction x of cells were designated as reserves, the adjusted mean fishing mortality in the remaining cells was set to fished/(1 − x). Fishing mortality continued to track fish abundance after reserve establishment.
Consequently, any spillover of fish from reserves to neighboring cells led to higher abundance and hence higher fishing mortality adjacent to reserves.
The simulations proceeded in yearly time increments with the following order of events: spawning-natural mortality-recruitment-fishing-movement. The initial frequency of the a allele (init.a) was set to 0.1 in all grid cells, and all grid cells started with the same initial population size (1,000 for skipjack and dogfish simulations, 500 for bluefin and great white shark simulations). The heritability at the beginning of each simulation (h 2 = V A /V P ) was estimated by randomly drawing movement distances in proportion to the genotype frequencies to obtain an estimate of total phenotypic variation (V P ) and using the onelocus model (Falconer & MacKay, 1996) to calculate additive genetic init.a, Table 1), and α is the average allelic effect: α = (μ AA − μ aa )/2 for the additive case considered here (dc = 1/2). We note that alternative genetic models, such as quantitative genetic models that assume constant heritability, may be less appropriate, especially for evolution in small populations, where evolutionary responses might be strongly constrained/limited to the dynamics of a few loci. Population size was allowed to stabilize across all cells for 50 years prior to the commencement of fishing. We then simulated 50 years of fishing prior to the establishment of reserves. After reserve establishment at year 100, we ran the simulations for 50 years (except where noted), which is a time frame of relevance for resource managers (at least, more so than longer evolutionary time frames). The efficacy of reserves was evaluated by comparing mean within-reserve population density (i.e., individuals per cell) at year 150 (i.e., 50 years postestablishment) to population density at year 150 in control simulations without any reserves. We ran ten replicates of each parameter combination. For each species (skipjack tuna, bluefin tuna, dogfish, and great white shark), for both levels of fishing mortality (low or high), we evaluated the effect of reserve number (or reserve size in the case of great white sharks) and movement distance on within-reserve population density and frequency of allele a 50 years after reserve establishment. We used Dunnett's test of multiple comparisons (Dunnett, 1955) to compare every parameter combination (i.e., each unique combination of reserve number or size, μ AA and μ aa ) to the control case (with no reserves). All simulations and analyses were performed using the R programming language (Hothorn, Bretz, & Westfall, 2008; R Core Team 2014). The R code used to run the simulations is provided in Supporting Information.
Given the lack of information in the literature that would inform values of the fry or pup survival parameter (s b ), we based the value for this parameter on two a priori assumptions. First, we assumed a single T A B L E 2 Range of selected life history characteristics (age at maturity, asymptotic length, and growth rate) among tuna and shark species (from fishbase.org; Compagno, 1984;Fromentin & Fonteneau, 2001;Compagno, 2002) Species

| Skipjack tuna
The majority of skipjack tuna catch worldwide occurs in the "tuna belt", which is an area extending approximately 2,000 NM from Wake Island in the north to New Caledonia in the south, and 4,000 NM from Indonesia in the west to the Marquesas in the east. We, therefore, used a 20 × 40 grid for our simulations involving skipjack tuna ( Figure 1). Sibert and Hampton (2003) Sibert and Hampton (2003). Heritability at the beginning of these simulations varied from approximately 0.01 to 0.06, except when μ AA = μ aa (no heritability). For skipjack tuna, we compared the effects of establishing one, two, or four reserves each of size 400 × 400 NM (Figure 1).

| Bluefin tuna
Bluefin tuna and other tuna in the genus Thunnus, such as bigeye tuna (T. obesus) and albacore (T. alalunga), are highly migratory species (FAO 1994  productive waters off the west coast of North America (Bayliff, 1994).
Although total ranges for these species (including transoceanic migrations) are very large, individuals may spend substantial periods of time in much more restricted areas off the coasts of North America or in the western Pacific (Bayliff, 1994;Boustany, Matteson, Castleton, Farwell, & Block, 2010). Therefore, in our simulations for bluefin tuna, we modeled migration and movement as two separate processes. All individuals migrated from the spawning ground to the feeding ground and back prior to maturity (Figure 1). The spawning and feeding grounds were simulated as separate 20 × 20 grids linked by migration ( Figure 1). The rate of movement was allowed to evolve within each area, but not the rate of migration between spawning and feeding grounds, which was determined by the mig and maturity.age parameters (Table 1). After birth, individuals migrated from the spawning to the feeding ground at (on average) maturity.age-mig, then migrated back to the spawning ground after (on average) mig years in the feeding grounds. The average time spent in the feeding and spawning grounds (and, hence, the value of the mig parameter) was drawn from Bayliff (1994). Migration between spawning and feeding grounds occurred in each time step prior to spawning (such that the yearly series of events was: migration-spawning-natural mortality-recruitmentfishing-movement). The earliest an individual could have spawned (if they migrated to the feeding ground in year 0 and returned to the spawning ground in year 1) was at age 2. Adults left the spawning ground in the time step immediately after spawning and then spent (on average) mig years in the feeding grounds before returning once again to the spawning grounds. Movement was modeled as a process that only occurred within the spawning or feeding grounds. We selected values for μ AA of 800, 600, and 400 NM and values of μ aa of 800, 600, 400, and 200 NM (Table 1). We contend that this range of values likely brackets the real annual displacement distances for bluefin tuna within either their feeding or spawning grounds (Bayliff, 1994;Block et al., 2005;Boustany et al., 2010). Heritability at the beginning of these simulations varied from approximately 0.006 to 0.06, except when μ AA = μ aa (no heritability). For bluefin tuna, we compared the effects of establishing one, two, or four reserves, each of size 400 × 400 NM, within the spawning grounds ( Figure 1).

| Spiny dogfish
Given the implications of climate change models for marine biodiver-  (Carlson et al., 2014;McFarlane & King, 2003). Heritability at the beginning of these simulations varied from approximately 0.03 to 0.07, except when μ AA = μ aa (no heritability). For spiny dogfish, we compared the effects of establishing one, two, or four reserves, each of size 400 × 400 NM (Figure 1).

| Great white sharks
Great white sharks are typically a coastal species, with major centers of abundance off the coastal waters of Baja California, Australia-New Zealand, and South Africa (Bonfil et al., 2005;Compagno, 2002).
White sharks are not considered a highly migratory species (FAO 1994), but they are known to make regular transoceanic migrations, perhaps for the purposes of mating (Bonfil et al., 2005). Indeed, white sharks have been observed to make the 20,000 km round trip between South Africa and Australia in less than nine months (Bonfil et al., 2005). Great white sharks also move extensively within their home range. For example, a tagging study revealed that approximately 25% of individuals tagged off the Western Cape of South Africa travelled 2,000 km to KwaZulu-Natal (near the border of Mozambique), and 12.5% made the return trip within a single year (Bonfil et al., 2005). In our simulations for great white sharks, we used a 16 × 3 grid to represent the coastline off South Africa and Mozambique, extending 300 NM from shore ( Figure 1). We biased movement within this coastal region to be along the coast (as opposed to seaward) by converting θ to the root of the function x − cos x × sin x − within the interval 0 < x < 2π (leading to >80% of movements along the coast).
We selected values for μ AA of 800 and 600 NM and values of μ aa of 800, 600, 400, and 200 NM (Table 1), chosen to bracket the likely annual displacement distances for great white sharks within their natal F I G U R E 2 Evolution of movement rate and changes in within-reserve population density for skipjack tuna (a), bluefin tuna (b), dogfish (c), and great white shark (d). Each arrow corresponds to one parameter combination and shows the evolution of movement rate (horizontal axis) and the associated change in population density (vertical axis) 50 years after reserve establishment (at the arrow head) relative to prereserve conditions (circle or diamond symbols at the tail of the arrow). Diamond symbols indicate prereserve conditions for scenarios with no heritable variation in movement distance (and hence no possibility of evolution  (Bonfil et al., 2005). Heritability at the beginning of these simulations varied from approximately 0.006 to 0.06, except when μ AA = μ aa (no heritability). Note that, using a negative binomial distribution for movement with a mean of 800 or 600 NM and a variance of μ + μ 2 /2, approximately 25% of individuals would travel at least 1,000 NM, which is on the order of observed movement distances along the coast of South Africa (Bonfil et al., 2005). We also allowed long distance (transoceanic) migration to another coastal population (e.g., off the Australian coast). All individuals spent an average of two years in their natal population (i.e., South Africa) before migrating, and migrants always returned to their natal population within a year.
For great white sharks, we compared the effects of protecting 8.8%, 16.6%, or 25% of the total grid, placed solely within South Africa's EEZ, corresponding to a reserve that spans 400, 800, or 1,200 NM of coastline, respectively (Figure 1).

| Skipjack tuna
Within-reserve population densities of skipjack tuna were significantly higher than population densities without reserves (Table 3,  Figs S1 and S2). In the scenarios with high fishing mortality, the effect of marine reserves was much more substantial than in equivalent scenarios with low fishing pressure (Figure 2a, Figs S1 and S2). Regardless of fishing mortality, all combinations of reserve number and movement distances resulted in significant increases in within-reserve population density relative to the control case without reserves (Table 3).
Regardless of fishing mortality or number of reserves, marine reserves were more effective when movement distances (μ) were lowest ( Figure 2a, Figs S1 and S2). Within-reserve population densities of skipjack tuna were higher than population densities outside reserves, and marine reserves also lead to increases in population density in nonreserve grid cells adjacent to marine reserves (Figure 3a).
Consistent with our predictions, decreased movement rate evolved following the establishment of marine reserves (Figures 2a   and 3b, Figs S3 and S4). There were significant increases in the frequency of allele a (the low-movement-rate allele) in all cases except when there was no heritable variation (i.e., when μ AA = μ aa ) or with the lowest nonzero heritability and the lowest selection pressure (i.e., when fishing mortality was low, there was one reserve, μ AA = 600 NM, and μ aa = 400 NM; Table 3). Movement evolution was much more extensive with high fishing mortality than with low fishing mortality (Figure 2a, Figs S3 and S4). Greater genetic changes were seen with higher heritable variation (i.e., when μ AA − μ aa was large; Figure 2, Figs S3 and S4). These evolutionary reductions in movement rate (i.e., increases in the frequency of allele a) were associated with substantial increases in within-reserve population sizes over the course of the 50 years following reserve establishment ( Figure 2a). The spatial pattern of allele frequency change indicates that the evolution of decreased movement was concentrated in and around reserves (Figure 3b).

| Bluefin tuna
Within-reserve population densities of bluefin tuna were significantly higher than population densities without reserves (Table 4, Table 4). Consistent with our prediction, however, there was substantial evolution of movement rate when fishing pressure was high (fished = 0.75; Figure 2b, Fig. S8), and greater evolutionary changes were seen with higher heritable variation (i.e., when μ AA − μ aa was large; Table 4; Figure 2b, Fig. S8). There were significant increases in the frequency of allele a (the low-movement-rate allele) in all cases except when there was no heritable variation (i.e., when μ AA = μ aa ) or very low heritable variation (i.e., when μ AA = 800 NM, and μ aa = 600 NM; Table 4). The spatial pattern of allele frequency change indicates that the evolution of decreased movement was concentrated in and around reserves (Figure 4b). Varying the length of time spent in the feeding grounds prior to returning to the spawning grounds (mig = 2 or 6) did not have a qualitative effect on the patterns described above, except that a longer residency in the feeding grounds (mig = 6) reduced the efficacy of marine reserves, which were situated in the spawning grounds.

| Dogfish
When fishing mortality was low, within-reserve population densities of dogfish were significantly higher than population densities without reserves (Table 5, Figure 2c, Fig. S9). When fishing mortality was high, population density continued to decline after reserve establishment, and within-reserve population density reached zero within 50 years in many instances, especially with little heritable variation in movement rate (Table 5, Figure 2c, Fig. S9). There were, however, significant increases in population density after reserve establishment when movement distances were low and there was sufficient heritable variation in movement distance (i.e., when μ AA = 200 NM and μ aa = 25 NM, or when μ AA = 100 NM and μ aa = 50 or 25 NM; Table 5). Marine reserves led to minimal increases in population density in nonreserve grid cells adjacent to marine reserves (Figure 5a).
Consistent with our predictions, there were significant increases in the frequency of the low-movement-rate allele following the establishment of marine reserves in all cases where there was heritable variation and population density did not go to zero (Table 5, Figures 2c and 5b, Figs S11 and S12). Movement evolution was more extensive with high fishing mortality than with low fishing mortality, and greater evolutionary T A B L E 5 Population size and frequency of low-movement allele for dogfish 50 years after establishment of marine reserves (average over 10 replicate simulations, unless otherwise indicated). Values that are significantly higher than control values (with no reserves), according to Dunnett's test of multiple comparisons to a control, are shown in bold and italic type (p < .05). Allele A was assumed to act additively throughout (columns in gray indicate cases with no heritability in migration rate) b Average of two replicate simulations (population density went to zero in eight of ten replicates); SD = 0.0569.

Number of reserves
c Average of five replicate simulations (population density went to zero in five of ten replicates); SD = 0.0367.
d Average of eight replicate simulations (population density went to zero in two of ten replicates); SD = 0.0296.
e Average of nine replicate simulations (population density went to zero in one of ten replicates); SD = 0.0281. changes were seen with higher heritable variation (i.e., when μ AA − μ aa was large; Table 5, Figure 2c, Figs S11 and S12). These evolutionary reductions in movement rate were generally associated with substantial increases in within-reserve population sizes over the course of the 50 years following reserve establishment (Figure 2c). The spatial pattern of allele frequency change indicates that the evolution of decreased movement was concentrated in and around reserves (Figure 5b).

| Great white sharks
Marine reserves did not cause significant increases in reserve population density except in some scenarios with the largest reserve size and low fishing mortality (Table 6, Figure 2d, Fig. S13). For scenarios with high fishing mortality (fished = 0.2), population density continued to decline after reserve establishment and was driven to zero by year 125 (25 years after reserve establishment) in all cases. The establishment of a reserve did not prevent this extinction (Figure 2d, Fig. S14).
When fishing pressure was low (fished = 0.05), marine reserves caused a slight increase in within-reserve population density relative to areas outside reserves when the marine reserve was large (Figure 6a). There was a trend of increased within-reserve population density with decreasing movement distance (Fig. S13).
There was very little evolution of movement rate in our simulations of great white shark. There were significant increases in the  Table 6; Fig. S15). These slight increases in allele frequency were localized within the marine reserve ( Figure 6b).

| DISCUSSION
We predicted that the evolution of decreased movement rate following the establishment of marine reserves would lead to higher within-reserve population size relative to situations with no heritable variation in movement rate. Our results suggest that the evolution of decreased movement rate can augment the efficacy of marine reserves and lead to higher within-reserve population density, especially for species, such as skipjack tuna, with relatively short generation times. Even when movement rates did not evolve substantially (e.g., given limited time relative to generation length and/ or with little heritable variation), marine reserves were an effective tool for the conservation of fish populations when movement rates were low or when reserves were very large (e.g., in the case of great white sharks).
We also predicted that higher fishing mortality (e.g., overfishing) would strengthen selection pressure and lead to faster evolution of movement rate. These effects were seen for every case where overfishing did not lead to extinction. Given the observed association between the evolution of decreased movement rate and increased within-reserve population size, we contend that marine reserves are an even more effective tool as an insurance policy (e.g., against inaccurate estimates of the maximum sustainable yield) than considered previously, without accounting for evolutionary changes.
Detailed and accurate descriptions of fish population sizes and movement patterns based on tag and recapture data (e.g., Sibert & Hampton, 2003) are crucial to understanding fish stocks, and fisheries scientist have made great strides in identifying the appropriate methods and types of data required for reliable fisheries management (e.g., Punt, Akselrud, & Cronin-Fine, 2017;Wetzel & Punt, 2015). It is generally the case, however, that the models of fish population dynamics upon which fisheries management is based do not include heritability of life history parameters (such as growth, maturation, or movement rates) and, hence, do not consider the effects of fisheries-induced evolution (Fraser, 2013;Palkovacs, 2011). Fish stocks do evolve in response to fishing pressure, and there is evidence that this evolution can have significant economic effects (Dunlop, Eikeset, & Stenseth, 2015). Among the studies that investigate the effects of fisheries-induced evolution, the focus has generally been on the consequences for growth, maturation,  (Dunlop et al., 2009). If movement rate is heritable, marine reserves may result in locally adapted within-reserve populations with reduced movement rate that are protected from selection for small maturation size (Miethe et al., 2009).
Importantly, our study suggests that the evolutionary effects of marine reserves are also applicable even to large pelagic fish such as tuna and sharks. In addition, we show that increasing the size or number of marine reserves enhances their overall efficacy for conservation purposes and that the efficacy is augmented by the evolution of movement rate regardless of reserve size or number. Our finding regarding the benefits of multiple reserves is in alignment with previous analyses that suggest the need for multiple closely spaced reserves to maintain connectivity and preserve populations of fishes, especially in light of increasing sea surface temperatures (Andrello, Jacobi, Manel, Thuiller, & Mouillot, 2015;Gerber, Mancha-Cisneros, O'Connor, & Selig, 2014). Overall, we contend that movement rate, like many traits, is more likely to be heritable than not and consideration of the evolutionary effects of marine reserves and fishing should be incorporated into fisheries management and spatial conservation planning.
In our simulations, the evolution of decreased movement rate within reserves was likely restricted by gene flow from outside of the reserves. Although movement rate is neutral outside of the reserves, and there is, therefore, no selection favouring the high movement rate allele outside the reserves, the strength of selection, averaged across the whole population, will still be diminished by gene flow. Hence, in order for adaptation to proceed rapidly in the reserves, selection must overwhelm the effect of gene flow from the surrounding populations; otherwise, gene swamping will occur (Gomulkiewicz & Holt, 1995;Gomulkiewicz, Holt, & Barfield, 1999;Holt & Gomulkiewicz, 1997).
Selection was least effective when movement rates were so high that swamping was substantial and few fish would remain within the reserve, especially if genetic differences in movement rate were not substantial (e.g., compare thickest to thinnest arrows in Figure 2). As a consequence, species with lower movement distances (e.g., dogfish) tended to exhibit greater evolutionary responses over the 50-year time window than highly migratory species (e.g., bluefin). We did not include any larval drift (between cells) in our models, and the amount of swamping may, therefore, be underestimated in our study. An analysis of the effect of the distance and direction of larval drift on the evolution of movement rates is a potential avenue for future research on marine reserves.
Whereas the focus of our research was on the conservation implications of marine reserves, we acknowledge the potentially important benefit of marine reserves to fisheries from spilloverincreases in catch adjacent to reserves (Roberts, Bohnsack, Gell, Hawkins, & Goodridge, 2001). This is an interesting topic in the context of evolving movement rate. We have shown that high heritability of movement rate results in locally adapted philopatric populations with high population densities within reserves. Whereas increased population density within reserves should result in more spillover (Halpern & Warner, 2002;Roberts et al., 2001), reduced movement should result in less per-capita spillover. The net effect (of increased density and decreased movement) on spillover is unclear, and a quantification of these effects on catch rates adjacent to marine reserves would constitute another interesting avenue for future research.
We did not model any particular fitness benefit associated with high movement distance and did not account for the advantages of movement to the survival and fitness of pelagic fish. One could hypothesize the existence of several such benefits. For example, high movement distances might have evolved due to a need to exploit spatially and/or temporally fluctuating resources. We assumed resources were fixed and homogenous across the simulation grid. In addition, we did not model any particular fitness cost associated with low movement distance, but such costs might exist. For example, low movement distances may cause negative density effects on populations (e.g., increased adult mortality or lower fecundity), as more stationary adults would tend to accumulate in particular locations (such as within reserves). We assumed that only fry survival was density dependent.
Future work on marine reserves could investigate the effects of additional density dependence (e.g., density-dependent adult survival or movement distance), as well as the effects of fluctuating resources.
Although such factors would undoubtedly have a quantitative effect on the results, the qualitative impact of the strong selection on migration due to high fishing pressures is likely to augment the effectiveness of marine reserves through the evolution of movement rate across a broader range of assumptions than considered here.
While our simulations are highly simplified representations, our study makes two important points that are generalizable. First, whereas the effectiveness of marine reserves may vary for different species of pelagic fishes depending on their life history characteristics, this effectiveness is most certainly underestimated if the potential for evolution of movement rate is neglected. We note that, in our simplified model of movement evolution based on a one-locus two-allele system, heritability was relatively low (maximum heritability was 0.07), and hence our estimates of the time scale of evolution are conservative relative to cases with higher heritability (e.g., when multiple genes contribute to movement rate). Second, if evolution of movement rate occurs due to fishing mortality outside of reserves, the rate of evolutionary change will be highest when mortality is highest (i.e., when selection is strongest), and marine reserves can therefore act as an important in-