Modeling the Effect of Environmental Factors on the Ricker Stock-Recruitment Relationship for North Pacific Albacore Using Generalized Additive Models

This study investigates how the inclusion of El Niño episodes affects the modeling of recruitment in North Pacific albacore. The relationship between spawning stock size and recruitment for the 1976 2004 period, including environmental variables (Southern Oscillation Index and sea surface temperature anomalies), was conducted by general linear and generalized additive models (GAM). The results indicate that the Southern Oscillation Index and the interaction between spawning stock size and sea surface temperature anomalies have significant effects on the recruitment of North Pacific albacore. GAM fit the original data better than general linear models, according to the Akaike information criterion. The recruitment declined when El Niño episodes occurred, but increased when La Niño episodes occurred. The recruitment also increased when the spawning stock size density was above 11.21 million tons with a -0.3°C sea surface temperature anomaly. El Niño was observed to have both long-period (above 10 15 years) and short-period (1 2 years) effects on North Pacific recruitment.


INTRODUCTION
Albacore tuna (Thunnus alalunga) is a highly migratory species that is distributed mostly in the temperate waters of the Pacific, Indian and Atlantic Oceans and the Mediterranean Sea.In the North Pacific Ocean, adult albacore are predominantly found along the outer margin of the Kuroshio Current, the North Pacific Transition Zone (NPTZ) and the California Current.Their spawning grounds are in tropical and subtropical waters (Ueyanagi 1969;Otsu and Sumida 1970;Nishikawa et al. 1978).Generally, mature albacore migrate annually from the east to the west Pacific in an anticlockwise direction.The immature albacore have a similar migration route, extending along the Kuroshio Current from 25° -35°N and from 130° -180°E (Kimura et al. 1997).Due to their migration pattern, albacore have been exploited by fishermen from Japan, Taiwan, Canada and the U.S. since the 1950s.As a result of this long history of exploitation, the biomass of the North Pacific albacore stock has been in decline (Sibert et al. 2006).The population declined steadily from 1975 to the late 1980s, but recovered to roughly 1.2 million tons in the 2000s.Despite the recent increase, however, the spawning stock biomass (SSB) has shown fluctuations since the 1970s (Takeuchi 2004), and the recruitment level is not well known.The recent increase in abundance was assumed, without any affirmative evidence, to be the result of a strong 2001 year-class recruits (Crone and Conser 2002) which probably cannot last for more than 10 years.
Several studies (e.g., Fogarty et al. 1991;Santiago 1998;Fromentin and Restrepo 2001) indicated that the mortality rates of larvae and juveniles are strongly affected by temperature and food availability.The population abundance and the number of recruits are related to the abundance of larvae and juveniles and the environmental variability.However, the stock recruitment relationship is not well known, and the link between recruitment and environmental parameters is also poorly understood and has not been modeled.Studies on environmental variations related to population abundance (Myers 1998) and recruitment success (e.g., Planque and Frédou 1999;Begg and Marteinsdottir 2002;Nishida et al. 2007) have became much more common in recent years.Among those studies, sea temperature (Planque and Frédou 1999;Begg and Marteinsdottir 2002;Nishida et al. 2007) and large-scale climatic oscillations (Yatsu et al. 2005;Arregui et al. 2006) were frequently investigated and discussed.
In the case of South Pacific albacore stock, the recruitment of young albacore was found to fluctuate in accordance with the El Niño/Southern Oscillation (ENSO) strength (Fournier et al. 1998), in which the recruitment significantly declined during El Niño episodes, but recovered during La Niña phases.However, there have been few attempts to assess the environmental influence on the recruitment of North Pacific albacore.
To examine the stock-environment-recruitment relationship properly, various models, such as linear regression models (Sparholt 1996), generalized additive models (GAM) (Daskalov 1999;Megrey et al. 2005), principal component analysis (PCA) (Yatsu et al. 2005), artificial neural networks (ANN) (Arregui et al. 2006) etc. have been applied.The most important regression technique for this purpose is the non-parametric GAM.Because of its flexible construction, the GAM can smooth explanatory variables separately and identify non-linear and non-monotonic relationships (Guisan and Zimmerman 2000;Guisan et al. 2002;Thuiller 2003).Although GAM is rarely used in recruitment modeling today, it has been widely applied in studies of species-environment interactions, especially when little is known about the relationships between factors and observations (Guisan et al. 2002;Olivier and Wotherspoon 2005).To clear the use of GAM in recruitment modeling, in this paper we apply different statistical methods (GLM and GAM) to build a stock-recruitment model incorporating some possible environmental factors for the North Pacific albacore stock.

Data Source
The study area in the North Pacific Ocean is defined in Fig. 1.Recruits of 1-year-old albacore (Fig. 2) were assumed and adopted from the Report of the Nineteenth North Pacific Albacore Workshop (Stocker 2004).Recruits from 1976 -2004 were estimated by adaptive virtual population analyses (VPA) in which catch-at-age data were summed from various fisheries and the CPUE was used to tune VPA (Stocker 2004).The SSB (S) (Fig. 2) was also estimated from mature oogive product [the proportion of sexually mature individuals at age (m a, x )] and the abundance (N) at age a and sex x.The proportions of sexually mature North Pacific albacore were proposed by Hsu and Chen (2005) as follows: where L is the mean fork length of each length class interval in cm, calculated by von Bertalanffy's growth equation proposed by Clemens (1961): . L e 135 6 1 and Chen 2005) are equivalent to age 4 for both males and females.The weight (W) in kilograms for all age groups can be estimated by the weight-length equation W = 4.936 × 10 -5 L 2.98 (Clemens 1961).A number of studies have shown differences in the length and weight of male and female albacore over the age of sexual maturity (e.g., Santiago 1998;Chen et al. 2010).
Worldwide monthly sea surface temperature (SST) data were collected from the Hadley Center Global Sea Ice and Sea Surface Temperature (HadISST) data set (Rayner et al. 2006).The annual mean SST in the North Pacific Ocean for each year was recalculated and compared with the mean SST from 1975 -2004.The Southern Oscillation Index (SOI), which is the difference in sea level pressure between Tahiti and Darwin, is an effective indicator of El Niño events in the North Pacific Ocean (Niebauer 1988).The SOI data obtained from the National Oceanic and Atmospheric Administration Climate Prediction Center (NOAA/CPC) were also transformed into anomalies as SST.

Statistical Analysis and Modeling
The well-known Ricker model (Ricker 1954(Ricker , 1975) ) was applied to assess the relationships between recruits, SSB and environmental conditions, because among various stock-recruit models the Ricker model is the best fitting model (Simmonds et al. 2011).Ricker models have been modified and widely applied in previous analyses of stock-recruitment relationships (Sakuramoto 2005;Yatsu et al. 2005;Arregui et al. 2006;Nishida et al. 2007).Ricker models are based on the assumption that the mortality rate of juveniles depends on the initial cohort size.A fundamental and extended multiple linear regression family of Ricker models is generally given as follows: R Se in which the function of f(X) is to estimate the effects of the environment on different factors, X s ; and ε t is the error structure of the relationship with a lognormal distribution.Therefore, the stock-recruitment relationship with environmental variables included can be expressed as where X t, 1 , X t, 2 , …, X t, p , are p's environmental variables that could influence the stock-recruitment relationship in year t, and γ 1 , γ 2 , …, γ p are the corresponding parameters for the environmental variables.Although multiple linear regression facilitates the incorporation of new variables, it may lead to strong collinearity because of variable S, the spawning stock size.Eberhardt (1970) provided a diagnostic method for the correlation between ln(R t+1 /S t ) and S t .The equation is where Y t is equal to ln(R During the analysis, the environmental factors that influence recruitment were examined by conducting a deviance analysis using a stepwise regression procedure, in which each step was performed by both the generalized linear model (GLM) and the generalized additive model (GAM).
A GLM is a more flexible extension of a regression model.A GLM is constructed by considering the mean of the response variable through a link function together with a simple combination of variables.A GAM (Hastie and Tibshirani 1990) is a non-parametric extension of a GLM.The statistical properties of GLMs are preserved in GAMs; moreover, by using a smoothing function instead of a link function, more flexible non-linear associations between the explanatory variables and the response can be estimated with the appropriate data.Thus, the general GAM form can be represented as where g(μ) is the expected response variable, and f j is the smoothing function for explanatory variable X j .In this analysis, we use a locally weighted regression called loess (Cleveland 1979), with cubic splines (Hamming 1973) and B-splines (de Boor 1978) as smoothers for estimating the stock-environment-recruitment relationships.The Akaike information criterion (AIC) (Akaike 1974) was used to select the best sub-model from the full set of models.The model with the smallest (most negative) AIC value indicates the best fitting model for the data.All of the statistical analyses were performed using R software (R Development Core Team 2011).

RESULTS
The different time series patterns of recruitment, spawning stock size, SOI and SST anomaly are compared in Fig. 2. The recruitment varied between years but followed a general trend.In 1976 -1990, the recruitment fell off gradu-ally.After 1990, the recruitment recovered but with severe variation between years.Comparing the recruitment with the SOI, there was a clear decline in El Niño years (SOI < 0) such as 1982 and 1997, but a recovery in La Niña years (SOI > 0) such as 1981 and 1998.The major SST anomalies after 1990 were positive and showed that the SST was slightly higher than before.However, there was no direct relationship between recruitment and SST anomalies.The spawning stock did not show a distinct association with either recruitment or SST anomalies.
The Ricker spawner-recruitment curve is shown in Fig. 3.The different recruitments estimated by the four processes are compared with a primitive data set in Fig. 4.Among the density-dependent assumptions, the traditional Ricker model (R/GLM) of North Pacific albacore in relation to spawning stock size was R Se . .
(R 2 = 0.01, P = 0.606).The maximum recruitment threshold of the R/GLM was 52.78 million generated by a spawning stock size of 55.56 million.The R/GAM recruitment estimation presented the peaks of original recruitment, such as in 1977 and 1982, better than the R/GLM (Fig. 4).However, the R/GAM still failed to depict the dramatic variation after 1993 and produced a steady curve between peaks and troughs (Fig. 4).On the temporal scale, unlike the primitive recruitment data, the recruitment predictions estimated by the R/GLM remained steady at around 20 -30 million from 1976 -2004 (Fig. 4).Under the GAM process, the best fitting Ricker model constructed by cubic splines (R/GAM) showed better flexibility for the estimation.The smoothing curve is shown in Fig. 5 (P = 0.046).Values greater than zero on the y-axis indicate a positive influence on the recruitment estimation, whereas values less than zero indicate a negative influence.The spawning stock effect indicated that recruitment was higher when the spawning stock was around 12 million and decreased when it was around 10 million.The quantitative results for the influence of environmental variables are presented in Table 1.The analysis of deviance table indicated that the recruitment of North Pacific albacore was significantly influenced by SOI and the interaction between spawning stock size and SST anomalies (p < 0.05).The interaction between spawning stock size and SST anomalies is illustrated graphically in Fig. 6, which shows that the stock-recruitment relationship gradually changed as the SST anomaly increased.When the spawning stock was greater than 11.21 million the recruitment grew rapidly with the spawning stock when the SST anomaly was -0.3.As the SST anomaly increased the recruitment growth rate declined and hardly recovered once the SST anomaly increased to 0.2.In addition, as the SSR anomaly increased to 0.2 the S-R curve looks more like the Ricker curve.However, the opposite pattern was found when the spawning stock was below 11.21 million.
The Ricker model including SOI and the interaction between S and SST (RE/GLM) can be presented as R Se .After combining the environmental influences, the recruitment estimated by the RE/GLM fit the original recruitment data better than the R/GLM and the R/GAM and was able to account for the substantial variation after 1993.The same factors were included and constructed using the GAM (RE/GAM).The best fitting model was smoothed by cubic spline and loess.The RE/GAM smooth curve (Fig. 7) shows a negative effect on recruitment when SOI was negative, but the effect shifted to positive when SOI increased.In terms of the interaction between S and SST anomalies, the smooth curve indicates that there was a slight decline when SST anomalies were greater than zero and the spawning stock size was around 12 million.The RE/GAM prediction was similar to that of the RE/GLM but provided a better fit to the original recruitment pattern.
Table 2 presents the four models and the AIC values for the different modeling methods.The AIC of the R/GLM was 31.07.After constructing the GAM, the AIC of the R/GAM was reduced to 28.74.Similar results were found after the SOI and SST anomaly effects were included.The AIC of the RE/GLM was 28.59 and the AIC of the RE/GAM was reduced to 24.46 using the GAM.The results demonstrate that the flexible GAM produces a better fit to the original data and thus improves understanding of the recruitment mechanisms with environmental effects.

DISCUSSION
The spawning size and recruitment of North Pacific albacore was in this study adopted from a previous stock assessment using VPA (Stocker 2004).The accurate spawning size and recruitment estimation by VPA was affected by age composition, natural mortality and fishing mortality of old ages (Hildén 1988;King 2007).However, these factors  were not considered in this study which may produce some bias in the results.
Adult North Pacific albacore spawn in March and April (Chen et al. 2010).The spawning area extends from the western North Pacific Ocean to the Hawaiian Islands between about 15° -25°N (Nishikawa et al. 1978).After hatching, the young albacore gradually move to a higher latitude (25° -40°N) and recruit to the stock.Environmental influences might be one of the key factors affecting the egg or larval survival rates of marine fish (Pepin 1991).In our case, these effects were demonstrated by the unobvious relationship between spawning stock size and recruitment (Table 1).ENSO episodes are a crucial factor in the recruitment of albacore and caused the size of recruitment to decline in our study.A similar result was also observed in South Pacific albacores by Lehodey et al. (2003).However, different from their study, we found that the El Niño effects on North Pacific albacore recruitment were not only long-term (above 10 -15 years) but also occurred over a short period (1 -2 years).
Environmental influences might be one of the key factors affecting the egg or larval survival rates of marine fish.In our case, these effects were demonstrated by the unobvious relationship between spawning stock size and recruitment (Table 1).A similar result was also observed in South Pacific albacores by Lehodey et al. (2003).However, different from their study, we found that the effects of El Niño on North Pacific albacore recruitment were not only long-term (above 10 -15 years) but also occurred over a short period (1 -2 years).
ENSO episodes are widely known to cause serious effects on various aspects of the marine ecosystem; for example, a change in primary production can lead directly to a decline in pelagic species populations.Previous studies have noted that primary productivity in the equatorial Pacific (10°S -10°N) decreases during El Niño periods and recovers during La Niña periods due to the reduction in  nutrients (Blanchot et al. 1992;Chavez et al. 1999).Abnormal variations in chlorophyll-a concentrations, such as the movement of the transition-zone chlorophyll front in El Niño and La Niña periods have also been observed in the subarctic-subtropical transition zone (30° -45°N) (Polovina et al. 2001;Sasaoka et al. 2002).There is still insufficient information on the spawning area of albacore in the North Pacific (15° -25°N) to assess the effect on primary productivity during El Niño periods.However, the influence of primary productivity variations could lead to a decrease in the fecundity and survival rate of eggs and larvae, which could explain the recruitment decline in El Niño years.
Several studies have noted that the dynamics of albacore recruitment may be influenced by SST fluctuations (Uosaki and Uehara 2002).Temperature changes may have direct effects on the survival rate of eggs and larvae of albacore.Pepin (1991) reported that an increase in temperature leads to more rapid daily development rates in the early life stages of marine fishes, but also increases daily mortality rates during both egg and postlarval stages.However, in our study, the influence on recruitment was represented as an interaction between SST and spawning stock size.As shown in Fig. 6, the relationship between recruitment, spawning stock size and SST anomaly revealed opposite results for high and low spawning stock size.A possible explanation for such a discrepant response is the food abundance variation under different temperatures.This supposition might be a new topic to be explored and more detailed studies are necessary.
The GAM is rarely used in stock-recruitment research, but is increasingly applied in spatial abundance estimations.The technique is highly flexible for fitting and estimating the variation appropriately when there is a lack of information (Daskalov 1999;Piet 2002;Venables and Dichmont 2004).Megrey et al. (2005) compared various modeling methods and concluded that the GAM is useful for identifying the relationships between recruitment and influential factors and can be applied first to the predictor and response variables to ascertain the functional form empirically from the data.In our results, the performance of the Ricker stock-recruitment model was enhanced, with better flexibility in the smoothing functions in the preliminary examination of recruitment dynamics for North Pacific albacore.Stock-recruitment studies are often limited by insufficient time series data and information about influential factors.To overcome the data limitations a pilot study conducted by GAM could be a possible solution.
The environmental influence on the stock-recruitment model is a central point of contention.Some researchers consider that the environmental effects are overemphasized ( 1998), whereas others argue the opposite.A review by Needle (2002) suggested that including environmental influences on recruitment in recruitment estimation is beneficial to the integrity of stock-recruitment models.In our study, a similar conclusion can be drawn from Fig. 4, which showed a great improvement after environmental factors were included.To estimate the recruitment accurately and to design a suitable management strategy, it is necessary to incorporate environmental factors into stock-recruitment models and more attention should be paid to this topic.
It is known the SST may affect marine organisms including fish on distribution and mortality (Hsu 1987;Loher and Armstrong 2005;Hsieh et al. 2009).Albacore is a highly migratory species.The juveniles usually move northerly at temperate waters to feed.The mature spawners migrate southerly to warmer waters to spawn.Because ages at the catches (or catch at size) of albacore ranged from around 1 year old to almost 9 years old, which those fish caught distribute almost all waters as Fig. 1 indicated.Moreover, the hatched region indicated in Fig. 1, inferred from market sampling survey for reproduction study (Chen et al. 2010), is the only known spawning region where the albacore resided include almost all sizes.As a consequence, to clearly identify where albacore are staying with exactly locating the environmental factors, mainly SST seems difficult without pursuing long term survey.This is an encouraging issue for further study on the relationship between stock recruitment and climate changes.

CONCLUSION
The results obtained in this study could in part explain the relationship between recruitment and spawning stock size in relation to oceanic environment factors such as ENSO episodes.The information may be used to improve the stock assessment in understanding the spawnerrecruitment relationship of North Pacific albacore.The SST anomaly may not be the only environment factor influencing fish recruitment.The current study may encourage the investigation of fish populations in relation to oceanic environmental changes by incorporating density-dependent factors into the relationship, as in the Ricker model in the present study.However, as the input data were mainly from 2004, the results obtained in this study may not represent the current population status.
Fig. 1.Study area.The black line shows the defined boundary of the North Pacific Ocean.The shadow area indicates the spawning location of North Pacific albacore.

Fig. 2 .
Fig. 2. The time series pattern of recruitment, spawning stock, Southern Oscillation Index (SOI) and Sea Surface Temperature (SST) anomalies.A negative SOI indicates the occurrence of El Niño, whereas the positive SOI indicates La Niña.The solid lines are recruitment in year t + 1 and the bars are SOI, SST anomalies or spawning stock size.

Fig. 3 .
Fig. 3. Ricker stock-recruitment curve of North Pacific albacore.The scatterplot shows the relationship between annul recruitment (millions) and spawning stock (millions).

Fig. 4 .
Fig. 4. Comparison between the original recruitment and recruitment estimated by the Ricker model (R/GLM), Ricker model constructed by GAM (R/GAM), environmental Ricker model (RE/GLM) and environmental Ricker model constructed by GAM (RE/GAM).

Fig. 5 .
Fig. 5. Smooth spline function of spawning stock used in the R/GAM.The x-axis indicates the spawning stock of North Pacific albacore, and the y-axis indicates the effect on recruitment.

Fig. 6 .
Fig. 6.The interaction between spawning stock size and sea surface temperature (SST) anomalies.The x-axis indicates the spawning stock of North Pacific albacore, and the y-axis indicates the recruitment.The different lines present the spawner-recruitment relationship with increased SST anomalies.

Fig. 7 .
Fig. 7. Smooth spline function of the Southern Oscillation Index (SOI) and the interaction between spawning stock size and sea surface temperature (SST) anomalies used in the RE/GAM.The SOI was smoothed by loess, and the spawning stock and SST anomalies were smoothed by cubic splines.The x-axis indicates the SOI, SST anomalies and spawning stock of North Pacific albacore separately, and the y-axis indicates the effect of these factors on recruitment.The 95% confidence intervals (dashed lines) and fitted values are also shown in the figure.

Table 1 .
The deviance table analysis for the environmental Ricker model (RE/GLM) showing the change in deviance after spawning stock size (S), Southern Oscillation Index (SOI) and sea surface temperature anomalies (SST) were introduced sequentially.