CHANGE IN BIOLOGICAL REFERENCE POINTS UNDER DIFFERENT BIOLOGICAL , FISHERY , AND ENVIRONMENTAL FACTORS

Background. In this study, the relation between the biological reference points (BRPs) associated with Maximum Sustainable Yield (MSY) and the relevant biological, fisheries, and environmental factors were investigated. This knowledge is crucial to build the capacity for timely adaptation of management to the changes in an ecosystem. The research highlights the considerations that need to be taken into account when estimating and using BRPs in practice, and can thus lead to avoidance of overfishing. The BRPs were represented by MSY and the related spawning stock biomass (BMSY), fishing mortality (FMSY), and FMSY proxies. Materials and methods. To obtain the BRPs, the method of Horbowy and Luzeńczyk (2012)—that combines yield-per-recruit and spawning stock-per-recruit analyses with stock-recruitment relation models (Beverton and Holt (B&H) 1957, Ricker 1975)—was used. Data from the three biggest Baltic Sea stocks were used to test the sensitivity of BRPs to the model input data, and the influence of regime shifts or dynamics of the biological and fisheries variables on the BRPs. Results. The analyses show that an increase in maturity and weight at age generally led to an increase in BRPs. The opposite effect is observed in the case of natural mortality (stronger when the Ricker stock–recruitment (S–R) relation was used) and selectivity (proportion of the F of partially recruited age groups to the mean F of fully recruited age groups). An increase in the steepness of S–R models increases FMSY and its proxies as well as MSY, and had a decreasing influence on BMSY. Furthermore, the BRPs are very sensitive to the different data range divisions of the input data. Conclusion. When estimating BRPs, the time period of the input data should be selected with caution, and the appropriate time period should be based on strong biological, ecological, and environmental knowledge as, according to this study, all of these factors influence estimates of BRPs.


INTRODUCTION
Maximum sustainable yield (MSY) is a fisheries management concept that originated in the 1930s.Over the years, this concept was alternately popular and criticized (e.g., Larkin 1977, Sissenwine 1978, Mace 2001).The definition of MSY has evolved; for example, the perception of fishing mortality associated with MSY (F MSY ) has developed from a target that could routinely be exceeded into an avoidable limit (Anonymous 1995, Mace 2001).In 2002, during the World Summit of Sustainable Development, the participants declared that all stocks should be maintained or restored to levels that could produce MSY by not later than 2015 (Anonymous 2002).As a result, MSY took on increasing importance as a worldwide management goal implemented by several jurisdictions, including the International Council for the Exploration of the Sea (ICES), which is the main advisory body of the EU and the member states.
ICES recommendations on future catches are now based on the MSY concept (Anonymous 2016a).If an estimate of F MSY is not available, a series of alternatives to F MSY reference points are used as conservative proxies (but not so that they result in substantial foregone yield).These include F max , which maximizes equilibrium yield per recruit; F 0.1 , the point at which the slope of the yieldper-recruit (YPR) curve is one-tenth the slope of the curve at its origin; F 20%-40%SPR , which produces 20%-40% of the virgin spawning stock-per-recruit (SPR) (i.e., SPR at F = 0); and F MSY , from production models (Anonymous 2013a).The newest ICES advice is to select the reference points for exploitation based on the lower confidence intervals of F MSY -rather than at the lower range of an interval, such as F 0.1 (Anonymous 2016a).
To estimate the MSY reference points, both the biological and fisheries characteristics of the stock are needed.The biological parameters are dependent on environmental factors, which in an indirect way influence the MSY reference points and can change them substantially.According to Möllmann et al. (2008), after analysis of 52 biotic and abiotic variables (including spawner biomass, recruitment, weight, fishing mortality of cod, herring and sprat, phytoplankton, zooplankton, salinity, and deepwater oxygen), the Baltic Sea had a major reorganization at the end of the twentieth century.They identified two states between 1974 and 2005, separated by a transition period during 1988-1993, called a regime shift.One of the definitions used to describe the regime shift is a change in the biological ecosystem in response to physical drivers (Collie et al. 2004, deYoung et al. 2004).In addition to a climatic regime shift, overexploitation is another very important factor, which can lead to surprisingly large changes in marine ecosystems (Scheffer and Carpenter 2003).In the late 1980s in the Baltic Sea, the regime shift was observed in the fish community when the system that was dominated by cod (Gadus morhua) changed into one with sprat (Sprattus sprattus) dominance (Anonymous 2013b).The factors responsible for that shift were initiated by climate-induced changes in the abiotic environment, which were further modified by anthropogenic impacts via fisheries-induced feedback loops in the food web (Möllmann et al. 2008).Reorganization in the system has important management implications, including changes in biological reference point values.Management that increases ecosystem resilience is not always possible, but building the capacity for quick adaptation to new situations is also beneficial (Crépin et al. 2012).
In this paper, the method of Horbowy and Luzeńczyk (2012) was used to obtain the biological reference points (BRPs) MSY, B MSY , F MSY , and its proxies for eastern Baltic cod (Subdivisions (SD) 24-32), herring in the Central Baltic (SD 25-29 and 32, excluding Gulf of Riga), and sprat in the whole Baltic Sea .First, the sensitivity of the biological reference points to the following input data to the model were tested: maturity (Mat), weight at age (W), natural mortality (M), exploitation pattern (Sel), and stock-recruitment (S-R) relation.Based on the sensitivity analysis results, the BRPs estimated for different data ranges (division based on the regime shift and dynamics of the biological and fishery variables) were investigated to find the linkage between them and the environmental factors, which could help to predict the BRPs by tracking environmental changes.

MATERIALS AND METHODS
The method of Horbowy and Luzeńczyk (2012), which combines YPR and SPR analysis with S-R models (Beverton andHolt (B&H) 1957, Ricker 1975), was used to derive the biological reference points.The equations for equilibrium spawning biomass (B eq ) and yield (Y eq ), using the B&H S-R relation, are where B eq (F) and Y eq (F) are the equilibrium spawning biomass and yield, respectively, at a given for the Ricker S-R relation.It is easy to use these equations to derive the fishing mortality that maximizes the equilibrium yield.The equilibrium yield can be maximized with respect to F to obtain biological reference points (Horbowy and Luzeńczyk 2012).Based on this method, MSY, B MSY , F MSY , and its proxies (F 40%B , F 95MSYlower , F 95MSYhigher ) were calculated.F 40%B is analogous to those from SPR analysis, but refers to the spawning stock equilibrium biomass; it is defined as the value of F for which the biomass is equal to 40% of the unfished spawning biomass.The variables F 95MSYlower and F 95MSYhigher are defined as the values of F that provide a catch of 95% of MSY; there are two fishing mortalities that comply with this condition, namely one lower than F MSY and one exceeding F MSY .
First, the basic runs to estimate BRPs were performed.They include the S-R relation and mean Mat, W, M, and Sel at age from the full available data series for each species (i.e., data from 1974-2015 for Baltic sprat and herring (Anonymous 2016b), and data from 1966-2012 for cod (2012 was the last year with an estimated SSB and R; Anonymous 2013b, 2014).Next, the sensitivity of the BRPs to the model input data was assessed using 10 additional options, including reducing or increasing the steepness of the S-R relation and the means at age of Mat, W, M, and Sel by 20 percentage points (pp).For cod, two different weights at age were used: the weight at age in the catch (Wc) and the weight at age in the stock (Ws); in the case of herring and sprat, those two biological variables are assumed as equal (Anonymous 2013b(Anonymous , 2016b)).
Additionally, runs to estimate BRPs with the highest and the lowest observed values for Mat, M, W, and Sel as well as +/-50% of h for each of the stocks were performed.For some of the cases (Mat for sprat and herring and M for cod), constant values were assumed for whole years of the stock assessment data, and no further investigation could be performed.
The next step was to use the knowledge about the sensitivity of BRPs to input data to interpret changes in the estimated BRPs for periods of differing environmental conditions.Two data series (before and after the regime shift) were used, with the division according to Möllmann et al. (2008), that is, before 1988 and after 1993, with a transition period of 1988-1993.In the case of Baltic sprat and herring, the first series used the input data from 1974-1987, and the second from 1994-2015 (Anonymous 2016b).For cod, the data series used for the analysis was divided into 1966-1987 (before the regime shift period) and 1994-2012 (after the regime shift; Anonymous 2013b, 2014).Another division of the period after the regime shift was assessed with the data from 1994-2005 for each species, as the changes in the last ten years could be significant and could influence the results.
Additionally, the general linear model (GLM) was used to detect changes in the time series input data (Mat, W, M, and Sel).Each input data series was separated into year and age effects according to the equation: dependent_variable = year + age + error, where error was log-normally distributed.The first year effect for each of the input data series was assumed as zero; the effects of the other years are the differences between each of the years and the effect in the first year.The significance of the differences and covariance were tested.Based on the GLM results, another division of the time series for each of the three species, taking into account the dynamics of the biological and fishery variables, was implemented to check the influence on the BRPs.

RESULTS
Growth and natural mortality of cod, herring, and sprat declined markedly in the periods considered.The growth declined by 40-60 pp over the full time period, while natural mortality was smaller by half in recent years than in the 1980s and 1990s.Sel for the partially recruited age groups fluctuated over the years, especially in the case of sprat.
The curves obtained by fitting S-R relations to the full available data series for each species (i.e., data from 1974-2015 for Baltic sprat and herring, and data from 1966-2012 for cod) with different levels of steepness were quite similar for both the B&H and Ricker S-R models (Fig. 1).Steepness varies between the species, and was equal to 0.71, 0.60, and 0.43 for cod, herring, and sprat, respectively, in the B&H model; it was 1.08, 0.56, and 0.39 in the case of the Ricker S-R relation.S-R models explain around 65% of the variance for cod, 40% for herring, and 30% for sprat.There were only two cases when the S-R model could not be fitted (cod +50% of h of B&H S-R relation, and sprat -50% of h of Ricker S-R relation).
The BRPs that were estimated by using different options for the input data were presented as proportions of the basic runs (the ones using the S-R relation and mean Mat, M, W, and Sel values from the full available data series for each species; Fig. 2).F MSY and its proxies demonstrated similar pattern.Sometimes, the highest and lowest observed changes in the input values were smaller or equal to 20%.The results show that the increase in values of Mat and W lead to an increase in BRPs, and for h to an increase in MSY, F MSY and its proxies.The opposite effect was observed in the case of M (stronger when the Ricker S-R relation was used) and Sel, that is, an increase in these variables led to a decline of the BRPs.The biggest differences in BPRs are visible for -70% of Sel for sprat, but only for F reference points; Sel does not have a big influence on MSY and B MSY when using both the B&H and Ricker S-R models.Additionally, the h of the Ricker S-R relation had a very small impact on B MSY .For the B&H S-R model, changes in B MSY were more visible in analyses with variable h, but for both S-R relations B MSY decreased while the other BRPs increased with increasing h.For the rest of the input data, the changes in all BRPs were not proportional, but had the same trend.
For cod and sprat, the best fits of both S-R relations were obtained for the full data set (correlation coefficient equal to 0.64 and 0.30, respectively).Cod showed a good fits also for data after the regime shift until 2012 (correlation coefficient of 0.43).For the sprat, it was possible to fit both S-R models only for the full data set and additionally for the period before the regime shift in the case of Ricker S-R model, and for the herring-for the full data set and after the regime shift until 2015 (correlation coefficient 0.48, slightly better than for the full data set).
In the case of cod, the period before the regime shift was characterized by a higher Ws and a lower Wc (Table 1).For herring and sprat, W and M were higher for the period before the regime shift and lower for the period thereafter (with one exception for sprat during 1995-2015, with the value slightly above 1).Selectivity has different patterns for both of the species; it is higher before the regime shift for herring, and lower for sprat.The value of h was lower after the regime shift period for herring, and higher for cod and sprat.The highest changes in the input data compared to the basic run were observed in W and h, and in the case of sprat for M as well.
The BRPs varied depending on the year range of the data used in the analysis (Fig. 3).For cod, the fishing mortality reference points were higher after the regime shift period for both of the S-R models (especially when 1995-2005 data were used), whereas MSY and B MSY were lower after the regime shift period.In the case of herring, F reference points estimated for data after the regime shift were lower than the basic run, and much higher for sprat.
GLM analyses were performed to investigate significant changes in the time series input data (Fig. 4).The year effects from the GLM were transformed into quartiles, as presented in Fig. 5. Dark grey indicates the

Table 1
The index of change of the model input data is defined as the ratio of the mean value of a given variable within the given time range to the mean value of that variable for the basic run for cod, herring, and sprat  For cod Ws, Wc, and Mat, three periods could be distinguished.First, there was a 'poor' period with a low fish condition until 1989 (significant changes in 1990 for Wc and Mat at levels of significance of 0.1 and 0.001, respectively, were detected).This was followed by a good period with high W and Mat during 1990-2006 (significant change was observed for Ws in 2007 and for Wc in 2009 at the level of 0.05).The last period was from 2007-2012, when low fish condition and Mat resumed.Cod Sel changed substantially over time, but high selectivity could be observed until 1975 and in the period 2000-2003, and low values in 2004-2009.For herring, a significant difference (at the level of 0.05) could be observed in W in 1986, and smaller values of M also started to appear from this period onwards.The value of M increased significantly in 2009 relative to the 20 previous years.According to changes in Sel, all data after 1983 are significantly different from the first year.In the case of W for sprat, one of the significant changes was observed in 1996.High values could be observed for M until 1989 and from 2008 onwards, with low values in the interim.The Sel data series changed in 1981 and in 2008, with levels of significance of 0.01 and 0.1, respectively.
Based on this analysis, two additional runs were performed to estimate BRPs for cod; one merged data from    1974-1989 and 1990-2007, with high and low M values, respectively.Both the S-R models that were fitted to the data from the periods defined above do not explain much of the recruitment variance.The best is the S-R relation for cod during 1966-1987, merged with 2007-2012 (correlation 0.45).For herring, the S-R model fitted best when the data from 1986 until 2015 were used, and slightly worse for 1986-2008 (with correlations of 0.34 and 0.29, respectively).In the case of sprat, the correlation for 1974-1995 was the highest (0.33), while the remaining periods did not exceed 0.1.
Due to the lack of fit of the both S-R relations for herring 1974-1985and sprat 1996-2015, and Ricker S-R relation for sprat 1974-1989and cod 1990-2006, the BRPs could not be estimated.The merged data for cod gave lower F reference points compared to estimates using 1990-2006 data (Fig. 6), and the trends in BRPs are similar to the runs that used the before and after regime shift approach.In the case of herring, the period that included the earliest data could not be analysed, as it had two data points less than the period before the regime shift (1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987), which did not allow the S-R models to be fitted.For sprat, the analysis shows that the BRPs are very sensitive, even to small changes in the data range.All BRPs were lower for 1974-1995 than for 1990-2007 (B&H S-R model), and fishing mortality reference points were higher for 1974-1989 than for 1990-2007 (Ricker S-R model).

DISSCUSION
For the B&H S-R model, an increase in h caused a decrease in the parameter a, which resulted in R approaching the asymptotic value for the low SSB.This means that the stock can produce high recruitment from the lower biomass, and can sustain higher F reference points and lower B MSY (e.g., cod 1996-2005 comparing to cod 1966-1987).For the B&H S-R relation, the maximum R-values are slightly decreasing with an increase in h (because of an increase in parameter b).For an increase in the h value of the Ricker S-R model (which results in an increase in parameters a and b), the peak of R is at a higher value; better recruitment ability in that case therefore also leads to the higher F reference points.Investigating the sensitivity of BRPs shows that, for the Ricker S-R relation, SSB is not as sensitive to changes in h as R, and the B MSY is therefore quite stable.
When the known data only reflects part of the recruitment curve, it is very hard to determine which type of model correctly defines the S-R relation (the Ricker or B&H type; Horbowy and Luzeńczyk 2012).However, in the analyses performed, the BRPs estimated by using different S-R models had rather similar trends.The h values for both models varied depending on the time series used, and it was difficult to establish the case in which the S-R relation was correctly reflected.The fits, which are far from the true S-R relation, could have resulted from using insufficient data, but also from having two sets of data showing different trends being treated as one data series (for example, the cod during 1966-2015).In some cases the S-R model could not be fitted, therefore BRPs were not estimated (e.g., herring from 1974 until the middle of the 1980s).Different approaches to get approximate S-R relations might be used (e.g., using estimates of h from other periods); however, to interpret changes in the BRPs estimated for differing environmental conditions, input data from the particular periods are needed.
The main difference between the species was that the period before regime shift was characterized by lower values of F reference points and higher B MSY and MSY in the case of cod, whereas the opposite trend was observed for sprat.The results could be explained by the higher mean biomass and recruitment of cod before the regime shift compared to the mean for the full data series (the same could also be observed for herring).Recruitment was relatively well correlated with an abundance of Pseudocalanus spp. in the spring for both species (cod and herring) and with salinity at 100 m depth in the Gotland Basin, which were both higher before the regime shift (Möllmann et al. 2003, 2008, Anonymous 2016b).The similarity in the biomass and recruitment trends for both It is in contradiction to the simulations of Holmgren et al. (2012), where the assumption of the planktonic food supply for herring reaching the levels present in the early 1980s caused large increases in F MSY and MSY.However, in this paper different S-R functions and models to estimate BRPs were applied.The almost complete disappearance of the predator cod from the northern Baltic proper (Casini et al. 2011) resulted in a stronger influence on the biomass of sprat in the northern areas, recruitment was correlated with the temperature at 60 m depth in spring and the sea surface temperature in summer in the Bornholm Basin, which was much lower before the end of the 1980s (Anonymous 2016b).This explains the opposite trend of the BRPs for sprat, with lower F reference points and higher B MSY and MSY after the 1990s (sprat 1990-2007 compared to 1974-1989 and 1974-1989).
The remaining input data, such as W and Mat, cause a rise in stock production; this results in a higher B MSY , and the stock can be fished with higher intensity (higher values of MSY and of F reference points), whereas the increasing M is reducing the ability of the stock to produce SSB; BRPs are therefore showing the opposite trend.For Sel, the interpretation is not as obvious.A decrease in Sel could be an indicator of low fishing pressure; on the other hand, it could result from a decrease in the amount of young fish in the catch, which could be the effect of a decrease in W.
Before the regime shift period, Ws were higher than the long-term mean for each of the three investigated species (Anonymous 2016b); after the regime shift, the situation was the other way round.The decrease in W in the case of herring and sprat was probably caused by a densitydependent effect (Horbowy 1997, Cardinale andArrhenius 2000) and the higher rate of food competition after 1990 (Casini et al. 2004).The sprat BRPs, as estimated by a model that accounts for density-dependent growth were higher compared to a model that assumes constant growth (Horbowy and Luzeńczyk 2016).In the case of herring, W has shifted within the last 30 years from being mainly driven by hydro-climatic forces (i.e., salinity, with high salinity increasing the abundance of the copepod Pseudocalanus spp., one of the main prey for herring) to an interspecific density-dependent control caused by an increase above the threshold biomass value for sprat, which is the main food competitor of herring (Casini et al. 2010).Sprat growth is dependent on temperature, and is based on the model of Frisk et al. (2015); the maximum body size is reduced with increasing temperature, because it affects respiration and activity costs.Pelagic fish production and high water temperature in autumn stimulated somatic growth of cod in 1980 (Uzars et al. 2000).The decrease in cod W after that period has a complicated background, and is closely connected with the worsening condition of the fish since 1990 (Anonymous 2016b).An even more drastic decrease in cod W in recent years is also visible in the BRPs when two of the periods (1996-2005 and 1996-2015) are compared.The F reference points are higher for the period, which does not include the most recent decrease in W.
In the case of the value of M for herring and sprat, the situation changed completely after the regime shift.It was caused by the relatively high biomass of cod until the end of 1980 (the biggest predator of sprat, and to some extent of herring) and the decrease of cod biomass, together with the M values of herring and sprat, from 1990 onwards (Anonymous 2013b).When the two periods for herring were reviewed (1986-2015 and 1986-2008, which exclude the most recent data when M was increasing again), the decreasing influence of higher M on F reference points is apparent.It is similar with the herring F MSY value, which substantially declined in the simulation, when cod stocks recovered (Holmgren et al. 2012).Similar results were obtained in the case of sprat when F reference points declined (approximately linearly) with the size of the cod stock (Horbowy and Luzeńczyk 2016).The M values for cod, and the values for Mat in the case of sprat and herring, were only available as constants for the full time series (Anonymous 2013b(Anonymous , 2016b)); the influence of these parameters on BRPs could therefore not be investigated.
When the two groups of input data, which diluted their influence on BRPs (e.g., high W, Mat, and h and low M and Sel), changed alternately for one species, then the estimates of BRPs could be quite similar.Moreover, fishing mortality reference points do not necessarily have to decrease when the environmental conditions are getting worse (e.g., cod before and after the regime shift).F is related to the ratio of the catch to the stock size.When, as in the case of cod, the stock size is more or less stable and Ws is decreasing because of environmental factors, then Sel can be decreased to maintain Wc at relatively the same value to prevent fishing on small and economically unprofitable fish.Consequently, a decrease in both components (catch and stock) could minimize the deteriorating effect of the environment on the F reference points.
Management decisions should consider the difficulties of achieving MSY simultaneously for multiple species and taking into account the trade-offs that result from interactions between species, mixed fisheries, and the multiple objectives of stakeholders (Kempf et al. 2016).For instance, including changes in W and M, caused by interaction between species, may significantly influence BRPs.The effect of omitting density-dependent growth in the estimation of MSY parameters, admittedly, leads to an underestimation of F MSY and MSY, but not including density-dependent natural mortality when it is an important driver of stock dynamics, causes overestimation of F MSY (Horbowy and Luzeńczyk 2016).Additionally, taking into account the regime shift while estimating MSY could prevent overfishing.Management plan during the rapid transition of environment should be even more cautious in the case of lack of knowledge of the ability of stocks to survive a number of years of adverse environmental conditions (Jiao 2009).To maintain the biodiversity, genetic diversity, and reduction of bycatch and waste, a more substantial reduction in fishing mortality may be necessary (Mace 2001).
F, and a and b are S-R model parameters that are established by fitting the recruitment equation for B&H S-R relation to the recruitment (R) and spawning stock biomass (SSB), estimated by using the stock assessment models(Anonymous 2013b(Anonymous  , 2016b)).A convenient way of parameterizing the S-R models is to use h (the concept of steepness,Mace and Doonan 1988), which is defined as the fraction of the recruitment from the unfished stock, R 0 , when the stock biomass is 20% of the unfished level, B 0 .Steepness characterizes the S-R relation; it is easier to compare these relations between stocks as opposed to using parameters a and b.The equations for a and b of the B&H S-R relation, including steepness

Fig. 1 .Fig. 2 .
Fig. 1.S-R relations fitted (black line) to the observed values -full available data series for each species (dots), and S-R relation with steepness decreased (dark grey) and increased (light grey) by 20% (solid line) and 50% (dashed line) in relation to steepness in the fitted S-R relation; 20% h line showed only if lacking 50% h; plots for cod (upper), herring (middle), and sprat (lower); for B&H (left panel) and Ricker (right panel) S-R relation; SSB in thousand tonnes and R in millions

Fig. 3 .
Fig.3.BRPs estimated in the runs for before and after the regime shift periods relative to BRPs estimated in the basic runs (with full available data series) for cod, herring, and sprat; estimated by using B&H (left) and Ricker (right) S-R relations

Fig. 4 .
Fig. 4. The year-effect from the general linear model (estimates and standard error) for input data used to derive the BRPs for cod, herring, and sprat

Fig. 5 .
Fig. 5. Estimates of the year-effects of the general linear model for different input data; time series transformed into quartiles; black represents the first quartile, dark grey the second quartile, lighter grey the third quartile, and light grey the fourth quartile for W and Mat, and the other way around for M and Sel, thus darker colours represent lower values of W and Mat and higher values of M and Sel; white represents no data or data assumed as constant

Fig. 6 .
Fig.6.BRPs estimated in the runs for specified periods (division based on the dynamics of the biological and fishery variables) relative to BRPs estimated in the basic runs (with full available data series) for cod, herring, and sprat; estimated by using B&H (left) and Ricker (right) S-R relations