Projecting the timescale of initial increase in ﬁshery yield after implementation of marine protected areas

,


Introduction
Marine protected areas (MPAs) are at the forefront of global marine conservation management and are being implemented widely to improve ecosystem health and fishery management (Sala et al., 2018).While controversial, MPAs are generally expected to achieve conservation goals, such as restoring demographic structure within MPAs, increasing spawning biomass, enhancing fisheries yield through "spillover", improving population stability and resilience, and protecting biodiversity (Pelletier and Mahevas, 2005;White et al., 2011), though the precise goals vary among MPAs.For example, in California, a network of 132 MPAs was implemented between 2003 and 2013 (Botsford et al., 2014) under mandated adaptive management.While the legislation enacting those MPAs did not include enhancement of fishery yield as a goal, there is public interest in how MPAs benefit fisheries.From a management perspective, it is important to know how soon such benefits could emerge in order to set monitoring expectations and plan decision-making.
Managers and stakeholders of MPAs are typically interested in knowing how soon expected effects of MPAs on fished populations will become evident inside and then outside MPAs (Halpern and Warner, 2002;Claudet et al., 2008).This interest is addressed when MPAs are managed through adaptive management, an approach to resource management that requires V C International Council for the Exploration of the Sea 2021.This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/ licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
quantifying expected benefits of management actions, then monitoring to see whether the expected results were attained (Walters, 1986;Parma et al., 1998).These results may suggest modifications to management.Thus, adaptive management of MPAs requires projection of the time-course of expected response and deciding what will constitute evidence that MPAs are successfully meeting their goals (Grafton and Kompas, 2005;White et al., 2011;Kaplan et al., 2019).
Setting adaptive management expectations for the fishery benefits of an MPA requires the ability to predict the dynamic response of fisheries yield to MPA implementation (note that hereafter and throughout the paper, we are referring specifically to MPAs that prohibit take of the species of interest, also known as no-take marine reserves).Empirical studies have shown a net movement of adults and larvae from no-take protected areas to harvested areas, the so-called spillover effect of larvae from the MPA to the fished area (Pelc et al., 2010;Di Lorenzo et al., 2020).Two types of spillover can influence yield: (i) dispersal of pelagic larvae spawned in the MPA to the protected area and (ii) movement of juveniles and adults from the MPA to fished areas.Larval spillover is difficult to quantify empirically because of the challenges in tracking pelagic larvae, but movement from MPAs to fished areas has been detected in a few locations using genetic parentage analysis (Harrison et al., 2012;Le Port et al., 2017).There has also been direct observation of adult spillover (e.g.Lowe et al., 2003) and subsequent catch of large adults outside protected areas (e.g.Roberts et al., 2001).Adult spillover is necessarily limited to occurring at the boundary between an MPA and a fished area (and can influence patterns of fishing there; Kellner et al., 2007), whereas larval spillover can originate anywhere within an MPA and contribute to the fished population over a broad spatial scale.Given this greater relevance of larval spillover to fishery dynamics at the regional scale, we focus exclusively on that type of spillover here.
To date, investigations of the dynamics of fished populations inside MPAs and the associated fishery yield outside MPAs have focused on longer-term, steady-state outcomes (McGilliard et al., 2011;White et al., 2011;Botsford et al., 2019), with a few notable exceptions analysing short-term responses of yield outside the MPAs (e.g.Babcock and MacCall, 2011).For example, Kerwath et al. (2013) noted an increase in fishery catch-per-unit-effort near a South African MPA and proposed an informal, verbal population model to explain the timescales associated with that pattern.Hopf et al. (2016a) formulated age and spatially structured models of Australian coral trout (Plectropomus leopardus) to understand the range of possible short-term responses of fishery yield to MPA implementation.Their model projections agreed with observed trends in coral trout fishery yield following the rezoning of the Great Barrier Reef Marine Park (Hopf et al., 2016b), but they did not explore the general question of how those patterns might vary for other species with different life histories.
An effective and meaningful way of assessing the timescales of the short-term, "transient" dynamics associated with MPA implementation is to focus on the time needed to recover or "fill in" the fished population age structure (White et al., 2013;Easter et al., 2020)."Filling-in" is the process by which the age structure that has been truncated by age-(or size-) selective fishing eventually returns to its unfished state as newly protected fish survive longer and re-populate the older age classes.Filling-in is the first measurable evidence that an MPA is reversing the effects of fishing, and it is a necessary precursor to most other projected benefits (i.e.those depending on increased abundance, biomass, or reproduction inside the MPA).Analysis of age-structured population models reveals that the ultimate increase in fish abundance inside the MPA, relative to pre-MPA abundance, is 1 þ (fishing mortality rate)/(natural mortality rate) (White et al., 2013).This corresponds to the intuitive idea that more heavily fished populations will exhibit a greater relative increase.The time required to reach that value is less with a lower natural mortality rate.These results have been applied to MPA planning and assessment processes in California, USA (CDFW and OPC, 2018) by setting expectations for how soon (if at all) detectable increases in fish abundance should be observed in different MPAs (Nickols et al., 2019) and by providing a basis for choosing which species to monitor for most easily detectable changes (Kaplan et al., 2019).
Here we address the response of fisheries yield as a second goal in the adaptive management of MPAs.Specifically, we characterize the timescale over which the filling-in of the population age structure within MPAs would lead to the production of larvae that could spillover to fished areas and increase yield there.We focus on the timescale because the eventual magnitude of any increase in yield from an MPA will depend on a variety of uncertain and location-specific factors, including patterns of larval connectivity and how the fishery is managed (as shown by Hopf et al., 2016a).However, as we show, the question of the timescale over which increases in yield would arise has a more general solution.We estimate the timing of filling-in and the consequent increase in fishery yield for a number of exemplary species, demonstrating how and when the increase in yield derived from MPAs will occur.For this study, we use a renewal equation formulation of an age-structured, two-patch model to describe the response of yield in the fished patch after MPA implementation.The 16 species we analysed are all managed as part of California's Nearshore Fishery Management Plan (CDFG, 2002).Our results show why the timescales of increases in fishery yield are longer than those of biomass increase inside the MPA and how the additional delay depends on fish life history and fisheries management.

Methods
For most demersal nearshore fishes in MPAs, the spatial scale of adult movement and the size of the MPA itself are small relative to the spatial extent of larval dispersal.Consequently, earlier theoretical investigations of transient dynamics in MPAs (White et al., 2013;Kaplan et al., 2019;Nickols et al., 2019) have considered the case of a single habitat patch (the MPA) embedded in a larger coastal meta-population.Most (or all) of the larvae recruiting to the MPA patch are produced elsewhere in the meta-population, so the patch has "open" demographics, and additions to abundance and biomass inside the MPA are due to filling-in of the age structure, rather than changes in reproduction.The goal of this approach was to examine the first-order, short-term dynamics inside individual MPAs, absent the uncertainty introduced by considering larval connectivity.
For an analysis of fishery yield, we extend this model framework to consider two adjacent patches embedded within a metapopulation: one patch that becomes an MPA and one that remains open to fishing (Figure 1).Because these patches are small relative to the size of the larger meta-population, we retain the open-recruitment dynamics, with a constant influx of recruits to both patches.We now add one additional process: as reproductive biomass in the MPA patch fills in, larval production in that patch also increases.Some proportion of those additional larvae will disperse to the fished patch, where they will grow and eventually enter the fishery (Figure 1).Our analysis here focuses on the timing of the resulting increase in yield in the fished patch.The magnitude of the increase will depend on larval connectivity between the MPA and the fished patch, which is highly uncertain, as well as the fishing mortality rate.Therefore, we focus solely on the timing of the increase which, as we will see, depends only on the life history.We also do not include the additional recruitment that would return to the MPA and contribute to further increases in biomass there, because this also depends on uncertain larval connectivity rates (in this case, the probability of self-retention).That effect would necessarily occur after the initial increase in yield due to filling-in and would require an additional lag of the age of first reproduction.Thus, it would be a second-order effect having only a small consequence on initial timing of the yield response.We focus solely on the first-order timing of increases in fishery yield attributable to filling-in inside the MPA.
A typical way to examine the dynamic consequences of the biomass increase due to filling-in of the age structure following MPA implementation is to use age-structured models (e.g.Leslie matrix models) to project forward the change in the abundance of individuals in each age class (White et al., 2013;Hopf et al., 2016b;Kaplan et al., 2019).However, here the increase in fishery yield (in the fished patch in Figure 1) will be driven by the new recruits produced in the MPA that disperse to fished areas, and their contribution to the fishery will be determined by the timelags associated with the growth and survival of the cohorts from these recruits as they enter the fishery.Fishery yield from a cohort depends on how large (and how old) the fish in the cohort are allowed to grow before being caught.It is, therefore, more illuminating not to use the Leslie matrix approach, but to model this process using an alternative age-structured model to describe these recruitment-driven, age-structured population dynamics.Rather than projecting forward all of the age classes each year (as in a Leslie matrix), we project the new recruits from the MPA patch entering the fished population in each year and their eventual fate in the fishery.This equivalent formulation is known as a renewal equation because it keeps track of new entrants to the population (Lotka, 1907;Botsford et al., 2019).

A renewal equation for yield
A renewal equation model is an age-structured model that avoids having to directly track the abundance in each age class.Instead, it keeps track of the abundance of each age class in a cohort in terms of the recruitment to that cohort (which can vary over time) and the proportion of recruits in each cohort that survive to a given age (which is a function of age and is assumed to not vary with time).For example, the abundance of fish of age a at time t (N a,t ) would be the number of recruits a years prior to t (R t-a ) multiplied by the probability of survival from recruitment to age a (S a ).
In this modelling framework, age-dependent population processes, such as reproduction and fishery yield (both of which depend on survival to a given age and the biomass of individuals at that age), can be expressed concisely in terms of past values of recruitment.The mathematical representations of these agedependent processes are named influence functions because they represent the influence a cohort of offspring has on some aspect of population dynamics in each successive age of the cohort (Botsford et al., 2019).
For example, the description of the effect of each recruit on yield as the recruits' cohort ages turns out to be the expression for the yield per recruit (YPR) commonly used in fisheries management (Beverton and Holt, 1957;Botsford et al., 2019).For a given age a, YPR a is the product of the probability of survival to age a (S a ), biomass at age a (W a ), and the fishing mortality rate F. The survival term is comprised of the probability of survival to the age of entry (a c ), when the population just experiences natural mortality M, and the probability of surviving both natural and fishing mortality afterwards: with YPR(a) ¼ 0 for a < a c .We refer to this as the yield influence function (All variables used in this paper are described in  Proportional size of the MPA patch Table 1).Here we assume that W a is proportional to length-atage cubed, and that length-at-age follows the von Bertalanffy relationship [see Botsford et al. (2019) for a description of the von Bertalanffy equation for the indeterminate growth of fishes].This influence-function approach to describing dynamics is useful in situations when recruitment changes substantially from year to year and variability in recruitment is greater than any change in post-recruitment processes.Thus, the inter-annual variability in recruitment can be translated into its eventual effects on population dynamics using the influence functions.
Applying the renewal equation approach to describe the yield in the fished patch that will result from the filling-in of the age structure in the nearby MPA patch requires: (i) a description of the number of additional new recruits that will initiate a cohort in the fished patch and (ii) the influence function for the yield that will be realized from that additional recruitment.The former (i) is based on the expression for how reproductive biomass (B t ) increases over time inside the MPA relative to its value at implementation (B 0 ) as the age structure fills in, because larval production is generally proportional to reproductive or spawning biomass in fishes.The second requirement, the influence function, is the expression for YPR a in Equation (1).Kaplan et al. (2019) provided an expression for the increase in biomass (of the fished age classes) and showed that it depends on the natural mortality rate M, the pre-MPA fishing rate F, the age of entry to the fishery a c , and the growth rate k (all in units of time À1 ).That expression is: where the indexing variable K i takes on values of 1, -3, 3, and -1 for i ¼ 0, 1, 2, 3 [this formulation is the expansion of the cubic polynomial that arises from assuming that biomass is proportional to length cubed, as in Beverton and Holt (1957); see additional explanation in Botsford et al. (2019) and Kaplan et al. (2019)].The growth rate k is a parameter in the von Bertalanffy growth function.
The premise of our analysis is that the increase in biomass inside the MPA translates into an increase in larval production, some proportion of which, denoted c U (to indicate a connectivity probability that is uncertain), is exported to recruit in the fished area.The term c U is intended to correspond to the probability of connectivity between patches that appears in most spatial metapopulation models.Our analysis (Equation 2) assumes that reproductive output is proportional to biomass and that biomass is proportional to length cubed.Those assumptions allow us to describe the system's dynamics using the relatively simple formulations of Equations ( 1) and (2), revealing how the results essentially depend on a few key life history parameters.In actual applications, when those assumptions do not hold and parameters are available to calculate larval production from the filling-in of biomass more accurately (e.g.Barneche et al., 2018), they can be accounted for in the model formulation.We also conducted full age-structured population simulations for each species (methods described in the Supplementary data) using speciesspecific life history and fishing-related parameters to confirm that the timescale results of the recruitment index obtained with Equation (2) were indeed consistent with the more complicated modelling approach (see Supplementary Table S1 and Figures S1  and S2 in the Supplementary data).In the Supplementary data, we also show how slight the differences are that arise when reproduction is not proportional to length cubed.
Given the expression in Equation ( 2) for the increase in reproductive biomass in the MPA, the resulting value of recruitment to the fished patch, in addition to what it would have received without the effect of filling-in in the MPA, is: (3) Because of the uncertainty in connectivity, we will refer to R t as a recruitment index to emphasize that, in practice, one would be unlikely to know its actual magnitude, but that it is useful to know whether and how quickly it is increasing over time.Here, as elsewhere in ecology, we use the word "index" to mean a quantity proportional to the value, with the constant of proportionality unknown (e.g.catch-per-unit-effort as an index of abundance).
We combine these two functions, recruitment (Equation 3) and yield-per-recruit (Equation 1), to predict the timing of yield increases.We note that, at any time, the total yield Y t will be the sum over all ages that are harvested of the yield from each harvested age a.The yield at time t from cohorts that are currently of age a is found by taking the recruitment a years ago (R tÀa ) and multiplying that by the yield from that cohort at time t (YPR a ).Summing this product over all ages that are harvested gives the total yield at time t: where A is the maximum age for the species.Here, Y t represents only the additional yield due to larval spillover during the fillingin process; recruits arriving from elsewhere in the metapopulation support the baseline level of yield that was present prior to MPA implementation.The magnitude of Y t is also uncertain because the value of c U is unknown because of our uncertainty about larval connectivity patterns.We will refer to Y t as the added yield index to indicate this.Note that this formulation does not account for any density-dependence in the system (e.g.density-dependent recruitment), though such effects would not be expected to alter the timing of the initial increase in added yield and could be considered some of the uncertainty represented by c U.
The relationship between the recruitment index and the added yield index is illustrated graphically in Figure 2. Figure 2a is an example of the increase in relative biomass (and the recruitment index) described in Equation ( 2) for an MPA that goes into place at t ¼ 0. Figure 2b is an example of the yield influence function YPRa the contribution per recruit to yield at each age (Equation 1).It is zero for the first several years before fish are large enough to enter the fishery, then rises to a peak and declines as individuals are caught by fishing and die from natural mortality.Figure 2c illustrates how the total fishery yield at any given time is the sum of the yield influence functions for every cohort present in the population at their current age a.In this example, cohorts that recruited in years 1, 4, 7, and 10 would be 9, 6, 3, and 0 years old at t ¼ 10; the corresponding yield influence functions for those cohorts are shown, with an open circle indicating the value of yield at t ¼ 10.The additional yield (provided to the fished patch) would be a sum of those values at t ¼ 10.The example in Figure 2c assumes that recruitment is constant over time, so each cohort has the same abundance.If recruitment is increasing over time (as it would be if there were spillover from an MPA), the yield influence function for each successive cohort will be multiplied by the recruitment index for that year.This is Table 2. Life history and fisheries management parameters and biomass change through time after no-take reserve implementation and resultant fishery yield lag timings from population modelling of 16 species of California nearshore species.illustrated in Figure 2d: the added yield at any time (the purple curve) is the sum (vertically) over all of the R t Â YPR functions (the blue curves) at that time.
A key observation is that when the values depicted in Figure 2d are summed over age at a particular time (e.g. at t ¼ 25), the subscript for recruitment goes back in time (t -a, for increasing values of a) as the YPRa subscript goes forward in age (increasing values of a).That is, a given cohort is older at a particular time t (with respect to the influence function) and also recruited at an earlier time.This means that one can determine the total added yield at a given time t by computing a weighted moving average of R t over all the prior years.The weightings will be the yield Figure 3. Calculation of timing of increase in added yield for four species, each denoted by colour: kelp rockfish (purple), blue rockfish (blue), gopher rockfish (red), and kelp greenling (yellow).(a) Expected increase in biomass in the MPA in the years following MPA implementation (relative to the biomass at the time of implementation); the additional recruitment to the fished patch is assumed to be proportional to this increase.Round symbols on each curve mark the time at which 95% of the asymptotic maximum is reached.(b) Yield influence functions for each species; the round symbols mark the mean age of each.(c) Index of added yield in the fished patch (the convolution of the recruitment index and the yield influence function).Round symbols on each curve mark the yield lag time at which 95% of the asymptotic maximum is reached.(d) Same as panel (c), but with a shorter time-horizon to reveal detail.

Figure 4.
Yield influence functions for the 16 California nearshore fish species.Vertical line in each indicates the mean age of the function.Differences in the shape of the functions correspond to differences in natural mortality and fishing mortality rates; the location of the left edge of the function corresponds to the estimated age at entry to the fishery.
influence function for the age at time t of each prior cohort.This explains the form of Equation ( 4).Because this is a weighted moving average, with the weightings biased to the right [because YPR(a) is zero until age a c and then has a long tail], the increase in yield will always lag the increase in biomass in the MPA.
For the mathematically inclined, we add that Equation (4) could also be understood by noting that the sums in that equation have the form of a convolution (Weisstein, 2020), which is a way of expressing the idea that the yield at a given time t is the sum of YPR for individuals of all possible harvested ages and an individual that is age a at time t started life a years before at time t -a.For additional details on this interpretation of the analysis, see the Supplementary data.

Application of the two-patch model
We use Equation ( 4) (Figure 2d) to calculate the expected timing of increase in fishery yield for the 16 species of nearshore fishes from California (as in Figure 2d).We based our calculations on the life history parameters available in the literature for each species and the estimate of fishing mortality rate (F) from fishery stock assessments (Table 2).Given both uncertainty in larval connectivity (c U ) and great variability among MPAs in the level of harvest prior to protection [and the expected magnitude of biomass increase; e.g.Nickols et al. (2019)], we focus primarily on the timing of the increase in yield rather than its magnitude.We compared the timing of the different species' expected increase in yield in terms the time required to reach 95% of its ultimate, asymptotic level of relative yield.We refer to this as the yield lag time.

Results
We calculated model results for 16 species of demersal nearshore groundfish, but for simplicity in the main text, we focus on the results for four representative species that show a range of possible outcomes; results for the others are shown in the Supplementary data.These four species are kelp rockfish (Sebastes atrovirens), gopher rockfish (Sebastes carnatus), blue rockfish (Sebastes mystinus), and kelp greenling (Hexagrammos decagrammos).In Figure 3a, we show the predicted response of biomass for each species to MPA protection, as calculated by Kaplan et al. (2019) based on stock assessment estimates of fishing mortality F. These range from an approximate threefold increase for blue rockfish to a less than two-fold increase for kelp greenling (Figure 3a).The four species included represent a range of the life history factors affecting the timing of yield response, including the age of entry to the fishery a c and the natural and fishing mortality rates M and F, respectively.These differences were reflected in the shapes of the yield influence functions (Figures 3b  and 4).As noted above in Equation (2), the contribution to yield at the first several ages (i.e. a < a c ) is zero, then the yield contribution jumps up to the biomass of an individual of age a c , and from then on the contribution was controlled by the interaction between increasing individual biomass and declining abundance due to fishing mortality and natural mortality.The range of shapes of this function for the 16 California fish species is shown in Figure 4.
The responses of the added yield index (Figure 3c) represent the combined effects of the different yield influence functions (Figure 3b) with the recruitment index (Figure 3b) for each species, as given in Equation (3).As such, the responses include the effect of the constant c U describing uncertainty in larval connectivity.In order to emphasize our interest in the timing of the yield response (but not its magnitude), we chose values of c U for each species so that the recruitment index had an asymptotic maximum value of 1.The shapes of the resulting transient added yield indices (Figure 3c and d) depended on two characteristics of the yield influence function.Their starting time depended on a c , the age of first capture, and the rate of increase in yield depended on the rate of increase in the biomass function (Figure 3b) and how much the width of the yield influence function (Figure 3a) extends that increase: the wider the yield influence function, the more it will slow down the response.This is because the yield obtained from a cohort is spread over a broader range of ages (a more mathematical explanation is that this effect arises because convolution is essentially equivalent to a weighted moving average, with the weightings being the YPR a functions flipped about age zero; thus a wider YPR a function smoothes out the moving average more).For example, the added yield index of kelp rockfish begins to increase first because it has the earliest value of a c , but its increase is very slow because its yield influence function YPR a (Figure 3b) is quite broad.Kelp greenling and gopher rockfish both had steeper responses in Figure 3c and eventually crossed over those of blue and kelp rockfish; despite greater values of a c , the former species having the most sharply peaked yield influence function in Figure 3b, so increases in yield were realized more quickly after a c .In Figure 3d, we repeated the plot of the yield indices over a shorter period to clarify the differences in when managers and stakeholders could expect the effects of MPAs on yield to begin to appear.Keep in mind that the results in Figure 3c and d reflect only the additional yield added to the system due to reproduction inside the MPA.The overall yield in the entire system would initially decline due to the MPA closure before eventually increasing, but that dynamic (of the entire system, including the yield in the area that is closed) is not shown here.
It is also useful for managers to know ahead of time when the effects of MPAs on yield will approach their maximum value and how much longer they will have to wait after seeing a response in biomass inside the MPA to see a response in yield outside.We plotted the yield lag times (the time required to reach 95% of the maximum value of the yield index) in comparison with the time to reach 95% of the maximum values of the MPA biomass ratio B/B 0 [equivalent to the values calculated for biomass by Kaplan et al. (2019)].The yield lag time ranged from 19 years (kelp greenling) up to as long as 62 years for the long-lived china rockfish (Sebastes nebulosus; Figure 5).
Knowing that the added yield index is the result of a convolution of the two functions for recruitment and yield-per-recruit suggests that the lags between the increases in recruitment and yield indices will be the dominant timescale of the YPR a function.To test this concept for the 16 species, we computed the correlation between the extra lag for yield effects (i.e.difference between values in each bar in Figure 5) and a measure of the width of the yield influence function, i.e. its mean age.The mean age was calculated as the sum over age of each age multiplied by the proportional yield influence function at that age; this is similar to the way a measure of generation time can be calculated from a reproductive influence function (Botsford et al., 2019).This would be interpreted as the mean age at which the biomass of a cohort is harvested.The correlation between the mean age of the YPR function and the yield time-lag was very strong (Pearson correlation, r ¼ 0.96; Figure 6).

Discussion
We developed a simple, two-patch model of the increase in yield from MPA implementation and used it to project the timescales of the increase in fishery yield after the implementation of an MPA.We demonstrated how the timing of increase in yield depended on the combined effects of: (i) the timing of increase in larval production inside the MPA following implementation and (ii) the way that fishery management affects the age structure of the fishery (i.e. via the age of entry to the fishery a c and the fishing rate F).The latter effect involves a concept, YPR, which is well-known in fisheries management (Beverton and Holt, 1957;Botsford et al., 2019).We found that the expected increase in fishery yield beyond MPA boundaries would take an additional 7-18 years beyond our previous projected increase in biomass following implementation (Kaplan et al., 2019) for 16 representative US Pacific coast nearshore groundfish fisheries.Importantly, these results can be applied to any system, because the length of the projected delays is predictable from the mean age of the YPR function (which is possible to calculate for most species, given mortality, and size-at-age data).However, though the timing of fishery benefits is predictable, their magnitude is not because of the uncertainty in larval export from MPAs to fished areas.From a conceptual point of view, our results improve our understanding of how spatial (i.e.MPAs) and conventional approaches to fishery management are fundamentally linked by the agestructured dynamics of fished populations (Botsford et al., 2019).There are a number of descriptions of how effects of MPAs and conventional fishery management on long-term yield are, in some sense, doing the same thing (limiting fishing effort) in different ways (Mangel, 1998;Hastings and Botsford, 1999;White et al., 2010;Botsford et al., 2019).Our results here address a different aspect, the timing of the short-term, transient response of yield to MPA implementation.Our analysis here adds an additional time-lag to the previous projected timescale of increase in MPA biomass following implementation, which was 9-43 years for the California species analysed here (Kaplan et al., 2019).The additional delay for yield ranges from 7 to 18 years.The duration of these additional projected delays was on the timescale of the mean ages of the yield influence functions (known in fisheries as YPR functions), which is possible to calculate for most species, given mortality, size-at-age data.However, we caution that the magnitude of eventual fishery benefits is more uncertain than those for biomass, given the uncertainty in larval connectivity from MPAs to fished areas.These results contain some element of the earlier results showing equivalence between MPA and conventional management (Mangel, 1998;Hastings and Botsford, 1999) in the sense that higher fishing rates prior to MPA implementation lead to greater increases in yield after MPA implementation, because there is greater scope for increase in biomass inside the MPA and thus spillover from the MPA will be relatively higher.
From a practical point of view, managers have found previous predictions of the amount of increase in MPA biomass and its timing useful in planning which species and locations to monitor (CDFW and OPC, 2018;Kaplan et al., 2019).However, our predictions of the timing of fishery yield provide a different type of information.The increase in fish abundance to the level of completely filling-in the age structure that had been truncated by fishing could be easily approximated by (1þ F/M) as long as the local value of F was known (White et al., 2013(White et al., , 2016;;Kaplan et al., 2019;Nickols et al., 2019).This means that locations and species with high fishing effort prior to MPA protection would be attractive as places and species to plan to sample (CDFW and OPC, 2018).Our extension of this theory to project increases in the timing of fishery yield is not intended to provide specific predictions of yield, but rather an index of the timing of yield increase.The actual level of fishery yield realized will depend both on highly uncertain larval connectivity pathways from the MPA to fished areas and on the behaviour of the fishing fleet itself (e.g.Kellner et al., 2007).Nonetheless, a key question is how soon increases in yield could be expected, and we have shown that the timing can be estimated if one knows the form of the yieldinfluence function in the location in question.
Our results are related to earlier studies that have empirically quantified the transient (i.e.short term) response of fishery yield to MPAs using age-or size-structured models.Specifically, Kerwath et al. (2013) examined the fishery for a South African temperate reef-associated fish (roman, Chrysoblephus laticeps) and noted that the observed delay in increase in yield after MPA implementation should at least be greater than the age of entry to the fishery.The age at first catch for that species is 5 years (Kerwath et al., 2008), so Kerwath et al. (2013) attributed all increases in yield before 5 years to a mechanism involving diffusion of adults from the MPA.In their 2013 study, Kerwath et al. found that in one no-take reserve (Goukamma MPA), yield increased to twice the initial amount after 10 years following implementation of the reserve.This is approximately what our model would have predicted for this species: a 5-year lag given the age at first capture, followed by a further $5-year lag driven by the distribution of the yield influence curve for this species.
Theoretical assessments of how fishery yield may respond to MPAs over short timescales are limited.In one such study, Hopf et al. (2016a) showed that in a closed population with well-mixed dispersing larvae, yield always decreased after MPA implementation, but then recovered to previous levels after a 5-to 10-year lag and continually increased afterwards.The long-term outcome of that model depends on their assumptions of geometric growth without density-dependence and fishery management but our model framework would produce similar results for the timescale of increases in fisheries yield, given the life history of the tropical grouper species modelled by Hopf et al. (2016a).Additionally, had we specified the relative sizes of the fished and MPA patches in our analysis so that we could plot the total yield of the entire system, the resulting plot would have shown the well-known initial drop in total yield following implementation described by Hopf et al. (2016a).
The use of the renewal equation to describe the dynamic behaviour of recruitment, with influence functions representing the consequences of recruitment at older ages, may be useful in analyses by others working on populations distributed over space and linked by larval transport.These models have been used in other non-spatial contexts in which variability in recruitment is a central issue (Botsford et al., 2019).Analysis using this approach is often easier and clearer when temporal variability is confined to a single variable like recruitment, and population dynamics then depends on the magnitudes and shapes of the influence functions.For example, influence functions have represented lifetime reproduction over age and cannibalism over age in crab populations (Botsford and Wickham, 1978) and have been used to formulate another influence function, a "net area function" in barnacle populations where recruitment depended on available free space and the influence function was size (and thus space occupied) as a function of age (Roughgarden et al., 1985).
Our study provides an expectation of the amount of time required to observe increased fishery yields due to larval spillover from an MPA.There is growing empirical documentation of larvae produced in MPAs that have dispersed to fished areas (Pelc et al., 2010;Harrison et al., 2012;Johnson et al., 2018;Baetscher et al., 2019), but until now there was not a rigorous expectation for how soon to expect evidence of that spillover to begin.Our mathematical approach, based on the yield-per-recruit concept, provides a framework for comparing species' expected responses, based on their life history and fishery management.These results can guide the adaptive management of MPAs (and their comanagement with conventional fisheries; White et al., 2011;Carr et al., 2019) and provide a stepping stone to future analyses of higher-order dynamics of spillover and yield.Such analyses will require a better understanding of the larval transport pathways from MPAs, as well as variability in those pathways [reviewed by White et al. (2019)], in order to predict both the timing and the magnitude of changes in yield.

Figure 1 .
Figure1.Schematic of the two-patch population model, including constant recruitment arriving in each patch (proportional to patch area).As biomass increases in the left patch after it becomes an MPA, additional larvae are spawned there, some of which disperse (purple arrow) to the fished patch, where they contribute to additional fishery yield.

Figure 2 .
Figure 2. Diagram of the model components in Equations (1-4).(a) Increase over time in biomass inside the MPA patch translates into an increase in the index of recruitment exported to the fished patch.(b) A yield influence function showing the contribution to YPR of each age class.(c) Schematic showing how the fishery yield at any point in time (here, t ¼ 10 years) is the sum of the added yield index and influence function (at the appropriate age) for each cohort present at that time.Here, four cohorts are shown, recruiting at years t ¼ 1, 4, 7, and 10, so that they are ages 9, 6, 3, and 0, respectively, at t ¼ 10.The yield at t ¼ 10 is the sum of the values in the open circles for each curve.This corresponds to the summation in Equation (3) of recruitment t minus age a years in the past times the yield influence function for age a (assuming recruitment is the same magnitude each year).(d) As the recruitment index increases over time, the product of recruitment and the yield function increases for each successive cohort (blue curves); summing at each point in time produces the expected trajectory of the added yield index (purple curve).

Figure 5 .
Figure 5.Time required to reach biomass time (i.e.time for "filling in" biomass in the MPA; lower bars) and yield lag time in the fished patch (upper, longer bars) for 16 California nearshore fish species.The species are arranged in order of decreasing natural mortality rate (M) from top to bottom.

Figure 6 .
Figure6.The relationship between the difference in the yield lag time and biomass lag time and the mean age (yrs) of the yield influence function for each of the 16 nearshore fish species.

Table 1 .
Description of variables used in this study.