Time-Varying Stock-Recruitment Model For Estimating Population Characteristics of Green Turtles

. The stock-recruitment relationship (SR), customarily used in fisheries assessment, can be used to analyze demographic data of sea turtles to infer changes in hatchling production ( R ) as a function of nester abundance ( S ), recruitment rates and the influence of environmental conditions on these population features. The SR Cushing model ( 𝑅 = 𝑎𝑆 𝑏 , where 𝑎 and 𝑏 are the model parameters) is well-suited for representing the dynamics of recovering populations, such as the green turtle ( Chelonia mydas ) in Campeche, Mexico. This study aimed to explore the SR Cushing model using a time series of the abundance of nesters and hatchlings (1984 – 2020). By applying local regressions (9-yr moving windows), we found that the time series of parameter b (the change in R as a function of S ) and the recruitment rate (hatchlings per nester) were inversely correlated with a 26-yr cycle of the Atlantic Multidecadal Oscillation – sea surface temperature (SST), over the Atlantic – ( 𝑟 2 = 0.83 and 𝑟 2 = 0.64 , respectively, at a 3-yr lag). Model diagnostics using the time-dependent Cushing model substantiated that the log-normal distribution of hatchlings of C. mydas in Campeche depends on the abundance of nesting females and on a low frequency SST signal ( 𝑟 2 = 0.98 ). The positive trend in nester numbers of green turtles in Campeche during the past 44 years may be the result of persistent conservation efforts, while the drastic and sporadic changes in the growth rate of annual arrivals and hatchling production are suggestive of population dynamics driven by low frequency, basin-wide environmental signals.


Introduction
The stock-recruitment relationship (SR) describes how the nonreproductive fraction of an animal population is determined by the abundance of adults (Iles and Beverton 1998). The recruits ( ) are the individuals who are about to enter the reproductive stock, the fishery or the general population, and the adults ( ) are the individuals who produced them. Particularly in marine turtles, hatchling numbers could constitute , and the abundance of nesting females . If the annual is plotted against the corresponding sizes of , through an ample range of stock sizes, it is possible to identify the shape of the SR of the population (Hilborn and Waters 1992). This concept has been extensively used in fisheries assessments, but seldom applied to determine key aspects of the population dynamics of sea turtles. This study focused on exploring the SR of the green turtle population, Chelonia mydas, of the Campeche coast, in the southern Gulf of Mexico.
In fisheries science, the SR has been modeled using a variety of curves, each presenting a different set of conditions. For example, the Beverton-Holt function (Beverton and Holt 1957) supposes that is constant at high levels of , while the Ricker equation (Ricker 1954) assumes that is maximized at high levels of , thereafter decreasing asymptotically. Since these two models contemplate the existence of a recruitment carrying capacity, an underlying assumption is that is measured across a wide range of values, if not all. However, when a population is recovering from low numbers, the SR is confined around small or intermediate values of and , and a power function may be a better-suited representation of such a situation (Illes 1994). The growing population of green turtles nesting along the Campeche coast, Mexico, and the Cushing SR model (Cushing 1971) possess these attributes.
The green turtle population in Campeche was depleted by poaching until the late 1970s, when the first Federal Government's chelonid conservation initiatives were implemented (Garduño-Andrade et al. 1999;Guzmán-Hernández et al. 2014); thereafter, it took more than 30 years for the population to begin recovering. During this elapsed time, the estimated number of nesting females increased from 6 in 1984 (1,267 hatchlings) to 2,728 (995,040 hatchlings) in 2017 (Guzmán-Hernández 2018). However, the growth rate of the nesting population and hatchling production have significantly changed. In fact, (VGH, PdML, MLC, Uribe-Martínez, Huerta-Rodríguez et al.), found that historical abundance data on local nesters and hatchlings show drastic and sporadic increments that surpass the steady population growth that is expected from a sustained conservation effort alone. Such increments divide the SR data of this population into three statistically distinguishable groups, each separated by an ~10-year gap: one from the mid-1980s through the mid-1990s, the second up to the mid-2000, and the last from that date until 2020. The authors do not give insights into the causes of these changes; however, Figure 1 is suggestive that they could be environmentally related.
As said, the SR data plot of the green turtle population nesting in Campeche shows no sign of a plateau, so it can be represented by a model that does not assume a carrying capacity effect, but still allows for estimating features such as the dependence of as a function of , and the annual recruitment rate (hatchlings per nester). Cushing's SR function is such a model, whose form is = , with and being the model parameters; in this case, the larger the value of , the faster grows as increases. Another advantage of the Cushing model is that environmental effects can be easily accommodated into its function, as has been done in other, more elaborate SR models (Marjomäki 2004;Olsen et al. 2011;Subbey et al. 2014). The present contribution aimed to explore the Cushing SR model using green turtle nesting females and hatchlings abundance data  from the Campeche coast, one of the largest green turtle populations in Southeast Mexico.

Fig. 1
Relationship between the (ln) number of nesters and (ln) the number of hatchlings of the green turtle (Chelonia mydas) from Campeche, in the southern Gulf of Mexico. The scatterplot shows three blocks of points discriminated by clustering analysis (K-means test). The time scale in the upper part of the graph represents the study period  and the circles are the positions of the elements of each block (differentiated by filling patterns) with respect to the calendar years. Behind the timeline in the upper part of the graph, depicted as a light gray area, a 26-yr cyclic component of the Atlantic Multidecadal Oscillation Index (sea surface temperature over the Atlantic Ocean in the Northern Hemisphere, AMO) is also displayed.

Stock-recruitment data
In this study, the number of females nesting along the Campeche coast was defined as the adult population or adult stock ( ) and the number of hatchlings they produced was assumed to be the recruitment ( ) for a given year, . Although these variables have been generated by two different methods, i.e., direct counts in the field  and indirect estimations from the nesting abundance (2013-2020), the results of both types of measurements have shown reliable compatibility. A detailed description of the methodology for data gathering and data treatment in this region can be found in Garduño-Andrade et al. (1999) and Guzmán-Hernández et al. (2015).
The time series of nester and hatchling abundances in Campeche are representative of other nesting sites along the Gulf of Mexico and the Caribbean ((VGH, PdML, MLC, Uribe-Martínez, Huerta-Rodríguez et al.), and have been generated uninterruptedly since 1984. In Campeche, green turtles nest more frequently at five beaches -where conservation camps have been installed and managed by the Federal Government. The camps are managed by the Comisión Nacional de Áreas Naturales Protegidas -for monitoring and conservation purposes: Chenkan, Sabancuy, Isla del Carmen, Isla Aguada and Cayo Arcas (Fig. 2). Hatchling and nester abundance data used in the present analysis are the sum total of these variables from all five camps.

Data exploration, modeling and fitting procedure
To investigate whether the parameters of the Cushing function were constant or varied during the 1984-2020 period, the model was fitted to the green turtle SR data in chronological order using 9-yr moving windows, assigning the estimated parameters to the year of the last data pair (the 9-yr window was chosen by balancing the number of observations for the regression and the statistical precision of each individual datapoint). Since and can be described equivalently by a logarithmic relationship in the Cushing function, the linearization ln = ln + ln was fitted to the first nine data pairs (1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992). Parameters and were estimated by least squares; then, the model was fitted to the next window (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993). This procedure was repeated until the function was applied to the last nine data pairs (2012-2020), resulting in a time series of estimated parameter values .
Standard linear regression also provides confidence intervals for each parameter. Since the local regressions are based on a constrained number of data points, these intervals may be wide, but they are nevertheless useful for exploratory analysis, as shown below. We tracked the value of from each fit because it is indicative of the type of dependence of with respect to . If > 1, increases "exponentially"; if = 1, increases linearly; and when < 1, increases "logarithmically". In other words, the greater the value of is, the less hatchling production is restricted by the abundance of nesting females.
By analyzing the time elapsed between the three data groups in the SR of the green turtle in Campeche found by (VGH, PdML, MLC, Uribe-Martínez, Huerta-Rodríguez et al.; Fig. 1), we propose that an environmental signal with a periodicity of between 20 or 30 years resembles this spacing. In other Chelonid species that share reproductive season with green turtles, such as the hawksbill turtle in the southern Gulf of Mexico, basin-wide oscillations in the sea surface temperature (SST) seem to affect the abundance of nesting populations (del Monte-Luna et al. 2012). It is then reasonable to suggest that parameter of the Cushing model applied to the SR data of the green turtle in Campeche, if shown to be time-variant, may also be related to an environmental covariable. With this empirical motivational framework in mind, we chose a specific SST cycle of ~26 years identified by del Monte-Luna, et al. (2015) in the AMO (Atlantic Multidecadal Oscillation: SST over the Atlantic Ocean in the Northern Hemisphere; Fig. 3), and examined these relationships with the rigor of statistical modeling.
The main postulated concept is that and are time-dependent, possibly through a 26yr SST cycle. The proposed model is: where has a normal distribution with mean 0 and (constant) variance 2 and SST is the 26-year cycle of the AMO in year . This is equivalent to positing that = exp( 0 + 1 SST ) and = 0 + 1 SST , where the (unknown) parameters 0 , 1 , 0 , 1 regulate the way in which the Cushing model varies as a function of the temporal SST. The possibility that SST bears no relationship with is also accounted for, by setting 1 = 1 = 0. In Model (1), is the median value of for a given and SST, which has been randomly perturbed by a multiplicative, log-normally-distributed factor. It is relevant to note that this stochastic structure implicitly accounts for heterogeneous (increasing) variance in as increases, as well as positive skewness. Both of these characteristics are apparent in the data (Fig. 4). The recruitment rate (number of hatchlings per unit nester, RR), obtained by differentiating med( ) as a function of , is: The parametrization in Equation (1) allows the model to be rewritten in an equivalent linear form, which involves log-normal distributions and has occurred in discussions (Iles 1994) regarding ways of incorporating randomness into SR relationships: Equation (4) was fitted by ordinary least squares (equivalent to maximum likelihood under the normality assumption), using all available observations ( , ) jointly without any grouping due to moving windows. This provided estimates for the model parameters in (4) as well as standard errors and confidence intervals. Estimates of and and RR as a function of SST and are derived by substituting estimates in (1) and (3). We initially postulated the linearization in Equation (4) merely for numerical convenience, giving rise to standard linear least-squares. However, diagnostic regression procedures (see Iles, 1994) based on the analysis of residuals subsequently proved that the independent, constant variance and additive normality assumptions implicitly stated in (4) are tenable for the data.

Results
The values of parameter that resulted from applying the Cushing model as local regressions to the green turtle SR data were empirically found to be inversely and maximally correlated to the 26-yr cycle of the AMO at a 3-yr lag ( 2 = 0.83; < 0.05; 2 is the coefficient of determination, Fig. 3). The corresponding confidence intervals are wide as anticipated, but the mere fact that some of them do not overlap in Figure 3 indicates that this parameter is not constant over the observed range.
An interpretation is readily visualized in Figure 3: regardless the nester abundance, as more females go out to nest during warming periods, total hatchling production decreases rapidly ( = 0.6). In contrast, during cooling periods, hatchling production is significantly higher in the presence of more nesting females ( = 1.2). Although parameter a holds no biological meaning, we found no evidence suggesting it has varied over time with SST. The Cushing model presented in Equation (4), extended by incorporating the 26-yr cycle in the AMO with a 3-yr lag, showed an acceptable performance, as the coefficient of determination relating the observed and estimated was 2 = 0.98 ( < 2.2e-16; Fig. 4). Table 1 (all model coefficients, except the one associated with parameter , were statistically significant at = 0.05). The estimate for the variance parameter ( 2 ) was 0.0678. The time series of RR values, derived from Equation (3), were also negatively correlated with the SST at a 3-yr lag, changing roughly twofold between the extremes of the AMO cycle, from 349 (in 2017) to 194 (in 2008) hatchlings per nester as a function of . This RR, as well as the parametrically estimated value of parameter are also superimposed in Figure 3. Figure  4 shows the observed values of plotted against 97.5% and 2.5% estimated quantiles obtained from fitting the regression model; the shaded regions in Figure 4 are thus the 95% prediction intervals for each year's .

Discussion
In this study we extended the use of a time-varying version of the Cushing SR model from fisheries assessment to the population dynamics of the green turtle in Campeche, assuming the nesting females were representative of the adult stock and the hatchlings were recruitment. The interpretation of the model parameters and the derived quantities must strictly be circumscribed by these two components of the population and within the observed range of and . The status of any other life stages, and the behavior of this SR outside those ranges should not be inferred from the functions used.

Fig. 3 Time series of (thin black line) a 26-yr cyclic component of the Atlantic
Multidecadal Oscillation Index (sea surface temperature over the Atlantic Ocean in the Northern Hemisphere, AMO); parameter of the Cushing stock-recruitment model (thick gray line), fitted by local regressions (9-yr windows) using hatchling and nester abundance data  of the green turtle (Chelonia mydas) population from Campeche, Mexico; and 95% confidence intervals of the moving-window parameter (light gray shaded area); estimated parameter of an extended Cushing stockrecruitment model (dashed line), explicitly incorporating the AMO 26-yr cycle into its specification, and the estimated recruitment rate (number of hatchlings per nesting female, empty circles) derived from the extended Cushing model. Time series of parameter from local regressions, and its 95% confidence intervals, spanning from 1992 to 2020 because the estimates of the first regression window (1984 -1992) were assigned to the last data pair of that period.

The SR from fisheries to sea turtles
Originally, the experimental evidence of SR originated from laboratory testing with terrestrial beetles (Watt 1955). Since then, its application has focused on marine fish and invertebrates, with little attention given to long-lived species. Particularly in sea turtles, Doi et al. (1992) applied the software DOIRAP, an age structured model complemented with a Beverton-Holt function, to estimate the maximum sustainable yield of Cuban hawksbill populations (Eretmochelys imbricata). The limitations of Doi's model are that recruitment was estimated under equilibrium conditions and survival and growth rates were assumed to be constants (Crouse 1999;Heppell and Crowder 1995). Therefore, two critical issues emerge when applying SR models to sea turtle demographic data: how recruitment is estimated and the assumptions of the equilibrium conditions (time-invariant parameters).
In fisheries science, "true" recruitment may be conceived as individuals who just entered the reproductive stock (Subbey et al. 2014). However, in this case, this information is difficult to obtain because it implies the existence of (1) sufficient data for estimating the abundance of each age class of the population and (2) the time series of the abundance of each age class, spanning at least the number of years of the age at first maturity plus the number of years (data pairs) over which the model will be fitted; for example, ten years. In green turtles, this means having an ~45-year-long dataset of the abundance of each age class, to yield an effective 10-year SR time series (assuming ~35 years as the age at first maturity; Zug et al. 2002), which makes fitting an SR curve almost a prohibitive task. In this case, we redefined as the hatchlings (new individuals added to the general population), which, along with the nesters, are the population segments easiest to count.
We do not suggest that hatchling abundance represents the number of individuals entering into the reproductive population. In sea turtles, the natural mortality between egg and adulthood is the highest among all life stages (Heppell et al. 2003); therefore, the strength of the age class that constitutes the "true" recruitment may not be proportional to that of the hatchling cohorts that produced it. Alternatively, we used demographic variables typically obtained in the field as direct inputs into a simple model for inferring basic population features, as recommended by Piacenza et al. (2016). Moreover, the more unrealistic assumptions of fishery models previously applied to sea turtles were overcome. First, by not having to estimate the "true" recruitment and instead reframing a stock-recruitment model exclusively in terms of the relationship between nesters and hatchlings; second, showing that parameter of the Cushing model is time-variant; and, third, enriching the SR model by setting its parameters as the functions of an environmental covariable.

Fig. 4
Observed number of hatchlings per year, with 95% prediction intervals constructed from fitting the extended Cushing model (Equation 4). The shaded regions outline 97.5% and 2.5% estimated quantiles of the log-normal distribution associated with for a given sea surface temperature (from a 26-yr cycle of the Atlantic Multidecadal Oscillation Index) and for each year. This evidence shows that variance increases as the median (or mean) increases. The inner graph shows how the model fits to the smallest values of the dataset.

Time-varying parameters and recruitment rate
The results indicated that as SST increases in a multidecadal scale over the North Atlantic Ocean, the recruitment rate and the total hatchling production of green turtles in Campeche decrease when more females go ashore to nest, in contrast to what happens during a cooling period. This finding is in line with the relationship that has been identified between hatchling success, nester abundance and environmental conditions (Broderick et al. 2001;Chaloupka et al. 2008). The pattern is that during a warming (cooling) period at the Yucatán Peninsula, the ecosystem productivity in critical habitats decreases (increases; Manzano-Sarabia and Salinas-Zavala 2008), thus impoverishing (enhancing) the quality and quantity of available food for green turtles. This environmental setting tends to be unfavorable (favorable) for the body condition and remigration times of sea turtles (Stubbs et al., 2020), causing a reduction (increment) in the recruitment rate -clutches laid-, hatchling production -hatchling success-, and annual arrivals of nesting females (Solow et al. 2002;Saragoça et al. 2020).
Decadal signals consistent with those found in this study have also been identified in the nester body length of green turtles in Hawaii, USA (Piacenza et al. 2016), as well as in nest counts of loggerheads (Caretta caretta) in Florida (FWRI 2021). However, we acknowledge that in the green turtle population of Campeche, the 3-yr lag in the correlation between the AMO cycle and the time series formed by parameter , and between the AMO and the recruitment rate, suggest that SST is not causative. Instead, all of these variables may respond to a large-scale, yet unidentified, environmental factor with a periodicity of ~26 years. The identification of decadal signals in demographic/autoecological variables of sea turtles and environmental indices, represents fertile ground for new research avenues.

Trends in abundance of hatchlings and nesters
We resorted to extensions of the Cushing model akin to Marjomäki's (2004) and Subbey's et al. (2014) in the sense of incorporating functions of covariates (environmental factors) that alter the SR relationship. However, this is done here by allowing each individual parameter in an SR model to depend on the covariates, instead of adding (Marjomäki 2004) or multiplying (Subbey et al. 2014) the whole SR function by a single term. The fitting procedure of the enriched Cushing model showed that the median of the log-normal probabilistic distribution of the annual number of green turtle hatchlings in Campeche, depended on the abundance of nesting females and a lowfrequency, basin-wide SST signal. It is not surprising that the annual number of nesters was a reliable predictor of hatchling production (Fig. 1), as has been demonstrated in other instances (Piacenza et al. 2016). Notably, by means of a dynamic modeling framework, it was possible to unveil a multidecadal signal in basic population features concealed behind a time series of demographic metrics showing an exponential trend and a high interannual variation.
The overall abundance trend of green turtle hatchlings and nesters in Campeche, may reflect the effects of sustained conservation efforts. In Southeast Mexico, since the 1980s, logistics and human resources for surveilling and protecting sea turtles have become increasingly efficient (Guzmán-Hernández et al. 2014. This has resulted in an exponential increase in the pooled regional turtle abundance, which is roughly akin to the abundance trend of individual species (Garduño-Andrade et al. 1999).
Although translocation activities of sea turtle nests officially stopped in 2012 (the management switched to mostly in situ monitoring), there were no significant effects in terms of nests lost (Guzmán-Hernández et al. 2015). Currently, green turtle numbers at Campeche suggest that the regional nesting population is successfully recuperating (VGH, PdML, MLC, Uribe-Martínez, Huerta-Rodríguez et al.), as are other populations of the species in distant regions of the world (Broderick et al. 2006).

Methodological generality
A byproduct of the statistical model is an explicit distribution for given and SST . It is log-normal with parameters 0 + 1 SST + 0 ln( ) + 1 SST ln( ) and 2 . This distribution is right-skewed, having median med( ) = exp( 0 + 1 SST ) 0 + 1 SST , mean given by med( )exp( 2 /2), and variance [med( )] 2 [exp( 2 ) − 1] [exp( 2 )] (notice that the variance is proportional to the square of its median, confirming that the variance increases as becomes larger). This distribution is potentially useful for implementing simulation studies involving and for evaluating scenarios based on simulated, plausible or estimated values and SST in future years. By taking quintiles, the distribution enables the construction of prediction intervals for for a given and SST, as shown in Figure 4. Probabilistic statements regarding also become accessible knowing the distribution. For example, the probability that will exceed (or fall short of) a given threshold of interest for conservation.
Our approach was motivated and driven by the green turtle in Campeche. However, the modeling procedure has wide-spread general application. It consists of two stages: (1) Begin with a tentative SR relationship and estimate its parameters locally as a function of time. To examine whether any patterns emerge relating the behavior of these estimated parameters with one or more time-stamped environmental variables, (2) let the parameters in the SR relationship individually depend parametrically on covariates (using simple linear functions, as we have done, is worthwhile as a first try). We formulated and fitted a regression model; our case was conveniently reduced to multiple linear regression, but it may turn out to be a nonlinear model (e.g. an equation from a generalized additive model). In any case, check for the validity of the statistical assumptions in the model to ensure that any calculated confidence intervals and hypothesis tests regarding the parameters are ultimately valid.