Contrasting life-history responses to climate variability in eastern and western North Pacific sardine populations

All procedures accorded to administrative provision of animal welfare of the Fisheries Research Education Agency Japan. All statistical tests used in this study are two-sided.

Otolith samples

From the western North Pacific, age-0 JP sardine were collected from samples taken during acoustic and sub-surface trawl surveys in the offshore Oyashio region conducted during 2006–2010 and 2014–2015. The surveys were conducted by Japan Fisheries Research and Education Agency every autumn since 2005 which aim to estimate the abundance of small pelagic species. The abundance of young-of-the-year sardine in the region in the season, approximately 10–15 cm in standard length (SL), is considered a proxy for the abundance of recruits of the Pacific stock and used to tune the cohort analysis in stock assessment⁴. As representatives of the young-of-the-year population in the region, 2–6 trawl stations each year that had relatively larger catch-per-unit-effort were selected (Supplementary Fig. 1), and 9–20 individuals were randomly selected from each station for otolith analyses (Supplementary Table 1). Age of fish was initially judged by SL (10–15 cm) and later confirmed by the counts of otolith daily increments.

From the eastern North Pacific, archived otoliths of CA sardine captured in cruise surveys and in the pelagic fishery of the Southern California Bight during 1987, 1991–1998, and 2005–2007 were collected. Fish in the size range of 10–16 cm SL were regarded as age-1 individuals born in the previous year, following Takahashi and Checkley⁵⁶. The number of individuals varied between year classes in the range of 4–20 (Supplementary Table 2).

Otolith processing, microstructure and somatic growth analysis

Sagittal otoliths were cleaned to remove the attached tissue in freshwater and then air-dried. Otoliths of JP sardine were embedded in epoxy resin (Petropoxy 154, Burnham Petrographics LLC) on slide-glass, while those of CA were glued to slide-glass using enamel resin and then ground and polished with sandpaper to expose the core. For some otoliths of CA sardine, the polished surface was coated with additional resin to facilitate identification of the daily increment width. Using an otolith measurement system (RATOC System Engineering Co. Ltd.), the number and location of daily increments were examined along the axis in the postrostrum from the core. Although daily increments were clearly observed until the otolith edge for JP sardine, it was difficult to do this for CA sardine probably because they had experienced winter when otolith growth slowed down. Therefore, the rings were counted as far as possible for CA sardine, which typically resulted in more than 150 counts. The first daily increment was assumed to form after 3 days post hatch (dph) for JP and 8 dph for CA sardine following Takahashi et al.²⁶ and Takahashi and Checkley⁵⁶. The otolith radius at each age was calculated by adding all the increment widths up to that age. Standard lengths at each age were back-calculated assuming a linear relationship between otolith radius and standard length using the biological intercept method³⁴ as follows:

where SL_n is the standard length at age n, SL_catch is the standard length at catch, SL_first is the standard length at the age of first daily increment deposition fixed at 5.9 mm for JP sardine and 5.5 mm for CA sardine following the previous studies^26,56, OR_n is the otolith radius at age n, OR_first is the otolith radius at the age of first daily increment deposition, and OR_catch is the otolith radius at catch. Based on rearing experiments of field collected eggs, Lasker⁵⁷ showed the SL of CA sardine at 6–8 dph ranged between 3.8 to 6.5 mm, and Matsuoka and Mitani⁵⁸ showed the total length at 2–4 dph ranged between 4.8 to 6.2 mm, corresponding to 4.7 to 6.1 mm in SL. To deal with these uncertainties regarding the size at the age of first daily increment deposition, we conducted Monte Carlo simulations (10,000 times) to estimate the uncertainties of back-calculated SL, assuming that the initial SLs fall between 3.8 to 6.5 mm for both sardines. Standard deviations of the temporal back-calculated SL at each age were presented as the uncertainty of each SL_n estimation, which varied between 0.51 and 0.73 at the end of larval stage (JP: 45 dph, CA: 60 dph), between 0.34 and 0.64 at the end of early juvenile stage (JP: 75 dph, CA: 90 dph) and between 0.20 and 0.53 at the end of late juvenile stage (JP: 105 dph, CA: 120 dph). These values were significantly smaller than the variability of estimated SL among individuals assuming initial sizes of 5.9 and 5.5 mm for JP and CA sardine, respectively (standard deviations: 4.2, 8.1 and 8.3 in JP sardine and 5.5, 9.1 and 10.3 in CA sardine for the end of larval, early juvenile and late juvenile stages, respectively), suggesting that the back-calculated SL is robust to variations of initial size. Nevertheless, the biological intercept method assumes a constant linear relationship between fish and otolith size within individual⁵⁹, which can vary depending on physiological or environmental conditions^60,61. Therefore, to examine the relationships between temperature and growth, we used both otolith growth, which contains fewer assumptions, and back-calculated somatic growth as growth proxies. Since the use of the two proxies did not show remarkable differences in the relationships between temperature and growth (Supplementary Figs. 11, 12), we mainly used the back-calculated SL in the discussion, which has a more direct ecological implication.

To more generally test whether growth trajectories are different between the western and eastern boundary current systems, otolith growth data of JP and CA sardines were compared with those of sardines in the east to south and west coasts of South Africa. The biological intercept method to back-calculate standard length could not be used in sardine from South Africa because the size at catch was large, some over 20 cm, and otolith radius and standard length were not linearly correlated for fish of this size. Therefore, the otolith radius and increment width were directly used as proxy for size and growth in this comparison, respectively. For visualisation (Fig. 2a), the means of year class mean otolith radii were estimated for JP and CA sardines. For CA sardine, otolith radii at ages were simply averaged within each year class. For JP sardine, to account for the variation in the number of individuals captured at the same station, otolith radii were first averaged within each station, and the station means were averaged within each year, weighted by catch-per-unit-effort. For South African sardine, data of otolith daily increment widths from hatch to 100 dph of 67 adults captured at six stations on the east to south coast (>22^oE), and 51 individuals captured at six stations on the west coast (<20^oE) from 2015 to 2017, published in Sakamoto et al.³⁰, were used.

Micro-milling, powder collection and isotope analysis

As the contamination of enamel resin used for CA sardine otoliths may potentially affect the isotope analyses, we removed the resin by dissolving it in acetone for the otoliths coated with the resin, and the otoliths were air-dried. They were then embedded in Petropoxy 154, which does not create any gas in the reaction with phosphoric acid, and ground again to expose the otolith surface. Test analyses using otoliths of horse mackerel demonstrated that this procedure (i.e., embedding with the resin and removing it with acetone) does not affect the subsequent isotope analyses. Before micromilling, the surface of the polished otolith was cleaned using an ultrasonic cleaner with Milli-Q water and air-dried for a few hours. The otolith portions deposited during hatch–30, 31–45, 46–60, 61–75, 76–90, 91–105, and 106–120 dph for the JP sardine and hatch–30, 31–60, 61–90, 91–120, and 121–150 dph for CA sardine, were milled out sequentially using a high-precision micromilling system Geomill 326. The difference in the temporal resolution was due to the slower growth rates of CA sardine. The milling depth was set to 50 μm for the area near the core and 100 μm for the rest. After each milling, the otolith was observed under a microscope to check for otolith fractions that had cracked outward from other milling area, and these were removed using a needle if present. The milled powders were then collected using a needle and a stainless-steel cup and poured into response vials. After each collection of powder, the otolith was cleaned with an air duster to avoid cross-contamination between the milling paths.

The δ¹⁸O and δ¹³C values of collected otolith powder were determined using the micro-volume analysis system MICAL 3c^62,63 at the National Institute of Technology, Ibaraki College, for the area nearest to the core, and the automatic system DELTA V + GAS Bench for the rest at the Atmosphere and Ocean Research Institute, the University of Tokyo. The otolith powders were reacted with phosphoric acid at 25 °C, and the released CO₂ was purified before being introduced into the mass spectrometer for the MICAL3c system. The response with phosphoric acid was performed at 72 °C for DELTA V. The δ¹⁸O and δ¹³C values were reported in δ-notation relative to the Vienna Pee Dee Belemnite (VPDB) scale and are given as a‰ value. Analytical precisions were better than ±0.10‰ and ±0.17‰ for δ¹⁸O, and better than ±0.10‰ and ±0.15‰ for δ¹³C, respectively. The acid fractionation factor of calcite was used to facilitate comparisons with isotopic values reported in previous studies⁶⁴. Because the difference between the acid fractionation factor of calcite and aragonite is temperature dependent⁶⁵, 0.09‰ was subtracted from the δ¹⁸O value determined by DELTA V to adjust for the different response temperatures.

Estimation of M_oto

Dissolved carbon in fish blood is derived from two sources: dissolved inorganic carbon (DIC) from the ambient water and metabolic carbon released from the respiration of food. Hence, the isotopic composition of dissolved carbon in fish blood, and consequently in the otolith carbonate, is a weighted mean of the δ¹³C values of DIC and metabolic carbon⁶⁶. As the δ¹³C values of DIC are generally higher than those of metabolic carbon, blood δ¹³C values decrease when respiration rates increase and the proportion of metabolic carbon in blood increases³⁸. The proportion of metabolically derived carbon in otolith carbonate (M_oto) from otolith δ¹³C can therefore be used as a proxy for the field metabolic rate of fish. Following Chung et al.³⁸, M_oto was estimated as:

where δ¹³C_diet and δ¹³C_DIC are the δ¹³C values of the diet and DIC in seawater, respectively. The ε term is the total net isotopic fractionation during carbon exchange between DIC and blood and between blood and endolymph in which the otolith is formed, which was set to 0 based on Solomon et al.⁶⁶ The δ¹³C values of the diet of JP and CA sardines were estimated from the δ¹³C of sardine muscle or zooplankton. The δ¹³C values of JP sardine muscle reported in previous studies were in the range of −17.5 to −20.0‰ (−17.5 to −20.0‰ (Yasue et al.⁶⁷: larvae and juveniles in the south coast of Japan); −18.4 ± 0.8‰ (Ohshimo et al.⁶⁸: larvae to adults around Japan)), and δ¹³C values of copepods in the Oyashio and Kuroshio-Oyashio regions were reported as −19.8 to −21.6‰⁶⁹. Considering that the diet-tissue enrichment of δ¹³C is approximately 1.5‰ in marine fish⁷⁰ we assumed that the δ¹³C_diet for JP sardine would fall in the range of −19.0 to −22.0‰. Because the reported δ¹³C values of the muscle of CA sardine in the California Current region were −17.0 to −20.0‰ (−19.8 ± 0.2‰⁷¹; −17.0 ± 0.8‰⁷²; −17.5 to −18.0‰⁷³), the δ¹³C_diet for CA sardine was assumed to be in the range of −18.5 to −21.5‰. The δ¹³C values of DIC in seawater were extracted from the World Ocean Database⁷⁴. As the δ¹³C_DIC is known to show temporal shifts due to the emission of anthropogenic CO₂ that results in a reduction of the δ¹³C of atmospheric CO₂ known as the ¹³C Suess effect⁷⁵, we extracted the δ¹³C_DIC data observed in the Kuroshio-Oyashio system (130–180 °E, 30-45 °N) during 2006–2015 and the California Current system (110–130°W, 30–40°N) during 1986–2006. As the value varied between +0.53 and +1.05‰ in the former system and −0.31 and +2.20‰ in the latter system (Supplementary Fig. 2), we assumed that the δ¹³C_DIC in each region fell within those ranges. To deal with these uncertainties regarding δ¹³C_diet and δ¹³C_DIC, we conducted Monte Carlo simulations 10,000 times to estimate the mean and standard deviation of individual M_oto values at the given life stages. The means of the M_oto values were used in subsequent statistical analyses. Standard deviations of the temporal values were presented as the uncertainty of each M_oto estimation, which varied between 0.01 and 0.03. The relationship between M_oto and oxygen consumption rate has been determined only for juvenile Atlantic cod³⁸, although the relationship can be different between species and life stages especially for the larval stage when seawater DIC may be incorporated through cutaneous exchange rather than through the gills. Therefore, we did not convert M_oto values to the oxygen consumption rate and analysed them directly as a metabolic proxy, and we compared them only within each life stage to investigate relationships between temperature and metabolic rate.

Measurements of seawater δ¹⁸O in the California Current region

As otolith δ¹⁸O in fish is affected by temperature and δ¹⁸O of ambient water, the seawater δ¹⁸O distribution is essential for estimating the temperature from otolith δ¹⁸O. Although the distribution of seawater δ¹⁸O in the western North Pacific has been relatively well studied (Supplementary Fig. 3a), data are limited in the eastern North Pacific especially off southern California where CA sardine grows. We therefore collected surface and sub-surface seawater samples for the isotope analysis that were collected during the 1708SR CalCOFI cruise survey conducted by the California Cooperative Oceanic Fisheries Investigations in August 2017. At every two stations on three line transects extending offshore from the Southern California Bight, seawater samples for δ¹⁸O analysis were taken from 10 m and 50 m depths using CTD-attached Niskin bottles and preserved in sealed glass vials to prevent evaporation. The sampling range covered from the inshore area of the Southern California Bight to the offshore California Current region which was assumed to represent larval and juvenile habitats off Southern California. After membrane-filtration (pore size: 0.45 μm, Toyo Roshi Kaisha, Ltd.), δ¹⁸O values were measured at the National Institute of Technology, Ibaraki College using the Picarro L2130-i system. Data were reported in δ-notation against the VSMOW (Vienna Standard Mean Ocean Water) with a precision better than ±0.05‰.

Seawater δ¹⁸O of the Southern California Bight did not show large variation either horizontally or vertically (Supplementary Fig. 3b, c), ranging between −0.42 and −0.20‰, with mean value of −0.32‰ and a standard deviation of 0.05‰. The salinity of these samples measured by CalCOFI ranged between 33.08 and 33.56, and a significant correlation was detected between seawater δ¹⁸O and salinity as follows:

The extent of potential inter-annual and seasonal variations in seawater δ¹⁸O was analysed based on variations in salinity, as seawater δ¹⁸O is generally correlated with salinity⁷⁶. A total of 14732 salinity measurements of bottle samples in the upper 50 m of CalCOFI cruises in spring and summer during 1986, 1990–1997 and 2004–2006, which horizontally encompassed 117.2–125.7 °W and 29.8–35.1 °N, were extracted from the CalCOFI Hydrographic database (https://www.calcofi.org/ccdata/database.html, accessed on 17^th April, 2020, currently changed to CalCOFI Data Portal (https://calcofi.org/data/)). Within the measurements, 69% (10194 measurements) were in the salinity range that occurred in our δ¹⁸O measured samples, and 95% was in the range between 32.93 and 33.69, corresponding to −0.44 and −0.23‰ in seawater δ¹⁸O based on the regression above (Eq. 3). This suggests that spatial and temporal variations of seawater δ¹⁸O in the main habitat area of CA sardine are limited and likely to fall in the range of −0.32 ± 0.12‰.

Conversion of otolith δ¹⁸O to temperature

Otolith δ¹⁸O is affected by both temperature and δ¹⁸O of the ambient water. In the habitat area of CA sardine, the seawater δ¹⁸O in the habitat area did not show large horizontal or vertical variation (Supplementary Fig. 3b, c), while seawater δ¹⁸O in the habitat area of the JP sardine is known to show considerable variation, ranging from −1‰ to +0.5‰ (Supplementary Fig. 3a). Therefore, the temperatures were calculated using different methods.

As seawater δ¹⁸O is generally correlated with salinity⁷⁶, otolith δ¹⁸O can be regarded as a 2-variable function of temperature and salinity when seawater δ¹⁸O variation cannot be ignored. If seawater temperature and salinity are completely independent of each other, it would be impossible to estimate both from just one otolith δ¹⁸O value, although they are often closely related. Therefore, using the relationship between salinity and temperature, which varies annually and seasonally, estimating both parameters from otolith δ¹⁸O would become possible with a certain range of error. To build formulas to calculate the temperature from otolith δ¹⁸O for each month of each year, the surface layer (< 30 dbar in pressure) temperature and salinity observed by Argo floats in the range of 130–180 °E, 30–45 °N but excluding nearshore areas and the Sea of Japan, from February to October 2006–2010 and 2014–2015, were extracted from the Argo float dataset Advanced automatic QC(AQC) Argo Data ver.1.2a distributed by JAMSTEC⁷⁷ (Supplementary Fig. 4a). The number of observations in each month was approximately 2200 on average, varying from 319 to 5427. For each temperature and salinity pair, the corresponding otolith δ¹⁸O was calculated using the seawater δ¹⁸O-salinity relationship in the Kuroshio-Oyashio system³⁹

and the otolith δ¹⁸O-temperature and seawater δ¹⁸O relationship for the JP sardine³⁷.

The temperature was plotted against the otolith δ¹⁸O for each month, which generally showed a curve-shaped relationship (Supplementary Fig. 4b). A quadratic function was fitted to the plots using the least squares method and was used as the formula to estimate the temperature from otolith δ¹⁸O. The root-mean-square-errors of the formulas can be regarded as proxies for the accuracies of temperature estimation. These were 1.0 °C on average, varying from 0.3 °C to 2.0 °C with a tendency to increase in summer months (Supplementary Fig. 4b, c; Supplementary Table 3). All otolith δ¹⁸O values of the JP sardine were converted to temperature using the formula made for the month in which the median date of each milled area belongs. These analyses were performed using MATLAB R2017a (The MathWorks, Inc., Natick, Massachusetts, United States).

The ambient temperature of CA sardine was calculated using otolith δ¹⁸O-temperature and seawater δ¹⁸O relationship for the JP sardine³⁷ (Eq. 5), with the fixed seawater δ¹⁸O value of −0.32‰. As the seasonal and inter-annual variations of seawater δ¹⁸O were mostly limited to the range of ±0.12‰, the error due to these variations is smaller than 0.7 °C.

Relationship between inter-annual variation of experienced water temperature and environmental indices

To understand how variations in habitat temperature are controlled, inter-annual variations of experienced water temperature were compared to environmental indices. As the distribution during juvenile stages can be significantly different between years due to the changes in dispersal and migration patterns³⁹, detecting the correlation between experienced water temperature in the stage and indices would be difficult. Thus, we focused on the temperature during 1 or 2 months from hatch. Year-class mean hatch dates of JP sardine varied between March and April. As there are no data for hatch date distribution for CA sardine, we assumed that CA sardine hatched in April, which is known as the peak spawning month of CA sardine (e.g., Lo et al.⁷⁸). Based on these assumptions, the mean experienced temperature of JP sardine from hatch to 60 dph was compared to mean PDO⁴¹ values from March to May, while the experienced temperature of CA sardine during 0-30 dph was compared to the PDO during April. The PDO index were downloaded from the webpage of National Centers for Environmental Information (https://www.ncei.noaa.gov/access/monitoring/pdo/, accessed on 10^th August 2022). As we did not find significant correlations for CA sardine, the temperature for CA sardine was compared to the Coastal Upwelling Transport Index⁴² (CUTI) off Southern California (mean of the indices at 33–36 °N). The CUTI index was downloaded from the MIKE JACOX webpage (https://mjacox.com/upwelling-indices/, accessed on 27th October 2020). It should be noted that because the number of data points (7 and 11 for JP and CA sardine respectively) are not large, the relationships observed here may not represent the true relationship between sardine nursery temperature and environmental indices, even if correlation coefficients may be high.

Statistical analyses for the differences of otolith increment widths among regions

To understand the mechanism creating the differences and similarities in otolith growths trajectories among JP, CA, and west coast, and south-east coast SA sardines (Fig. 2a), otolith increment widths during every 10 days between hatch and 100 dph were analysed (Supplementary Fig. 5). A linear mixed-effects model based on the R 4.1.3 and the libraries lmerTest, MuMIn and emmeans was used. Each 10-day increment width (IW) was modelled by two fixed factors (Region and Age) and a random effect (individual fish, Fish.ID) as lmer(IW ∼ Age*Region + (1 | Fish.ID)). Here, age was used as a factor. The diagnostic for the model showed a straight Q-Q plot and the normality of residuals (Supplementary Fig. 6). The pairwise comparison between regions at each age using emmeans showed that there were significant differences in some pairs in most of the age range, although increment widths of JP and the south-east coast SA sardines were not significantly different during 21–50 dph and those of CA and the west coast SA sardines were not significantly different during 31–60 dph (Supplementary Table 4).

Statistical analyses for the differences of M_oto and experienced seawater temperature between JP and CA sardines

To test for the differences in M_oto and experienced seawater temperature between JP and CA sardines, we used a linear mixed-effects model based on the R 4.1.3 and libraries lmerTest, MuMIn and emmeans. The M_oto and temperature were modelled by two fixed factors (Region and Age) and a random effect (individual fish, Fish.ID) as lmer(M_oto∼Age*Region+(1 | Fish.ID)) and lmer(Temperature ∼ Age*Region+(1 | Fish.ID)), respectively. Outliers of M_oto or the temperature detected using the boxplot function were excluded in each analysis. Here, age was used as a factor and values for every 15-day interval for JP sardine were grouped into the age groups for every 30-day interval for the test (0–30, 31–60, 61–90, 91–120 dph). The diagnostics for the models showed mostly straight Q-Q plots and normalities and homogeneities of the residuals (Supplementary Figs. 7, 8). The pairwise comparisons showed that M_oto was significantly higher in JP sardine in all age groups and the experienced temperature was significantly higher in age groups except for 91–120 dph than in CA sardine (Supplementary Tables 5, 6).

Statistical analyses for the effect of temperature on M_oto and M_oto variations

The relationship between the experienced water temperature and M_oto was investigated to test the effect of temperature on metabolic performance. As M_oto was strongly dependent on life stage (Fig. 2b), we used life-stage-averaged data for each individual (JP larva: 0–45 dph, early juvenile: 46–75 dph, late juvenile: 76–105 dph, CA larva: 0–60 dph, early juvenile: 61–90 dph, late juvenile: 91–120 dph). Outliers of M_oto detected using the boxplot function in R 4.1.3 for each stage and region were excluded.

First, to test the difference of M_oto between JP and CA sardines taking the effect of temperature into account, we applied a generalised linear model. Gaussian function family function was used as overall M_oto distribution showed no significant discrepancy from Gaussian distribution (Shapiro-Wilk test, w = 1.00, p > 0.05). Theoretically, the relationship between metabolism and temperature tends to show a linear trend after the metabolic rate is log-transformed⁷⁹. Thus, we applied “identity (data without transformed)” and “log (data transformed)” links to evaluate if model shows a better linearity with data transformation. Based on AIC, however, the result showed M_oto have a better linearity without data transformation (Supplementary Table 7). We, therefore, used “identity” links for the further model selection. Model selection base on AIC was performed for models including temperature, region (JP and CA sardines), life history stages (larvae, early juvenile and late juvenile) and interactions of these factors. The full model including all the interactions had the lowest AIC (Supplementary Table 7). As the diagnostic for the full model showed normality and homogeneity of residuals (Supplementary Fig. 9), we selected this model for interpretation. The CA sardine at the larval stage as the baseline, we found only JP sardine at early and late juvenile stages has relatively higher M_oto values, and the temperature-dependent slope is significantly gentler in JP sardine at early and late juvenile stages (Supplementary Table 8).

Next, the diversity of M_oto across temperature range was assessed to estimate the optimal temperature in each stage. The relationship between the maximum metabolic rate and temperature is known to be parabolic, while that between the standard metabolic rate and temperature is logarithmic^28,79. As the highest field metabolic rate would be constrained by maximum metabolic rate and the lowest field metabolic rate would be close to resting metabolic rate⁴³, fish would have the most diverse metabolic performance at the optimal temperature with the widest aerobic scope. Thus, we modelled the highest and lowest M_oto values in each 1 °C bin using a polynomial regression and a generalised linear model with Gaussian distribution and a log link for the 95th and 5th percentile values of each bin, respectively (Supplementary Fig. 10). The values of the bin that included less than four values were excluded from the regression analyses to reduce the uncertainty caused by under-sampled temperature bins. The gap between the two regression lines was considered as a proxy for the aerobic scope, and the temperature at which the gap reached the maximum was regarded as the optimal temperature.

Statistical analyses for the relationships between temperature and growth

To understand how variation in ambient water temperature affects early life growth of sardines, we compared back-calculated standard length at around the end of the larval stage (hatch–35 mm; JP: 45 dph, CA: 60 dph), the end of the early juvenile stage (35–60 mm; JP: 75 dph, CA: 90 dph), and the end of the late juvenile stage (60–85 mm; JP: 105 dph, CA: 120 dph) and the mean seawater temperature from hatch to the ages. Median of each sampling batch were used as minimal data unit. Pearson’s r and p-values were first calculated for each comparison (Supplementary Table 9). As the relationship between mean temperature and standard length of JP at 75 dph seemed to be dome-shaped rather than linear, we introduced quadratic term of temperature and tested whether the term increased explanatory power using a linear model and stepwise model selection based on AIC. The model selection showed that the full model (Standard length ∼ Temperature² + Temperature) was the best model, and the coefficients of the quadratic and linear terms were both significant (Supplementary Table 10). To account for these multiple tests, we corrected the p-values of the coefficients of the quadratic term in the linear model for JP sardine at 75 dph and of the Pearson’s r for the rest using the Benjamini-Hochberg procedure with α = 0.05, and selected the null hypotheses that could be rejected (Supplementary Table 9). To compare the temperature that allow maximisation of growth rate and optimal temperature derived from the analysis of M_oto for each stage, median somatic growth rate and otolith increment width in each 1 °C bin was calculated together with its 3-window running mean (Supplementary Figs. 11, 12).

Statistical analyses for the relationships between sea surface temperature and survival index

To test whether habitat temperatures during the first 4 months after hatch affect the survival of sardines in the first year of life on a multidecadal scale, satellite-derived sea surface temperature (SST) since 1982 and survival of JP and CA sardines were compared. The log recruitment residuals from Ricker recruitment models (LNRR)¹³, representing early life survivals taking into account the effect of population density, were calculated based on the stock assessment data for JP and CA sardines as follows:

where LNRR_t is the LNRR at year t, R_t is the recruitment of year-class t, S_t is the spawning stock biomass in year t, and a and b are the coefficients of linear regression of ln(R_t/S_t) on S_t. Pearson’s r between the LNRR and the mean SST values from March to June for JP and from April to July for CA sardine was calculated for each grid points in the western and eastern boundaries of the North Pacific basin, derived from a SST product based on satellite and in situ observations⁸⁰ (Global Ocean OSTIA Sea Surface Temperature and Sea Ice Reprocessed (https://resources.marine.copernicus.eu/product-detail/SST_GLO_SST_L4_REP_OBSERVATIONS_010_011/INFORMATION), accessed on 11^th August and 28^th October 2021). The correlations were generally negative and positive in the western and eastern regions, respectively (Supplementary Fig 13a, b). In particular, mean SST values in the area where eggs, larvae and juveniles of JP or CA sardines are mainly found in the months^{26,39,49,56,78,81,82} (dotted areas in Supplementary Fig 13a, b) were compared with LNRR values to test the relationship between habitat temperature and survival in the early life stages (Supplementary Fig 13c). It should be noted that the mean SST values were not significantly correlated with otolith-derived year-class mean temperatures of JP and CA sardines during the larval to late juvenile stages (JP: r = 0.01, p = 0.98, n = 7, CA: r = 0.29, p = 0.38, n = 11), likely due to the short periods analysed, patchy distribution and inter annual variation in larval and juvenile dispersal and migration patterns. Nevertheless, the regions included areas where SST showed weak to significant (p < 0.05) positive correlations with otolith-derived temperatures (Supplementary Fig 13a, b). In addition, mean SST values for the western and eastern regions were significantly correlated with PDO (1982–2016, March to June, r = −0.56, p < 0.001, n = 35) and Coastal Upwelling Transport Index (31–40 °N, 1988–2016, April to July, r = −0.50, p < 0.01, n = 29) in longer time scales, respectively. These are the indices that were able to reasonably explain otolith-derived temperatures during hatch to 30–60 dph (Fig. 2d, e). These suggest that the mean SST values in the regions at least partly reflect the habitat temperature of the sardines and can be used to test the relationship between habitat temperature and survival.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com