journal article
LitStream Collection
Model-free time series analysis detected the contributions of middle-age spawner biomass and the environment on Pacific bluefin tuna recruitment
Nakayama,, Shin-Ichiro;Fukuda,, Hiromu;Nakatsuka,, Shuya
doi: 10.1093/icesjms/fsz129pmid: N/A
Abstract The relationship between the biomass of Pacific bluefin tuna (PBF) spawners and the amount of recruitment (stock–recruitment relationship, SRR) is unclear. It is likely that environmental effects have masked the SRR of PBF. As the basis of constructing an effective SRR for PBF, we examined the effect of spawning biomass at different ages and the spatiotemporal patterns of environmental effects on the amount of recruitment, using a recently developed model-free nonlinear time series analysis method (empirical dynamic modelling, EDM). EDM revealed where, when, and how the environment affected the amount of recruitment. EDM also found a significant contribution of ages 8–9 spawners on recruitment dynamics and that the amount of recruitment plateaus with increase in ages 8–9 spawners. Based on knowledge obtained from EDM, we formulated several example SRRs that incorporated environmental effects (sea surface temperature). The newly developed SRR with information from EDM outperformed the SRR without this information. Finally, we interpreted the results based on preceding observational and experimental studies and discussed the potential of applying the combination of EDM and mathematical modelling towards the sustainable use of other stocks. Introduction Many pelagic fish stocks have major fluctuations in their recruitment (Checkley et al., 2009; Thorson et al., 2014). Fluctuating and unpredictable recruitment leads to uncertainty in stock assessments and management. Therefore, improving the understanding of recruitment dynamics is important for robust fisheries management. The stock–recruitment relationship (SRR) is the amount of recruitment as a function of spawning stock biomass (SSB) based on the density-dependent effect of recruitment, and it is a useful tool to shrink the uncertainty in recruitment. The amount of recruitment is proportional to SSB when there is no density-dependent effect, thus the SRR is a linear function (Quinn and Deriso, 1999). With density-dependent effects, the SRR is represented as a concave function. Determining when density-dependent effects are present is difficult in some cases, and thus formulating an SRR can be intractable. Pacific bluefin tuna (Thunnus orientalis; PBF) is a typical stock with an unclear relationship between SSB and recruitment. PBF is a large pelagic fish species that migrates across the North Pacific (Ishida et al., 2018). PBF has historically been exploited and is still commercially important. Because the biomass of PBF has decreased and is at very low level in recent years (ISC, 2018), sound stock assessment is essential to manage and rebuild the stock. PBF exhibits major fluctuations in its recruitment(Figure 1), as so many other pelagic fish species. Although Nakatsuka et al. (2017) fitted the Beverton–Holt (BH) and hockey-stick SRRs to PBF recruitment and SSB for 1980–2014 (in which the uncertainty on SSB was low relative to before 1980), the SRR for PBF is still not clear as a whole. Almost no deterministic relationship between recruitment and SSB has been assumed in the stock assessments of PBF: the current SRR gives only the average and deviation of recruitment regardless of the value of SSB (ISC, 2018). One possible reason for the ambiguous SRR is that the density-dependent effect is so intense that SSB has never become low enough to cause recruitment decrease. Another potential reason is that factors other than SSB contribute to recruitment and mask the relationship between SSB and recruitment. The effect of environment is a likely candidate that masks the SRR of PBF. Environmental change is likely to affect the recruitment of PBF (Harford et al., 2017; Ishida et al., 2018; Muhling et al., 2018), but there is no available information for the strength of the effect. Therefore, whether and how the density-dependent effect and environmental effects contribute to the dynamics of recruitment is essential information to understand the mechanisms of the fluctuations in recruitment and to establish an effective SRR for PBF. Figure 1. View largeDownload slide The time series of PBF recruitment (upper panel) and SSB at age (lower panels) from 1952 to 2016. Numbers at the top left of the lower panels indicate age. Because the maturities at ages 0, 1, and 2 are zero, only SSB at ages 3–20 are depicted. Figure 1. View largeDownload slide The time series of PBF recruitment (upper panel) and SSB at age (lower panels) from 1952 to 2016. Numbers at the top left of the lower panels indicate age. Because the maturities at ages 0, 1, and 2 are zero, only SSB at ages 3–20 are depicted. Investigating the correlation between potential driving factors and recruitment might not clarify the contribution of SSB and environment on the recruitment of PBF if there are nonlinear relationships between them. In a nonlinear system, cause and effect do not always correlate: the cause may upregulate, downregulate, or not contribute to the dynamics of an effect, depending on the status of the system and thus time (Deyle et al., 2016). Therefore, factors that affects recruitment may be missed if we focus only on correlation. Empirical dynamic modelling (EDM; Ye et al., 2015; Chang et al., 2017) is a series of model-free time series analysis methods that can identify whether a certain factor contributes to the dynamics of another factor (Sugihara et al., 2012) and estimates how the effect reacts against the changes in the identified causes (Deyle et al., 2016). These features are useful to clarify the structure of the target systems from a given time series data, particularly when the nature of the relationship is unknown. Our aim in this study is to clarify the effects of spawners and the environment on PBF recruitment as a basis for constructing an SRR. Using EDM, we first identified the spawner age classes that affected recruitment. Second, using sea surface temperature (SST) as an environmental index, we identified the areas and months in which the environment affected recruitment. These processes were essential to evaluate the potential density-dependent effects and environmental effects appropriately, because including non-contributing factors would obscure the relationship between the true cause and recruitment. Then, we tracked which direction and how strong these factors affected recruitment. Based on the results from the EDM, we finally formulated SRRs that incorporated the identified density-dependent effects and environmental effects. The newly developed SRRs outperformed the one currently used in stock assessment and captured the effect of an extreme decline in SSB. Based on these results, we discussed the plausibility of the detected density-dependent effects and environmental effects, and how future stock assessment and management of pelagic fish stocks should be conducted. Methods Basic principles of EDM EDM is a series of model-free time series analysis methods for nonlinear dynamic systems. EDM predicts the near-future of the target system (Sugihara and May, 1990; Sugihara, 1994), detects whether one factor affects another or not (Sugihara et al., 2012), and measures the strength of the effect (Deyle et al., 2013,, 2016). Assume a target nonlinear dynamic system consisting of n components, the target variable y = [y(1), y(2), …] (here, it is recruitment of PBF) and n−1 forcing factors xi = [xi(1), xi(2), …] (1 ≤ i ≤ n−1) that affect the target variable directly or indirectly (they might be SSB, water temperature, etc.). No stochastic processes are assumed. Taking each factor as a coordinate axis, the dynamics of the system are represented by a trajectory in an n-dimensional state space. In the trajectory, nearby points show behaviours similar to the system. This basic principle guarantees the model-free feature of EDM. For example, we are able to predict the state of the system in time t + 1 by observing where the neighbours of [x1(t), x2(t), …, xn−1(t), y(t)] go in the next time step. However, it is generally impossible to observe all of the forcing factors. In such a situation, state space reconstruction enables us to analyse the target system without observing all of the forcing factors in the system. Consider the time delay coordinates of y, Y(t) = {y[t–(E–1)τ], y[t–(E–2)τ], …, y(t)}, where E and τ denote the embedding dimension (E ≤ t) and time delay, respectively. In this study, we fixed τ at 1, according to a preceding study (Harford et al., 2017). Takens’ theorem (Takens, 1981) guarantees that the E-dimensional trajectory of Y holds the entire information of the original trajectory if E is selected appropriately. The appropriate E is found by examining the one time-step predictions of Y (Sugihara and May, 1990; Sugihara, 1994; Sugihara et al., 2012; see Supplementary Appendix A): conducting the prediction under various values of E, the best prediction will be given under the most appropriate E. Also, a part of the elements of Y is substitutable with arbitrary xi (Deyle et al., 2013). For example, the trajectory of Y*(t) = {x1(t), x2(t), y[t–(E–3)], y[t–(E–4)], …, y(t)} preserves the entire information contained in the original trajectory. State space reconstruction analyses the trajectory of Y or Y* instead of the original trajectory, therefore enabling us to analyse the target system based on only the observation of y (and some of xi), without unobservable factors nor the knowledge of how many factors drive y. In this study, we used mainly two techniques in EDM, convergent cross mapping (CCM; Sugihara et al., 2012) and multivariate S-map (MS; Deyle et al., 2016; Ushio et al., 2018). CCM is a method specialized for detecting whether a certain factor affects another. MS finds out how the driving factors affect the target variable. The basic principles and algorithms of the two methods are described in the following sections and Supplementary Appendices A and B. Overview of CCM CCM is applied to a pair of time series, a candidate cause and a candidate effect. The candidate effect must have a nonlinear structure and not be masked by stochastic processes. The methods used to check the nonlinearity of a time series are described in Supplementary Appendix A. CCM identifies whether a candidate cause time series directly or indirectly affects a candidate effect time series. Here, we explain the intuitive mechanism of CCM (see Supplementary Appendix A for a detailed algorithm) using the target system in the previous section. x1 affects the behaviour of y and is not affected by y. Therefore, x1 is the cause and y is the effect. Suppose that we do not know the true causal relationship and suspect x1 and y as the cause and effect, respectively. Since the candidate effect time series must be stationary (Chang et al., 2017), it should be first-differenced or logarithmic-differenced before analysis, as necessary. In this case, CCM detects the causality between the candidate cause and the increase or decrease in the candidate effect. Nevertheless, an increase or decrease in the candidate effect would inevitably affect its own dynamics. The candidate cause time series x1 can be reconstructed based on information from the candidate effect time series y (cross mapping), because the trajectory of Y contains the entire information of the system. When we estimate the value of x1(t), for example, we must first find the E + 1 nearest neighbours of Y(t). The temporally corresponding x1 of the nearest neighbours will be close to each other, because the similar behaviour of y is obtained from the similar behaviour of x1. The estimate of x1(t) is calculated as an average of the temporally corresponding x1 of the nearest neighbours. This method does not work if there is no effect of x1 on y, or if the candidate cause and effect are reversed (i.e. x1 is the candidate effect and y is the candidate cause). In the latter case, y cannot be reconstructed from the information of x1 because x1 does not contain any information about y. Therefore, by examining whether the candidate cause time series is reconstructed from the information of the candidate effect time series, we can judge whether the candidate cause actually affects the candidate effect. In this study, the accuracy of reconstruction was evaluated using the correlation coefficient of the original and reconstructed candidate cause time series. A problem occurs when x1 and y highly correlate, in spite of the absence of any causal relationships (e.g. spurious correlation). In this case, x1 might be reconstructed with high accuracy if x1 affects y. The difference between this type of correlation and a true (nonlinear) causal relationship appears when the number of available Y(L) is small. When x1 truly affects y, the accuracy of the reconstruction of x1 is low with a small L, and the accuracy improves as L increases (this behaviour is called “convergence”). In contrast, if x1 and y are strongly correlated without any causality, the accuracy of the reconstruction will be high, even with a small L. Moreover, this accuracy does not improve, even when L is large. Therefore, the existence of convergence is a necessary condition for detecting causality. To confirm convergence, we compared the accuracies of the reconstruction under small and large L values. We were able to find convergence when the accuracy of the reconstruction under a large L was significantly higher than that under a small L. In addition to confirming convergence, the accuracy of the reconstruction under a large L must be significantly high to judge whether x1 actually affects y. The significance of the accuracy of cross mapping is tested by the twin surrogate method (Thiel et al., 2006, see Supplementary Appendix C), which creates many “surrogate” time series by shuffling the original time series according to a certain rule. The surrogate time series have approximately the same structure (autoregression function and mutual information with the other time series) as the original one, but any causal relationships with the other time series are removed. Using this method, we created 100 surrogate time series for every candidate cause time series. Then, we reconstructed the surrogate time series using the candidate effect time series, in the same way that we reconstructed x1. The accuracy of the reconstruction is regarded as “significantly high” when it is greater than the 95th percentile of the correlation coefficients between the reconstructed and original surrogate data. We judged that the candidate cause actually influenced the candidate effect when there was convergence and the accuracy of the reconstruction under a large L was significantly high. Note for interpretation of CCM results When CCM detects the effect of x1 on y, it cannot indicate whether the effect of x1 is direct or indirect. For example, if we detect that SST affects the dynamics of a certain fish population, SST might not be the closest driver. There might be additional factors (e.g. food availability) between SST and the fish population dynamics (Sugihara et al., 2012). In addition, if a third factor, z, is very strongly coupled with x1 (and not coupled with y) and is indistinguishable from x1 (i.e. the dynamics of z are only determined by x1), CCM might detect an effect of z on y, even though this effect does not actually exist (Cobey and Baskerville, 2016). Therefore, in a complex system (such as an ecosystem), a detected driver might be a proxy for a group of other variables (Sugihara et al., 2012). Also, CCM does not tell about the strength of the detected effects. Because the significance of the surrogate data method is determined not only by the strength of the effect but also by the strength of stochastic processes (process and observation errors), the results from CCM indicate only the presence or absence of the effect (the strength of the effect is measured by MS). Overview of MS MS is applied to a set of time series, the target variable (effect) and the causes identified by CCM. MS evaluates the strength of the effect from each cause on the target variable. Because the strength of an effect changes overtime in a nonlinear system, MS tracks it for every time point. Again, consider the dynamic system described in the previous sections. Suppose that x1 and x2 were identified as driving factors of y. The strength of an effect from xi (i = 1 or 2) to y is defined as an increase or decrease in y against a unit change of xi. Therefore, the effects from x1 and x2 to y in time t are represented by the gradient of the trajectory of Y* in the direction of the x1 and x2 axes, respectively. MS obtains this gradient by locally linearizing the original trajectory around Y*(t) and calculating Jacobian elements, Jx1(t) = ∂Y*(t)/∂x1(t) and Jx2(t) = ∂Y*(t)/∂x2(t), of the locally linearized system. The details of the local linearization and Jacobian elements calculation are described in Supplementary Appendix B. Recruitment and abundance-at-age time series of PBF PBF have been subject to stock assessment by the International science committee for tuna and tuna-like species in the North Pacific Ocean (ISC). We used a time series of recruitment and abundance-at-age (Figure 1) obtained from a stock assessment conducted in 2018 using Stock synthesis (SS, Methot and Wetzel, 2013), an integrated stock assessment model. Integrated stock assessment models estimate the most plausible status of the target stock based on available information on catch-per-unit-effort (CPUE), size and weight of the exploited fish, and so on. SSB is calculated as the sum of the multiple of the abundance and maturity for each age (0, 20, 50, and 100% for ages 0–2, 3, 4, and over 5, respectively; ISC, 2018). The abundance of age 0 fish corresponds to recruitment. If a certain SRR is assumed for a given recruitment, there is inevitably an effect from the SSB on the estimated recruitment. In this situation, usually we cannot apply CCM to clarify the true structure of the SSB and recruitment using the results obtained by the stock assessment model. However, although the BH SRR (Beverton and Holt, 1957) was used in the PBF stock assessment, we are able to regard that there was no assumption on the SRR because the steepness parameter, which determines the degree of recruitment drop against the SSB decrease, was fixed so as to keep the recruitment level constant regardless of the SSB level (steepness = 0.999). Under this setting, BH gives almost only the average and standard deviation of recruitment (Nakatsuka et al., 2017; ISC, 2018) and the amount of recruitment is determined almost only from the input data. In the stock assessment of PBF, PBF is assumed to spawn on 1 April and to recruit on 1 July. Because SS is based on a fishing year that starts on 1 July, the abundance of age 0 fish in a certain fishing year can be affected by the SSB in the previous fishing year. Taking this into account, we used recruitment from 1953 to 2015 [R(1–63)] and SSB at each age (provided by ISC) from 1953 to 2014 [Sa(1–62), 3 ≤ a ≤ 20; note that ages 0–2 fish are immature]. Although the period subject to the latest stock assessment is 1952–2016, we excluded data of the first and last years according to a preceding study (Harford et al., 2017). Target region and SST Two PBF spawning grounds were identified (Figure 2). One is in the southwestern North Pacific, spreading from off east Taiwan to the Nansei Islands (Uematsu et al., 2018; Nansei spawning ground in this study). In this spawning ground, mainly age 9 or older fish spawn from April to June (Ashida et al., 2015). The other spawning ground is in the Sea of Japan, where mainly ages 3–5 fish spawn in July and August (Okochi et al., 2016). Juveniles are carried northward by the Kuroshio and Tsushima current before recruitment (Kitagawa et al., 2010; Masujima et al., 2014), then stay in the coastal area of southern Japan, the East China Sea, and the Sea of Japan until recruitment (Rooker et al., 2001; Figure 2). Some of the recruits are exploited by the troll fishery in the East China Sea (Figure 2), mainly from autumn to spring (Ichinokawa et al., 2014). Their CPUE data are used to inform the amount of recruitment in the stock assessment. We downloaded the SST data for the area from 121°–150°E to 20°–45°N in the North Pacific (http://ds.data.jma.go.jp/tcc/tcc/products/elnino/cobesst/cobe-sst.html, Ishii et al., 2005), which covers the two spawning grounds and the migration routes of the juveniles. We created monthly average SST time series from April to December (1954–2015) and January to March (1955–2016) for each one-degree block, obtaining 626 SST time series per month. Figure 2. View largeDownload slide The spawning areas and migration routes of PBF. The shadowed areas indicate the Nansei and Sea of Japan spawning grounds. The dotted area indicates the East China Sea age 0 troll fishing ground. Thin black and thick grey arrows show the migration route of larvae before recruitment and the path of Kuroshio, respectively. Depicted based on Nakatsuka et al. (submitted), Uematsu et al. (2018), and Muhling et al. (2018). Figure 2. View largeDownload slide The spawning areas and migration routes of PBF. The shadowed areas indicate the Nansei and Sea of Japan spawning grounds. The dotted area indicates the East China Sea age 0 troll fishing ground. Thin black and thick grey arrows show the migration route of larvae before recruitment and the path of Kuroshio, respectively. Depicted based on Nakatsuka et al. (submitted), Uematsu et al. (2018), and Muhling et al. (2018). Applying EDM to the recruitment, SSB at age, and SST data We checked the nonlinearity and estimated the appropriate embedding dimension (E) for R and its first-difference, Rdiff(1–62) (Supplementary Appendices A and B). Both time series showed nonlinearity. The estimated embedding dimensions for R and Rdiff were 6 and 2, respectively. We conducted a series of CCM analyses to identify the factors affecting the amount of recruitment. According to preceding studies (Sugihara et al., 2012; Nakayama et al., 2018; Ushio et al., 2018), we set Rdiff as the candidate effect in the CCM analysis. First, to find the month and region in which SST drives recruitment, we set the SST of every one-degree block and month as the candidate cause. Then, to clarify whether and which age of spawning biomass affected recruitment, we set Sa (3 ≤ a ≤ 20) as the candidate cause. Note that the maturity rate for each age used in the stock assessment (ISC, 2018) does not affect the results of the CCM, because maturity rate determines only the scale, not the pattern, of the dynamics. All candidate cause and effect time series were normalized before these analyses, so as to have means of zero and standard deviations of 1. We found significant causality when the correlation coefficient between the original and reconstructed candidate cause time series under L = 60 was significantly higher than that under L = 30 (i.e. convergence) and the 95th percentile of the confidence interval for the correlation coefficients between the original and reconstructed surrogate data. Then, we conducted MS to clarify how the identified factors affected the amount of recruitment. We reconstructed the six-dimensional state space using the time series of recruitment, SSB, and SST. The original (not first-differenced) recruitment time series [R(1–63)] was embedded in the six-dimensional state space and then the 5th and 6th dimensions were substituted by the time series of SST and SSB. Because the CCM analysis identified 15 combinations of month and region where the effect of SST to recruitment was significant (see Results section and Figure 3), the 5th dimension consisted of the spatially averaged SST from the 15 combinations of the regions and months, Ti(1–63) (1 ≤ i ≤ 15). From the CCM analysis we also saw a significant effect of the ages 8–9 SSB on recruitment (see Results section and Figure 4). We set the sum of the ages 8–9 SSB [S8–9(1–62)] as the last dimension. Therefore, we reconstructed 15 state spaces [R(6–63), R(5–62), R(4–61), R(3–60), Ti(6–63), S8–9(5–62)]. The effect of Ti and S8–9, JT = ∂R/∂Ti and JS = ∂R/∂S8–9, respectively, were approximately calculated by MS for each state space. Also, because the partial derivatives JT and JS are embedding-dependent, we confirmed their robustness by applying MS to state spaces that excluded the information of SSB [R(6–63), R(5–62), R(4–61), R(3–60), R(2–59), Ti(6–63)] and SST [R(6–63), R(5–62), R(4–61), R(3–60), R(2–59), S8–9(5–62)], respectively. In MS analyses, we did not normalize the concerned time series because we were interested in the concrete values of the SST and SSB that give certain effects on recruitment. JS might contain the density-dependent effect derived from the underlying BH assumed in the current stock assessment, although a steepness of 0.999, under which recruitment is hardly affected by the SRR, was assumed. To estimate the strength of this potential effect (JBH) and compare it with JS, we calculated the theoretical gradient of the BH against S8–9. JBH was obtained as the partial derivative of BH [parameters were set as the same as the latest stock assessment (ISC, 2018)] with respect to SSB, expanded by the ratio of total SSB (Stotal) and S8–9. JT and JS show how R has changed against the unit change of Ti and S8–9, respectively, but the relative strengths of these effects are also of interest. To compare the strength of the effect from Ti and S8–9, we defined and compared standardized effect strengths, JT* = JTSD(Ti) and Js* = JsSD(S8–9), where SD(Ti) and SD(S8–9) denote the standard deviation of Ti and S8–9, respectively. Figure 3. View largeDownload slide Results of the CCM applied to monthly averaged SST and PBF recruitment. The one degree mesh is coloured where significant (p < 0.05) effects of SST were detected. The colours indicate the significance of the detected effects, based on the colour bar at the bottom. Blue rectangles and the numbers inside are the four regions referred to in the text. Figure 3. View largeDownload slide Results of the CCM applied to monthly averaged SST and PBF recruitment. The one degree mesh is coloured where significant (p < 0.05) effects of SST were detected. The colours indicate the significance of the detected effects, based on the colour bar at the bottom. Blue rectangles and the numbers inside are the four regions referred to in the text. Figure 4. View largeDownload slide Results of the MS applied to the SST in Region 1, July (T3, shown by the dashed line) and the recruitment of PBF. The strength of the effect (JT, shown by diamonds) is plotted against year (a) and SST (b). Figure 4. View largeDownload slide Results of the MS applied to the SST in Region 1, July (T3, shown by the dashed line) and the recruitment of PBF. The strength of the effect (JT, shown by diamonds) is plotted against year (a) and SST (b). Formulating example SRRs based on the results from EDM We formulated example SRRs based on the observed density-dependent effect and environmental effect. MS analysis indicated a plateau against the increase in S8–9. Also, the SST from 130°–140°E to 26°–34°N in July (T3) showed the largest, constant positive effect on recruitment (JT*) in the MS analysis (see Results section). Therefore, we incorporated this information into the traditional BH SRR. The original BH is formulated as Rt=4hR0S(t-1)1-hSSB0+5h-1S(t-1), where S(t−1), R0, SSB0, and h are SSB in year t−1, recruitment in the absence of exploitation, total SSB in the absence of exploitation, and steepness, respectively. R0 and SSB0 are dependent on a constant, spawner per recruitment in the absence of exploitation (SPR0) SSB0=SPR0R0. We incorporated the effect of SST into the error term: Rt=4hR0S(t-1)1-hSSB0+5h-1S(t-1)eε, ε∼NaTi*(t), σ2, where a, T*i(t), and σ denote the effect of SST, the standardized original SST time series in year t (i.e. standardized Ti, 1 ≤ i ≤ 15), and standard deviation, respectively. We used T*3, the standardized T3, which had the largest JT* in the MS analysis (see Results section). In this model, S(t−1) denotes SSB at arbitrary ages in year t−1 and SSB0 is defined as the SSB at corresponding ages under absence of exploitation. To identify which modification improves the prediction performance, we compared the Akaike information criterion (AIC) of the following models. Model 1-1: S was set as S8–9 and a was estimated. Model 1-2: Identical to Model 1-1, but S was set as Smin–max, where min and max were searched to obtain the least AIC. Model 2: S was set as Stotal and a was estimated. Model 3-1: S was set as S8–9 and a was fixed at 0. Model 3-2: Identical to Model 3-1, but S was set as Smin–max, where min and max were searched to obtain the smallest AIC. Model 4: S was set as Stotal and a was fixed at 0 (identical to the ordinary BH). For Models 1-2 and 3-2, we searched all possible sequential combinations of ages that constituted SSB. This process was essential to develop a more predictive SRR, because CCM does not always propose all of the factors that contribute to the dynamics of recruitment. The parameters were estimated using a maximum likelihood method. As in the EDM analysis, we constructed these models with recruitment and SST until 2015 and SSB until 2014. In addition, we also estimated the amount of recruitment in 2016 and compared it with that obtained through the stock assessment, using these models with the SST in the corresponding region and month in 2016 and SSB at the corresponding age in 2015. Results Effect of SST and SSB We identified 15 major combinations of region and month in which SST affected the amount of recruitment (Figure 3). The effect of SST was detected mainly in four regions: 128°–139°E 29°–34°N (in May, June, July, and August; Region 1, corresponding spatially averaged SST time series were denoted by T1–4, respectively), 120°–130°E 20°–25°N (June; Region 2, corresponding spatially averaged SST time series was denoted by T5), 140°–150°E 32°–45°N (June, July, October, November, and December; Region 3, corresponding spatially averaged SST time series were denoted by T6–10, respectively), and 130°–140°E 32°–45°N (October, November, December, January, and February; Region 4, corresponding spatially averaged SST time series were denoted by T11–15, respectively). For these 15 combinations of region and month, the spatially averaged correlation coefficients between the original and reconstructed recruitment under L = 60 varied from 0.18 to 0.39 (Supplementary Table D1). MS revealed almost constant positive effects of SST on the amount of recruitment in Regions 1, 3, and 4, whereas a negative effect was observed in Region 2 (Figure 4 and Supplementary Figure E1). A significant contribution of SSB to recruitment dynamics at ages 8 and 9 was detected by the CCM analysis (Figure 5). The sum of SSB at ages 8 and 9 [S8–9(t)] had positive effects almost throughout the stock assessment period (Figure 6a). The effect of S8–9 diminished and became nearly zero with increasing S8–9 (Figure 6b), indicating a density-dependent effect under which the recruitment reaches a certain plateau when S8–9 becomes large. These tendencies were robust, regardless of the region and month from where the time series of SST were sampled (Supplementary Figure E2). The detected effect of SSB at ages 8 and 9 was obviously larger than the potential effect from the assumed BH, when S8–9 was <20 000 t (Figure 6b andSupplementary Figure E2). Even when we excluded the information of SSB and SST from the state space, we obtained similar results on the effect strength of SST and SSB, respectively (Supplementary Figures E3 and E4). The standardized effect strength of SST (JT*) was the largest in Region 1 in July. In this region and month, JT* was always larger than Js* (Figure 7), indicating a larger contribution of SST than the amount of age 8–9 spawners on recruitment dynamics. Figure 5. View largeDownload slide Results of the CCM applied to the SSB at each age (ages 3–20, indicated by the number at the top left of each cell) and PBF recruitment. The horizontal axis represents the number of embedded vectors used for cross mapping (L). The solid line represents the mean correlation coefficient between the original and reconstructed SSB at age time series (Sa, 3 ≤ a ≤ 20). The dashed line and shadow in each cell represents the mean and 95th percentile of the results of the CCM for the surrogate data vs. the original recruitment time series (Rdiff), respectively. Figure 5. View largeDownload slide Results of the CCM applied to the SSB at each age (ages 3–20, indicated by the number at the top left of each cell) and PBF recruitment. The horizontal axis represents the number of embedded vectors used for cross mapping (L). The solid line represents the mean correlation coefficient between the original and reconstructed SSB at age time series (Sa, 3 ≤ a ≤ 20). The dashed line and shadow in each cell represents the mean and 95th percentile of the results of the CCM for the surrogate data vs. the original recruitment time series (Rdiff), respectively. Figure 6. View largeDownload slide Results of the MS applied to the SSB at ages 8–9 (S8–9) and PBF recruitment. The strength of the effect (JS, solid circles) is plotted against year (a, S8–9, is also shown by the dashed line) and SSB at ages 8–9 (b). The theoretical effect derived from the BH (JBH) used in the stock assessment is shown as crosses. Figure 6. View largeDownload slide Results of the MS applied to the SSB at ages 8–9 (S8–9) and PBF recruitment. The strength of the effect (JS, solid circles) is plotted against year (a, S8–9, is also shown by the dashed line) and SSB at ages 8–9 (b). The theoretical effect derived from the BH (JBH) used in the stock assessment is shown as crosses. Figure 7. View largeDownload slide The relative strength of the effects from the SST in Region 1 in July (T3, diamonds) and the SSB at ages 8–9 (S8–9, solid circles) on recruitment. Figure 7. View largeDownload slide The relative strength of the effects from the SST in Region 1 in July (T3, diamonds) and the SSB at ages 8–9 (S8–9, solid circles) on recruitment. Performances of the formulated SRRs In all cases, the parameters were successfully estimated (the maximum absolute value of the gradients of the log-likelihood against each parameter <0.0001; the estimated values are shown in Supplementary Table F1). For Models 1-2 and 3-2, the range of ages (min and max) that obtained the smallest AIC was ages 8–15 for both models. The models with S8–9(t) (Models 1 and 3) led to small increases in AIC compared to the models with total SSB (1.1 and 0.8 increases in AIC relative to Models 2 and 4, respectively). However, the use of S8–15(t) (Models 1-2 and 3-2) improved the performances of the models (7.51 and 4.5 decreases in AIC relative to Models 2 and 4, respectively). In all models in which the effect of SST, a, was not fixed at 0 (Models 1, 1-2, and 2), a was estimated as a positive value, indicating a positive effect of SST (Supplementary Table F1 and Figure 8a–c). This result is consistent with the results from MS analysis. Inclusion of SST in the SRRs (Models 1, 1-2, and 3) led to a drastic drop of AIC (13.7, 17.0, and 14.0 decreases in AIC relative to Models 3, 3-2, and 4, respectively). All models other than Model 4 more or less captured the true increase in recruitment in 2016, relative to 2015 (Figures 1 and 7g–l). This occurred because of two reasons: SSB at ages 8–9 and 8–15 in 2016 increased by 47.6 and 43.5% relative to those in 2015, respectively, whereas total SSB increased only 3.5%; and the corresponding SST (Region 1, July) in 2016 was high relative to that in 2015. Discussion In this study, we identified the regions and months where SST affected the amount of PBF recruitment. We iterated hundreds of statistical tests using surrogate data, therefore some results might contain false positives associated with the multiple comparisons. However, the effects of SST were detected in patches and the effects of SSB were detected in sequential ages (Figures 3 and 5), thus the general patterns of the results are likely robust. Region 2 in June, where a negative effect of SST was observed (Figure 3 and Supplementary Figure A2), spatially and temporally overlapped with the Nansei spawning ground (Ashida et al., 2015). A rearing experiment (Miyashita et al., 2000) revealed that the hatching rate and normal hatching rate of PBF eggs have an optimum temperature around 26°C. Since the SST in this region in June was higher than 26°C (mean 28.2°C, standard deviation 0.4°C), the negative effect of SST obtained from MS (Supplementary Figure A1) was consistent with the physiology of hatching. In addition, there was no clear linear correlation between summer SST and recruitment in this region (Muhling et al., 2018). Therefore, focusing on the nonlinear behaviour of the ecosystem might be key to understanding the effect of environment on egg hatching and early life stage survival in PBF. Nevertheless, comparing the absolute values of JT, the effect of SST was larger in Region 1 than it was in Region 2 (Figure 4 and Supplementary Figure A2; although we showed the non-standardized values, JT* also showed the same tendency). Region 1 overlapped with the nursery area of small juveniles (Rooker et al., 2001; Ohshimo et al., 2017), and a positive correlation between recruitment and summer SST was also observed (Muhling et al., 2018). However, an archival tag tracking (Fujioka et al., 2018) found no small juveniles moved from Region 1 to the East China Sea, whose fishery CPUE partly informs the recruitment estimate in an integrated assessment model. Therefore, SST in Region 1 cannot affect the recruitment time series data we used in this study, although it might contribute to the dynamics of actual recruitment. It is likely that SST in Region 1 is a proxy of other environmental factors that affect recruitment dynamics. One plausible candidate of these factors is the path of Kuroshio. A numerical particle tracking study (Kitagawa et al., 2010) indicated that the Kuroshio path affects the transportation efficiency of PBF larvae and small juveniles to their nursery area in the East China Sea. Also, Kuroshio variability affects the survival of PBF larvae during transportation (Kimura et al., 2010). Therefore, Kuroshio is likely to directly affect the early survival and recruitment of PBF. Kuroshio then transports warm water to Region 1 (Figure 2), thus it is possible that SST in Region 1 works as a proxy index of the efficiency of larvae and juvenile transportation to the East China Sea. Validating this hypothesis is important, because SST in Region 1 largely improved the performance of the SRRs we formulated (Figure 8 and Supplementary Table F1) and a statistical model for recruitment prediction (Muhling et al., 2018). There are no physiological explanations for the effects of SST detected in Regions 3 and 4. The SSTs in Regions 3 and 4 are strongly correlated, and the effects of SST were the strongest in Region 4 in January (Supplementary Figure E2). Therefore, the driver closest to recruitment is likely to be the winter SST in Region 4. Because the age 0 PBF fishery in the East China Sea is operated during winter (Ichinokawa et al., 2014), it is possible that the SST in Region 4 affects the recruitment time series from the stock assessment, through the age 0 fish distribution and the fishery CPUE, rather than the actual amount of recruitment. Figure 8. View largeDownload slide The SRRs (a–f) and predicted historical recruitment (g–l) of Model 1 (a and g), Model 1-2 (b and h), Model 2 (c and i), Model 3 (d and j), Model 3-2 (e and k), and Model 4 (f and l). Because the SRRs of Models 1, 1-2, and 2 (a–c) included information on SST, the SSR corresponding to the 95th percentile (red solid line), mean (black solid line), and 5th percentile (blue solid line) of SST are shown with the 95% prediction intervals of the corresponding colours. In other cells, the SRRs (d–f) and predicted recruitment (g–l) are show by the black solid lines. Grey shadows indicate 95% prediction intervals of the SRRs (d–f) and predicted recruitment (g–l). The black solid circles indicate the annual recruitment estimated in the stock assessment. Figure 8. View largeDownload slide The SRRs (a–f) and predicted historical recruitment (g–l) of Model 1 (a and g), Model 1-2 (b and h), Model 2 (c and i), Model 3 (d and j), Model 3-2 (e and k), and Model 4 (f and l). Because the SRRs of Models 1, 1-2, and 2 (a–c) included information on SST, the SSR corresponding to the 95th percentile (red solid line), mean (black solid line), and 5th percentile (blue solid line) of SST are shown with the 95% prediction intervals of the corresponding colours. In other cells, the SRRs (d–f) and predicted recruitment (g–l) are show by the black solid lines. Grey shadows indicate 95% prediction intervals of the SRRs (d–f) and predicted recruitment (g–l). The black solid circles indicate the annual recruitment estimated in the stock assessment. We also detected an effect of SSB at ages 8–9 on recruitment (Figure 5). Because SSB at age is calculated as the sum of the matured spawner weight in April (ISC, 2018) and PBF are born from April to June (Nansei spawning ground) or from July to August (the Sea of Japan spawning ground), spawners near ages 9 and 10 were classified into ages 8 and 9 in this study. Ashida et al. (2015) reported that most of the spawners in the Nansei spawning ground were over 200 cm, which corresponds to an age over 9 according to the known von Bertalanffy growth curve for PBF (Shimose et al., 2009). Therefore, the age classes affecting the amount of recruitment corresponded to the youngest age classes in the Nansei spawning ground. In contrast, Okochi et al. (2016) reported that smaller (thus, younger) fish mainly spawn in the Sea of Japan spawning ground. The present result highlights the relative importance of the Nansei spawning ground for PBF recruitment. The effect of SST in Region 1 in July was larger than the effect of SSB throughout the period of the stock assessment (Figure 7). Nevertheless, the magnitude of the effect of age 8–9 spawners reached approximately half of the effect of SST, therefore the abundance of these age classes should be regarded as an important factor that influences the amount of recruitment. The combination of EDM and mathematical modelling was effective to find the “useful” explanatory variables and describe the complicated relationships among them. The performances of the SRRs developed in this study, which were built based on insight from EDM, were improvements on the standard BH (Model 4). Although the use of SSB at ages 8–9 detected by CCM did not improve the model performance, expanding the age classes of SSB to 8–15 obtained the best model. Spawners over age 9 have been observed at the Nansei spawning ground (Ashida et al., 2015). One likely reason for needing an age class expansion is that CCM might miss the contribution of spawners over age 10. When multiple age classes contribute to recruitment, they affect in the aggregate, not independently. In this situation, the reconstruction of single age class dynamics by cross mapping is intractable even if the age class actually contributed to the recruitment, and thus, the detection power of CCM might be lower than usual. It is likely that SSB at ages 8–9 were detected because the age classes well represented the dynamics of the total population contributing to recruitment. Because SSB gradually declines with increasing age (Figure 1), the younger age classes are more likely to represent the entire group of spawners if multiple ages of spawners contribute to recruitment. Therefore, different results might be obtained if we set all possible subsets of the age classes as the candidate cause. Nevertheless, the effects detected in this study are still guaranteed because the false positive rate of CCM is very low even with intense observation and process errors (Nakayama et al., 2015). Another possible reason for needing an age class expansion is that the BH in the models we formulated might not precisely express the true relationship between recruitment and SSB. We adopted BH simply because the obtained effect of SSB at ages 8–9 indicated a plateau in recruitment, but there might be a more appropriate formulation of SRR. Observing the effect strength of SSB with SST fixed at certain levels will give us more insight because the currently observed effect strengths of SSB (Figure 6) were calculated using various SST values. One of the EDM methods, “scenario exploration” (Deyle et al., 2013), is applicable for this purpose. It predicts the dynamics of recruitment under arbitrary SSTs, and thus, calculates the effect strength of SSB under constant levels of SST. Although we found a complicated spatiotemporal pattern effect of SST using EDM (Figures 3 and 4 andSupplementary Figure E2), we included it into the SRR in a simple way, adopting only a single SST time series. Again, the scenario exploration technique might provide some insight by observing the effect strength of SST under fixed levels of SSB. Also, even a simple weighted average of SST in several areas and months has a potential to improve the SRR. Developing such methods that incorporate the findings from EDM to mathematical models will be the focus of our future efforts. The process of formulating SRR demonstrated in this study may be applied to other fish stocks, which have obscure SRRs due to multiple entangled factors. Because recent meta-analysis studies applying EDM (Munch et al., 2018; Pierre et al., 2018) and statistical models (Vert-pre et al., 2013) to fishery stocks found that the effects of stock size and environment to the recruitment are general, the procedure we proposed in this study has a potential to improve the performance of SRR in many stocks. At the same time, understanding the limits of this method is important for appropriate applications. For example, EDM does not tell whether the detected effects of SST are direct or indirect: SST might directly affect recruitment, or indirectly via the amount of food available for PBF larvae, for example. Also, the processes causing the density-dependent effects are still unknown (cannibalism and/or starvation are likely candidates: Tanaka et al., 2008; Uriarte et al., 2019). Such information will need to be obtained by empirical and experimental study. Therefore, a combination of our methods and empirical/experimental studies will lead to a better understanding of ecosystems and the sustainable use of fishery resources. Acknowledgements We are grateful to Masayuki Ushio for providing helpful advice on the application of the EDM. We are also grateful to Masato Abe, Aigo Takeshige, Kenji Morinaga, Tetsuya Akita, Momoko Ichinokawa, Hiroshi Okamura, Yukimasa Ishida, and other laboratory members for providing useful comments on an earlier version of the manuscript. This research was partly supported by the Promotion Program for International Resources Survey of the Fisheries Agency of Japan. References Ashida H. , Suzuki N. , Tanabe T. , Suzuki N. , Aonuma Y. 2015 . Reproductive condition, batch fecundity, and spawning fraction of large Pacific bluefin tuna Thunnus orientalis landed at Ishigaki Island, Okinawa, Japan . Environmental Biology of Fishes , 98 : 1173 – 1183 . Google Scholar Crossref Search ADS WorldCat Beverton R. J. , Holt S. J. 1957 . On the dynamics of exploited fish populations, Fishery Investigations Series II, Vol. XIX, Ministry of Agriculture . Fisheries and Food , 1 : 957 . WorldCat Chang C. W. , Ushio M. , Hsieh C. H. 2017 . Empirical dynamic modeling for beginners . Ecological Research , 32 : 785 – 796 . Google Scholar Crossref Search ADS WorldCat Checkley D. , Alheit J. , Oozeki Y. , Roy C. 2009 . Climate Change and Small Pelagic Fish . Cambridge University Press , Cambridge . Google Preview WorldCat Cobey S. , Baskerville E. B. 2016 . Limits to causal inference with state-space reconstruction for infectious disease . PLoS One , 11 : e0169050. Google Scholar Crossref Search ADS PubMed WorldCat Deyle E. R. , Fogarty M. , Hsieh C. H. , Kaufman L. , MacCall A. D. , Munch S. B. , Peretti C. T. et al. 2013 . Predicting climate effects on Pacific sardine. Proceedings of the National Academy of Sciences, 110 : 6430 – 6435 . Deyle E. R. , May R. M. , Munch S. B. , Sugihara G. 2016 . Tracking and forecasting ecosystem interactions in real time . Proceedings of the Royal Society B , 283 : 20152258. Google Scholar Crossref Search ADS PubMed WorldCat Fujioka K. , Fukuda H. , Tei Y. , Okamoto S. , Kiyofuji H. , Furukawa S. , Takagi J. et al. 2018 . Spatial and temporal variability in the trans-Pacific migration of Pacific bluefin tuna (Thunnus orientalis) revealed by archival tags . Progress in Oceanography , 162 : 52 – 65 . Google Scholar Crossref Search ADS WorldCat Harford W. J. , Karnauskas M. , Walter J. F. , Liu H. 2017 . Non‐parametric modeling reveals environmental effects on bluefin tuna recruitment in Atlantic, Pacific, and Southern Oceans . Fisheries Oceanography , 26 : 396 – 412 . Google Scholar Crossref Search ADS WorldCat Ichinokawa M. , Okamura H. , Oshima K. , Yokawa K. , Takeuchi Y. 2014 . Spatiotemporal catch distribution of age-0 Pacific bluefin tuna Thunnus orientalis caught by the Japanese troll fishery in relation to surface sea temperature and seasonal migration . Fisheries Science , 80 : 1181 – 1191 . Google Scholar Crossref Search ADS WorldCat ISC. 2018 . Stock assessment of Pacific bluefin tuna (Thunnus orientalis) in the Pacific Ocean in 2018. In Report in the 18th Meeting of the International Scientific Committee for Tuna and Tuna-like Species in the North Pacific Ocean Yeosu, Republic of Korea. 11–16 July. Ishida Y. , Fukuda H. , Fujioka K. , Sakai O. , Hiraoka Y. , Oshima K. , Nakatsuka S. et al. 2018 . Long‐term changes in recruitment of age‐0 Pacific bluefin tuna (Thunnus orientalis) and environmental conditions around Japan . Fisheries Oceanography , 27 : 41 – 48 . Google Scholar Crossref Search ADS WorldCat Ishii M. , Shouji A. , Sugimoto S. , Matsumoto T. 2005 . Objective analyses of sea-surface temperature and marine meteorological variables for the 20th century using ICOADS and the Kobe collection . International Journal of Climatology , 25 : 865 – 879 . Google Scholar Crossref Search ADS WorldCat Kimura S. , Kato Y. , Kitagawa T. , Yamaoka N. 2010 . Impacts of environmental variability and global warming scenario on Pacific bluefin tuna (Thunnus orientalis) spawning grounds and recruitment habitat . Progress in Oceanography , 86 : 39 – 44 . Google Scholar Crossref Search ADS WorldCat Kitagawa T. , Kato Y. , Miller M. J. , Sasai Y. , Sasaki H. , Kimura S. 2010 . The restricted spawning area and season of Pacific bluefin tuna facilitate use of nursery areas: a modeling approach to larval and juvenile dispersal processes . Journal of Experimental Marine Biology and Ecology , 393 : 23 – 31 . Google Scholar Crossref Search ADS WorldCat Masujima M. , Kato Y. , Segawa K. H. 2014 . Numerical studies focusing on the early life stage of Pacific bluefin tuna (Thunnus orientalis) . Bulletin of Fisheries Research Agency , 38 : 51 – 55 . WorldCat Methot R. D. Jr , Wetzel C. R. 2013 . Stock synthesis: a biological and statistical framework for fish stock assessment and fishery management . Fisheries Research , 142 : 86 – 99 . Google Scholar Crossref Search ADS WorldCat Miyashita S. , Tanaka Y. , Sawada Y. , Murata O. , Hattori N. , Takii K. , Mukai Y. et al. 2000 . Embryonic development and effects of water temperature on hatching of the bluefin tuna, Thunnus thynnus . Aquaculture Science , 48 : 199 – 207 (in Japanese). WorldCat Muhling B. A. , Tommasi D. , Ohshimo S. , Alexander M. A. , DiNardo G. , and Handling editor: Manuel Hidalgo. 2018 . Regional-scale surface temperature variability allows prediction of Pacific bluefin tuna recruitment . ICES Journal of Marine Science , 75 : 1341 – 1352 . Google Scholar Crossref Search ADS WorldCat Munch S. B. , Giron‐Nava A. , Sugihara G. 2018 . Nonlinear dynamics and noise in fisheries recruitment: a global meta‐analysis . Fish and Fisheries , 19 : 964 – 973 . Google Scholar Crossref Search ADS WorldCat Nakatsuka S. , Ishida Y. , Fukuda H. , Akita T. 2017 . A limit reference point to prevent recruitment overfishing of Pacific bluefin tuna . Marine Policy , 78 : 107 – 113 . Google Scholar Crossref Search ADS WorldCat Nakayama S. I. , Abe M. , Okamura H. 2015 . An introduction to convergent cross mapping: a novel method for causality detection in ecological time series . Japanese Journal of Ecology , 65 : 241 – 253 (in Japanese). WorldCat Nakayama S. I. , Takasuka A. , Ichinokawa M. , Okamura H. 2018 . Climate change and interspecific interactions drive species alternations between anchovy and sardine in the western North Pacific: detection of causality by convergent cross mapping . Fisheries Oceanography , 27 : 312 – 322 . Google Scholar Crossref Search ADS WorldCat Ohshimo S. , Tawa A. , Ota T. , Nishimoto S. , Ishihara T. , Watai M. , Satoh K. et al. 2017 . Horizontal distribution and habitat of Pacific bluefin tuna, Thunnus orientalis, larvae in the waters around Japan . Bulletin of Marine Science , 93 : 769 – 787 . Google Scholar Crossref Search ADS WorldCat Okochi Y. , Abe O. , Tanaka S. , Ishihara Y. , Shimizu A. 2016 . Reproductive biology of female Pacific bluefin tuna, Thunnus orientalis, in the Sea of Japan . Fisheries Research , 174 : 30 – 39 . Google Scholar Crossref Search ADS WorldCat Pierre M. , Rouyer T. , Bonhommeau S. , Fromentin J. M. , and Handling editor: Mikko Heino. 2018 . Assessing causal links in fish stock–recruitment relationships . ICES Journal of Marine Science , 75 : 903 – 911 . Google Scholar Crossref Search ADS WorldCat Quinn T. J. , Deriso R. B. 1999 . Quantitative Fish Dynamics . Oxford University Press , New York . Google Preview WorldCat Rooker J. R. , Secor D. H. , Zdanowicz V. S. , Itoh T. 2001 . Discrimination of northern bluefin tuna from nursery areas in the Pacific Ocean using otolith chemistry . Marine Ecology Progress Series , 218 : 275 – 282 . Google Scholar Crossref Search ADS WorldCat Shimose T. , Tanabe T. , Chen K. S. , Hsu C. C. 2009 . Age determination and growth of Pacific bluefin tuna, Thunnus orientalis, off Japan and Taiwan . Fisheries Research , 100 : 134 – 139 . Google Scholar Crossref Search ADS WorldCat Sugihara G. 1994 . Nonlinear forecasting for the classification of natural time series . Philosophical Transactions of the Royal Society of London A , 348 : 477 – 495 . Google Scholar Crossref Search ADS WorldCat Sugihara G. , May R. 1990 . Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series . Nature , 344 : 734 – 741 . Google Scholar Crossref Search ADS PubMed WorldCat Sugihara G. , May R. , Ye H. , Hsieh C. H. , Deyle E. , Fogarty M. , Munch S. 2012 . Detecting causality in complex ecosystems . Science , 338 : 496 – 500 . Google Scholar Crossref Search ADS PubMed WorldCat Takens F. 1981 . Detecting strange attractors in turbulence. In Dynamical Systems and Turbulence, Warwick 1980: Proceedings of a Symposium Held at the University of Warwick 1979/80, 898, pp. 366 – 381 . Ed. by D. A. Rand and L. S. Young. Springer , Berlin, Heidelberg . Tanaka Y. , Satoh K. , Yamada H. , Takebe T. , Nikaido H. , Shiozawa S. 2008 . Assessment of the nutritional status of field-caught larval Pacific bluefin tuna by RNA/DNA ratio based on a starvation experiment of hatchery-reared fish . Journal of Experimental Marine Biology and Ecology , 354 : 56 – 64 . Google Scholar Crossref Search ADS WorldCat Thorson J. T. , Jensen O. P. , Zipkin E. F. 2014 . How variable is recruitment for exploited marine fishes? A hierarchical model for testing life history theory . Canadian Journal of Fisheries and Aquatic Sciences , 71 : 973 – 983 . Google Scholar Crossref Search ADS WorldCat Thiel M. , Romano M. C. , Kurths J. , Rolfs M. , Kliegl R. 2006 . Twin surrogates to test for complex synchronisation . EPL , 75 : 535. Google Scholar Crossref Search ADS WorldCat Uematsu Y. , Ishihara T. , Hiraoka Y. , Shimose T. , Ohshimo S. 2018 . Natal origin identification of Pacific bluefin tuna (Thunnus orientalis) by vertebral first annulus . Fisheries Research , 199 : 26 – 31 . Google Scholar Crossref Search ADS WorldCat Uriarte A. , Johnstone C. , Laiz-Carrión R. , García A. , Llopiz J. K. , Shiroza A. , Quintanilla J. M. et al. 2019 . Evidence of density-dependent cannibalism in the diet of wild Atlantic bluefin tuna larvae (Thunnus thynnus) of the Balearic Sea (NW-Mediterranean) . Fisheries Research , 212 : 63 – 71 . Google Scholar Crossref Search ADS WorldCat Ushio M. , Hsieh C-H. , Masuda R. , Deyle E. R. , Ye H. , Chang C.-W. , Sugihara G. et al. 2018 . Fluctuating interaction network and time-varying stability of a natural fish community . Nature , 554 : 360. Google Scholar Crossref Search ADS PubMed WorldCat Vert-pre K. A. , Amoroso R. O. , Jensen O. P. , Hilborn R. 2013 . Frequency and intensity of productivity regime shifts in marine fish stocks . Proceedings of the National Academy of Sciences , 110 : 1779 – 1784 . Google Scholar Crossref Search ADS WorldCat Ye H. , Beamish R. J. , Glaser S. M. , Grant S. C. , Hsieh C. H. , Richards L. J. , Schnute J. T. et al. 2015 . Equation-free mechanistic ecosystem forecasting using empirical dynamic modeling . Proceedings of the National Academy of Sciences , 112 : E1569 – 1576 . Google Scholar Crossref Search ADS WorldCat © International Council for the Exploration of the Sea 2019. All rights reserved. For permissions, please email: [email protected] This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)