https://orcid.org/0000-0002-9182-5854
Abstract
Larval abundances of Atlantic bluefin tuna (ABT) in the Gulf of Mexico are currently utilized to inform future recruitment by providing a proxy for the spawning potential of western ABT stock. Inclusion of interannual variations in larval growth is a key advance needed to translate larval abundance to recruitment success. However, little is known about the drivers of growth variations during the first weeks of life. We sampled patches of western ABT larvae in 3–4 day Lagrangian experiments in May 2017 and 2018, and assessed age and growth rates from sagittal otoliths relative to size categories of zooplankton biomass and larval feeding behaviors from stomach contents. Growth rates were similar, on average, between patches (0.37 versus 0.39 mm d−1) but differed significantly through ontogeny and were correlated with a food limitation index, highlighting the importance of prey availability. Otolith increment widths were larger for postflexion stages in 2018, coincident with high feeding on preferred prey (mainly cladocerans) and presumably higher biomass of more favorable prey type. Faster growth reflected in the otolith microstructures may improve survival during the highly vulnerable larval stages of ABT, with direct implications for recruitment processes.
INTRODUCTION
Atlantic bluefin tuna (ABT, Thunnus thynnus, Linnaeus 1758) is the largest species of the tuna family (Scombridae), reaching sizes up to 4.6 m and 900 kg (Anonymous, 2019). Adults of the western ABT stock migrate large distances from rich North Atlantic feeding grounds to spawn mainly in the oligotrophic waters of the Gulf of Mexico (GoM) between April and June (Block, 2019; Fromentin and Powers, 2005; Rooker et al., 2008), which is believed to put their offspring in an optimal habitat for survival (Bakun, 2013). Spawning appears to be mediated by sea surface temperature above ~24°C beyond the 200-meter shelf margin and along the outer edges of anticyclonic eddies (Domingues et al., 2016; Muhling et al., 2010). The planktonic larvae, restricted to the upper 25 m of the water column, are subsequently challenged to grow rapidly, or perish due to starvation or predation, especially during the most vulnerable first weeks of life (García et al., 2013; Kimura et al., 2010; Shropshire et al., 2021).
Larval ABT habitat has been extensively studied and modeled in the GoM (Alvarez et al., 2021; Domingues et al., 2016; Laiz-Carrión et al., 2015; Lindo-Atichati et al., 2012; Muhling et al., 2010), as well as in the Mediterranean Sea (MED), where the eastern stock spawns (Alvarez-Berastegui et al., 2016; Ingram Jr. et al., 2017; Reglero et al., 2019). In the MED, García et al. (2006) examined DNA and protein ratios, and Ingram Jr. et al. (2017) incorporated environmental parameters into larval condition indices. The influence of habitat quality on larval growth has yet to be examined for the GoM. Higher growth rates generally lead to enhanced larval survival (Bergenius et al., 2002; Kimura et al., 2010). However, growth may be limited by the variable distribution of preferred zooplankton prey in a heterogeneous environment (McGruck, 1986; Shiroza et al., 2021), and in warm oligotrophic systems like the GoM, high metabolic demand can easily lead to starvation under food-limiting conditions (Shropshire et al., 2021).
Otolith microstructure is routinely used in fisheries management and ecological modeling to quantify age and growth (Begg et al., 2005; Campana and Jones, 1992). Specifically, otoliths provide a historical record of larval growth, with increments deposited daily and increment widths (IW) being proportional to somatic growth (Clemmesen, 1994; Gleiber et al., 2020b; Pannella, 1971; Robert et al., 2007). Similar to the early stages of other top pelagic predators like marlin, sailfish and swordfish, somatic growth rates of ABT larvae are highest during the first week of life (0.86 mm d−1 in days post hatch; dph) and decrease in the following week (0.62 mm d−1), but remain rapid compared to most fishes (Malca et al., 2017). Larval daily increment formation has been validated for the closely related Pacific bluefin tuna (Foreman, 1996; Itoh et al., 2000) and inferred for Southern bluefin tuna (Jenkins and Davis, 1990). Elevated temperature has been shown to enhance growth of larval tuna (García et al., 2013; Kimura et al., 2010; Tanaka et al., 2006) and is the main abiotic driver of tuna distribution and recruitment (Alvarez et al., 2021). However, given the sub-tropical conditions of western ABT spawning grounds, prey availability may be more important than temperature for limiting growth (Jenkins et al., 1991; Shropshire et al., 2021; Tanaka et al., 2006).
Quantifying year-to-year differences in larval growth is relevant for annual estimates of ABT spawning biomass based on larval length distributions (Ingram Jr. et al., 2010; Ingram Jr. et al., 2017; Scott et al., 1993) and may also shed light on environmental drivers of larval survival. While substantial temporal variability in mean growth has not been observed in western ABT larvae aged previously from 2000 to 2012 (Malca et al., 2017), those analyses were for uneven sample sizes of specimens averaged over mixed locations and without knowledge of feeding conditions. The pairing of larval diet with otolith biometrics for the conditions experienced in discrete patches of larvae might more usefully reveal successful feeding strategies or point to factors that determine relative quality of feeding and growth environments (Gleiber et al., 2020a). In addition, while ABT larvae are considered successful feeders in the sense of having daytime feeding incidences of up to 98–100% of at least one prey item in their stomach contents (Catalán et al., 2011; Llopiz et al., 2015; Shiroza et al., 2021; Tilley et al., 2016; Uriarte et al., 2019), there is little understanding of how the variability of their growth rates might be affected by their preferences for and availability of specific prey types.
Here, we report age-at-body size estimates for GoM-spawned ABT larvae and compare full otolith growth trajectories between two cohorts, each tracked and sampled for several days in Lagrangian experiments. Larvae were collected as part of the BLOOFINZ-GoM Program (Bluefin Larvae in Oligotrophic Ocean Foodwebs: Investigating Nutrients to Zooplankton in the GoM) during 2017 and 2018 cruises in the peak ABT spawning month of May (Gerard et al., this issue). We couple somatic growth estimates from otolith microstructure to previously determined dietary analyses from the two nursery habitats (Shiroza et al., 2021) and to modeling outputs relating to food limitation and zooplankton biomass structure (Landry et al., 2021; Shropshire et al., 2021). We also examine differences in somatic growth for early and late stage larvae in relation to their shifting dietary preferences. Our study goals were to determine (i) how biometrics change with larval ABT ontogeny, and (ii) to determine the environmental drivers most important for regulating growth. We also test the hypothesis that larval growth would be faster in the 2018 nursery area that was found to be richest in, and with the highest feeding incidences on the preferred prey (cladocerans) by the more advanced postflexion larvae (Shiroza et al., 2021).
METHOD
Larval sampling and processing
General study area for BLOOFINZ in the northeastern Gulf of Mexico. Area for net tow sampling during cycle experiment NF1704-C1 (○) and NF1802-C5 (●).
Temperature (°C), salinity (psu) and fluorescence (volts) profiles for the upper 25 m were determined from CTD casts conducted throughout C1 and C5. Fluorescence (volts) was uncalibrated, but close to mean Chlorophyll a values determined from extracted samples (mg m−3; Selph et al., 2021).
Otolith analyses
The largest otolith pairs (sagittal) were extracted from the cranial cavity of ABT larvae using sharpened glass probes and allowed to air dry. Otoliths were fixed with one drop of mounting medium (Flo-TexxTM), distal side facing up. One randomly selected sagitta was chosen for age estimation. If the selected sagitta was cracked or otherwise damaged, the corresponding sagitta was aged. We examined all otoliths at 400–1000× in immersion oil (transmitted light) with a compound microscope (Zeiss A.1), digital camera and calibrated image analysis software (Image Pro Plus 7 or Image Pro Plus 10). Daily increments were counted along the otolith radius (OR) twice by a single experienced reader. Since most ABT otoliths do not have a clear hatch mark, a minimum distance of 7 μm from the central core was designated as the starting point (García et al., 2013; Jenkins and Davis, 1990), and estimated increments were added along the OR until reaching the first observed increment. The last increment is often obscured by light microscopy artifacts and is less discernable (Campana and Jones, 1992), particularly for larger otoliths. To address this issue in age determination, a terminal increment was added if there was space for one along the OR. Measurement precision was calculated using Chang’s (1982) coefficient of variation (CV) adjusted for age. Hereafter, age refers to days post hatch (dph).
Recent otolith growth (μm) was calculated from the average of the last three completed increments, which has been shown to characterize larval conditions prior to collection (Clemmesen, 1994; Gleiber et al., 2020a; Shulzitski et al., 2015). Incomplete increments were excluded from further analyses.
Larval gut contents and prey availability
Larvae selected for aging were also examined for gut contents as reported in Shiroza et al. (2021). ABT larvae examined in this section were preserved, measured and staged in the same manner as mentioned previously. Briefly, the gut contents from daytime caught larvae were examined and characterized into nine taxonomic categories of prey: Ciliophora, Podonidae, Copepoda nauplii, Calanoida, Corycaeidae, Other Copepoda, Appendicularia, Acanthopterygii fish larvae and other prey that included crustacean zooplankton. Ingested carbon weights (mg C) for these taxa were then estimated from length-dry weight conversion factors in Supplementary Table S1 of Shiroza et al. (2021) for each individual larva. Gut content data also provide estimates of prey size range as a function of larval length, which can be used to estimate total prey biomass from measured zooplankton biomass. For all aged larvae, prey size length ranged from 80 to 914 μm.
During the BLOOFINZ cruises, bulk mesozooplankton biomass (μg C m−3) were measured in multiple size classes (0.2–0.5, 0.5–1, 1–2, 2–5, >5 mm) (Landry and Swalethorp, 2021). Measurements were made for each cycle during day and night hours; however, only day-time tows are considered here because larval ABT are visual feeders. A 1-m diameter ring net (0.2-mm mesh) was towed obliquely through the euphotic zone (100 m at C1, 135 m at C5), and the cod end contents were size-fractioned through nested Nitex screens. Here, we examine mesozooplankton biomass from 0.2 to 1 mm given the overlap with prey size range of aged larvae. A bongo-20 was towed in the upper 25 m of the water column with 0.05 and 0.2 -mm nets to record bulk biomass of smaller zooplankton, and identify, count and measure prey in both nurseries (see Shiroza et al., 2021). Together, both net samplings provide measurements of zooplankton biomass from 0.05 to 1 mm which fully encompasses the planktivorous prey fields for all aged larvae (see Table IV).
To provide an additional estimate of the small zooplankton biomass (SZB, 0.002–0.2 mm range), observed mesozooplankton biomass in the 0.2–1.0 mm range (herein referred to large zooplankton biomass, LZB) was adjusted to the top 25 m of the water column and scaled using the ratio of SZB to LZB estimated by a three-dimensional biogeochemical model NEMURO-GoM (Shropshire et al., 2020). The NEMURO-GoM model was designed to simulate zooplankton biomass distribution in the GoM and has been extensively validated using remote and in situ measurements including over two decades of zooplankton biomass measurements collected by the Southeast Area Monitoring and Assessment Program (SEAMAP). The ratio of LZB to SZB was determined from daily climatologies generated by the model at each sample location and day of the year over the 20-year simulation (1993–2012). An advantage of the LZB and SZB was that they could be estimated for each ABT larval sampling location, as bulk zooplankton net tows did not accompany every larval net tow. For more information on NEMURO-GoM see Shropshire et al. (2020). Two additional variables were calculated from SZB and LZB that estimated the respective prey biomasses defined as a function of ABT-larval length for small prey biomass (SPB) and for large prey biomass (LPB), respectively (see Table IV).
Lastly, to evaluate prey availability, a Food Limitation Index (FLI, Shropshire et al., 2021) was computed for all aged larvae. The index is defined as the ratio of metabolic requirement to assimilated ingestion where values >1 indicate prey limitation. The ingestion formulation includes many terms but is primarily a function of sensory radius and prey biomass. Prey biomass is estimated using in situ LZB and estimated SZB along with estimates of prey size range as a function of larval length. Sensory radius is modeled as a function of prey size and larval length using a recently determined anatomical relationship for visual acuity (Hilder et al., 2019). Metabolic requirement is estimated primarily as a function of larvae age and temperature. Specifically, metabolic requirement is derived from the first derivative of the age to weight relationship with assumptions regarding gross growth and assimilation efficiencies. In situ temperature values were obtained from the nearest sampling station for each cycle and are used to scale estimates of metabolic requirement. For more information on FLI formulations see Shropshire et al. (2021). Collectively, FLI values along with estimates of prey biomass and in situ measurements of mesozooplankton biomass are the metrics used to investigate drivers of differences in larval growth potential.
Data analysis
Data and statistical tests were carried out in R 4.0.3 (R Core Team, 2020). Least squares regressions were calculated for best fits of the following metrics to age: length, OR, body depth (mm), dry weight (mg), mean IW (μm), recent growth (μm) and residuals for age-at-length and OR at age. Analyses of covariance (ANCOVA) were carried out using age as a continuous covariate and log-transformed biometric variables when required to meet normality assumptions; overall growth was examined comparing the slopes between C1 and C5. Finally, recent growth between stages and cycles were analyzed in separate ANCOVAs with age as a continuous covariate to test for any significant interactions.
To test for differences in larval growth between C1 and C5, we analyzed the otolith microstructure following full growth trajectories with a linear mixed model using the nlme package (Pinheiro et al., 2014). Following the design outlined in Swalethorp et al. (2016) and Malanski et al. (2020), the model was fit to the data with otolith IW as the dependent variable and cycle as independent. Increment number was nested by individual larvae and included as a random effect, cycle included as a fixed effect. We split otolith growth trajectories in two groups using 7 dph as a cutoff because this is when the IWs for C1 and C5 intersect. The random effects were applied to both model intercepts and slopes. To correct for autocorrelation and nonindependence of our consecutive otolith IW measurements (Campana, 1996; Chambers and Miller, 1995), the model was refitted with an autocorrelation structure with increment number as the continuous time covariate using the corCAR1 function (Fox and Weisberg, 2015). Since larvae differed in number of growth increments (unbalanced design), we used the maximum likelihood to estimate slopes and model significance (Plant, 2012).
Modeling approach
We examined the effect of variables (diet, prey availability, and temperature) on recent otolith IW using generalized additive models (GAMs) with a Gaussian distribution. GAMs are nonparametric, generalized linear models with flexibility to handle both linear and complex relationships between the explanatory and response variables within the same environment (Wood, 2004, 2017). All models were constructed using the mgcv library in R. Variables included in model selection are listed in Table II. The model selected recent IW as the dependent variable for the subset of larvae (n = 139) with all considered dietary metrics measured or estimated.
The dietary metrics considered are the abundance and corresponding carbon (C) mass of all ingested prey and of preferred ingested prey. Preferred prey are the sum of the two prey categories, copepod nauplii and podonids, most preferred across all larval sizes. Concurrent, zooplankton prey biomass (small and large), FLI, and hydrographic variables were also included.
To account for potential correlations between explanatory variables, correlations (Spearman’s correlation matrix, ρ > 0.6) between all potential explanatory variables were identified and strongly correlated variables were modeled against the response in single-variable GAMs. The Akaike Information Criterion (AIC) (Akaike, 1974) of the single-variable GAMs were compared between correlated variable pairs, and the variable with the lowest AIC was included in the final model selection process. After the set of noncorrelated explanatory variables was identified, overall multicollinearity was assessed using the variance inflation factor (VIF) with 3 as cutoff. Smoothing functions were applied to continuous predictor variables restricted to 4 knots to avoid overfitting.
To select a final model, the restricted maximum likelihood (REML) method was used as it applies a double penalty to smooth terms and allows for removal of variables with minor predictive values (Marra and Wood, 2011). Model diagnostics and residuals were checked for potential deviations from normality and homogeneity of variance.
RESULTS
In total, 30 and 38 bongo-90 tows were taken during cycles C1 and C5, and ABT larvae were aged from 18 to 20 stations, respectively (Fig. 1). Aged larvae represented similar times during the peak month of the ABT spawning season in both years. Back-calculating spawning from observed ages indicates that adult ABT spawning events occurred between 24 April through 4 May 2017 for C1, and 29 April through 13 May 2018 for C5. Four larvae between 13 and 15 dph, one from C1 and 3 from C5, were found to be piscivorous, the first time that larval ABT piscivory has been established from direct observations.
The most salient difference between the two cycles was that larval abundances at C5 were almost 14 times higher than C1. Despite very narrow ranges, both temperature (24.19–25.95°C at 25-m depth) and salinity (35.6–36.43 psu) differed significantly between cycles (Wilcox, P < 0.001), with C1 having slightly cooler and more saline conditions consistent with an offshore environment. Fluorescence (also measured at 25-m depth) did not differ between sites (t-test = −2.921, df = 14.164, P = 0.995) although C1 had overall lower values (see abiotic variables in Table IV).
Shrinkage for freshly measured ABT larvae was measured for the first time and found to be 9.24% ± 3.5 (average ± SD) when preserved in 95% ethanol (Equation 1). This average includes larvae caught during both surveys, within C1 and C5 as well as outside the cycles to extend the fresh size range (2.45–9.64 mm SL) of larvae collected.
Biometrics and growth
In total, 198 daytime ABT larvae were aged from both cycles, with a subset of 158 of the same larvae examined for stomach contents. For C1, the 98 aged larvae ranged between 3 and 17 dph with one larvae at 19 dph. For C5, the 99 aged larvae were 3–16 dph. CVs (mean ± SD) for aged larvae was 3.28% ± 3.42 for C1 and 3.84% ± 4.47 for C5. We excluded four C5 larvae from further analysis that did not meet aging precision criteria (CV > 10%, E. Malca unpublished). Developmental stage distribution of the aged larvae were 20 preflexion, 30 flexion and 49 postflexion larvae for C1 and 32 preflexion, 30 flexion and 37 postflexion for C5 (Fig. 2a,b respectively). Preflexion larvae ranged from 3 to 10 dph (6.9 ± 2.0 dph). The youngest larva to reach flexion was 7 dph and the oldest flexion larva was 14 dph (10.2 ± 1.6 dph). Postflexion larvae were 9 dph and older (11.5 ± 1.0 dph).
Histogram for ABT larvae measured from (a) C1 and (b) C5 in the GoM. The boxplot indicates larval developmental stages preflexion, flexion and postflexion with corresponding larval length (SL, mm) from ethanol-preserved specimens (gray). Daytime aged ABT larvae with gut contents examined are also indicated.
Larval biometrics from the two cycles overlapped both with respect to length and OR (Fig. 3a, b). We found no difference in mean larval somatic growth rates between C1 and C5 (0.358 vs 0.368 mm SL d−1, Table I). Body depth was measured only on undamaged larvae (n = 157) and is reported here for western ABT larvae for the first time. On average, C1 larvae were slightly heavier, had greater body depth (Supplementary Fig. S1 a,b) and were older compared to C5; however, older C5 larvae (>12 dph) were slightly heavier than C1 larvae of the same age. Least squares regressions for length relative to age and for OR-at-age residuals relative to length-at-age residuals did not differ between cycles (Supplementary Fig. S1c, Table I), meaning that the body length relative to otolith growth patterns were consistent between C1 and C5 larvae, on average. Body depth, weight and OR-at-age all showed exponential relationships (Fig. 3b, Supplementary Fig. S1) and had linear fits (significant positive slopes) when log-transformed (Table I). Although C5 larvae appeared thinner compared to C1, neither body depth, weight or OR differed significantly between cycles for the larvae examined (Table I).
Upper panel provides ANCOVA summary for larval ABT biometrics
| ANCOVA | F | R2 | P |
|---|---|---|---|
| Overall somatic growth: length ~ age | 959.4 | 0.830 | <0.001* |
| Cycle C1 vs C5 | 0.411 | 0.829 | 0.5223 |
| Overall otolith growth: ORLog ~ age | 1510.8 | 0.885 | <0.001* |
| Cycle C1 vs C5 | 0.5617 | 0.885 | 0.455 |
| ORLog at age residuals and length ~ age residuals | 314.02 | 0.614 | <0.001* |
| Cycle C1 vs C5 | 0.0604 | 0.610 | 0.806 |
| WeightLog ~ age | 535.77 | 0.738 | <0.001* |
| Cycle C1 vs C5 | 4.460 | 0.738 | 0.036 |
| DepthLog at age | 648.12 | 0.805 | <0.001* |
| Cycle C1 vs C5 | 0.367 | 0.803 | 0.545 |
| Linear Mixed Model | |||
| Increment width early growth, C1 (102) vs C5 (103) | 17.056 | NA | <0.001* |
| Increment width late growth, C1 (93) vs C5 (64) | 4.336 | NA | 0.039* |
| ANCOVA | F | R2 | P |
|---|---|---|---|
| Overall somatic growth: length ~ age | 959.4 | 0.830 | <0.001* |
| Cycle C1 vs C5 | 0.411 | 0.829 | 0.5223 |
| Overall otolith growth: ORLog ~ age | 1510.8 | 0.885 | <0.001* |
| Cycle C1 vs C5 | 0.5617 | 0.885 | 0.455 |
| ORLog at age residuals and length ~ age residuals | 314.02 | 0.614 | <0.001* |
| Cycle C1 vs C5 | 0.0604 | 0.610 | 0.806 |
| WeightLog ~ age | 535.77 | 0.738 | <0.001* |
| Cycle C1 vs C5 | 4.460 | 0.738 | 0.036 |
| DepthLog at age | 648.12 | 0.805 | <0.001* |
| Cycle C1 vs C5 | 0.367 | 0.803 | 0.545 |
| Linear Mixed Model | |||
| Increment width early growth, C1 (102) vs C5 (103) | 17.056 | NA | <0.001* |
| Increment width late growth, C1 (93) vs C5 (64) | 4.336 | NA | 0.039* |
Lower panel shows a linear mixed model results considering full otolith increment width history from all larvae. Asterisk indicates significant results at the 0.05 level. For the linear mixed model, early growth are 1–7 dph, late growth 8–15 dph. Number inside parenthesis indicates number of ABT considered for each cycle.
Upper panel provides ANCOVA summary for larval ABT biometrics
| ANCOVA | F | R2 | P |
|---|---|---|---|
| Overall somatic growth: length ~ age | 959.4 | 0.830 | <0.001* |
| Cycle C1 vs C5 | 0.411 | 0.829 | 0.5223 |
| Overall otolith growth: ORLog ~ age | 1510.8 | 0.885 | <0.001* |
| Cycle C1 vs C5 | 0.5617 | 0.885 | 0.455 |
| ORLog at age residuals and length ~ age residuals | 314.02 | 0.614 | <0.001* |
| Cycle C1 vs C5 | 0.0604 | 0.610 | 0.806 |
| WeightLog ~ age | 535.77 | 0.738 | <0.001* |
| Cycle C1 vs C5 | 4.460 | 0.738 | 0.036 |
| DepthLog at age | 648.12 | 0.805 | <0.001* |
| Cycle C1 vs C5 | 0.367 | 0.803 | 0.545 |
| Linear Mixed Model | |||
| Increment width early growth, C1 (102) vs C5 (103) | 17.056 | NA | <0.001* |
| Increment width late growth, C1 (93) vs C5 (64) | 4.336 | NA | 0.039* |
| ANCOVA | F | R2 | P |
|---|---|---|---|
| Overall somatic growth: length ~ age | 959.4 | 0.830 | <0.001* |
| Cycle C1 vs C5 | 0.411 | 0.829 | 0.5223 |
| Overall otolith growth: ORLog ~ age | 1510.8 | 0.885 | <0.001* |
| Cycle C1 vs C5 | 0.5617 | 0.885 | 0.455 |
| ORLog at age residuals and length ~ age residuals | 314.02 | 0.614 | <0.001* |
| Cycle C1 vs C5 | 0.0604 | 0.610 | 0.806 |
| WeightLog ~ age | 535.77 | 0.738 | <0.001* |
| Cycle C1 vs C5 | 4.460 | 0.738 | 0.036 |
| DepthLog at age | 648.12 | 0.805 | <0.001* |
| Cycle C1 vs C5 | 0.367 | 0.803 | 0.545 |
| Linear Mixed Model | |||
| Increment width early growth, C1 (102) vs C5 (103) | 17.056 | NA | <0.001* |
| Increment width late growth, C1 (93) vs C5 (64) | 4.336 | NA | 0.039* |
Lower panel shows a linear mixed model results considering full otolith increment width history from all larvae. Asterisk indicates significant results at the 0.05 level. For the linear mixed model, early growth are 1–7 dph, late growth 8–15 dph. Number inside parenthesis indicates number of ABT considered for each cycle.
Sagittal growth trajectories examined using linear mixed model, showed that average growth (assessed as IW) was significantly larger in early (1–7 dph) stage larvae for C1 compared to C5 (Fig. 4, Table I). Thereafter, the relationship shifted, with IW being significantly larger in 8–15 dph larvae in C5. Similarly, recent IW (last three increments) also diverged around 7 dph (not shown).
Explained growth
We explore available biotic and abiotic variables that could regulate growth using a GAM approach (Supplementary Fig. S3 a–d). Since IW were the only growth estimate found to differ between cycles (Table I), and since differences in the possible explaining variables listed in Table II were generally greater between than within cycles, the GAM analysis was carried out using recent IW (average of the three most recent IWs). Recent IW-at-age was highly correlated to both length-at-age and body depth-at-age underlining its appropriateness as a historical and recent growth metric (Supplementary Fig. S2a,b). Because diet was an important variable in our analyses, four larvae with damaged or empty guts were not included in this analysis. Consequently, diet and age were analyzed for 59 and 95 daytime-collected and aged larvae from C1 and C5, respectively.
Names and descriptions of variables included in the GAM selection process and their data sources
| Variable name | Variable definition | Data sources |
|---|---|---|
| Dietary metrics | ||
| Ingested prey | Sum of all prey in larval gut | Shiroza et al. (2021) |
| Ingested preferred prey | Sum of prey preferred by all larval sizes in gut (copepod nauplii and podonids) | |
| Ingested prey C | Sum of all prey carbon (C) weights in gut | |
| Ingested preferred prey C | Sum of preferred prey C weights in gut | |
| Prey habitat | ||
| Large Zooplankton Biomass (LZB) | Mesozooplankton biomass (0.2–1 mm, μg C m−3) from upper 25 m | Field estimates Landry et al. (2021) |
| Large Prey Biomass (LPB) | Portion of LZB within larval-specific prey size range defined as a function of larval length SL (mm) | Field estimates Landry et al. (2021). SL from Shiroza et al. (2021) |
| Small Zooplankton Biomass (SZB) | Microzooplankton biomass (0.002–0.2 mm, μg C m−3) estimated from measured LZB multiplied by the ratio of SZB to LZB estimated from a biogeochemical model | Field estimates and Shropshire et al. (2021) |
| Small Prey Biomass (SPB) | Portion of SZB within larval specific prey size range defined as a function of larval length SL (mm) | Field estimates and Shropshire et al. (2021). SL from Shiroza et al. (2021) |
| Food Limitation Index (FLI) | Ratio of metabolic requirement to assimilated ingestion, values > 1 indicate prey-limited habitat | Field estimates and Shropshire et al. (2021) |
| Abiotic variables | ||
| Temperature | Temperature °C, 0–25 m | Gerard et al. (this issue) |
| Salinity | Salinity psu, 0–25 m | |
| Fluorescence | Fluorescence v, 0–25 m | |
| Cycle | 2 levels: C1, C5 | |
| Variable name | Variable definition | Data sources |
|---|---|---|
| Dietary metrics | ||
| Ingested prey | Sum of all prey in larval gut | Shiroza et al. (2021) |
| Ingested preferred prey | Sum of prey preferred by all larval sizes in gut (copepod nauplii and podonids) | |
| Ingested prey C | Sum of all prey carbon (C) weights in gut | |
| Ingested preferred prey C | Sum of preferred prey C weights in gut | |
| Prey habitat | ||
| Large Zooplankton Biomass (LZB) | Mesozooplankton biomass (0.2–1 mm, μg C m−3) from upper 25 m | Field estimates Landry et al. (2021) |
| Large Prey Biomass (LPB) | Portion of LZB within larval-specific prey size range defined as a function of larval length SL (mm) | Field estimates Landry et al. (2021). SL from Shiroza et al. (2021) |
| Small Zooplankton Biomass (SZB) | Microzooplankton biomass (0.002–0.2 mm, μg C m−3) estimated from measured LZB multiplied by the ratio of SZB to LZB estimated from a biogeochemical model | Field estimates and Shropshire et al. (2021) |
| Small Prey Biomass (SPB) | Portion of SZB within larval specific prey size range defined as a function of larval length SL (mm) | Field estimates and Shropshire et al. (2021). SL from Shiroza et al. (2021) |
| Food Limitation Index (FLI) | Ratio of metabolic requirement to assimilated ingestion, values > 1 indicate prey-limited habitat | Field estimates and Shropshire et al. (2021) |
| Abiotic variables | ||
| Temperature | Temperature °C, 0–25 m | Gerard et al. (this issue) |
| Salinity | Salinity psu, 0–25 m | |
| Fluorescence | Fluorescence v, 0–25 m | |
| Cycle | 2 levels: C1, C5 | |
Names and descriptions of variables included in the GAM selection process and their data sources
| Variable name | Variable definition | Data sources |
|---|---|---|
| Dietary metrics | ||
| Ingested prey | Sum of all prey in larval gut | Shiroza et al. (2021) |
| Ingested preferred prey | Sum of prey preferred by all larval sizes in gut (copepod nauplii and podonids) | |
| Ingested prey C | Sum of all prey carbon (C) weights in gut | |
| Ingested preferred prey C | Sum of preferred prey C weights in gut | |
| Prey habitat | ||
| Large Zooplankton Biomass (LZB) | Mesozooplankton biomass (0.2–1 mm, μg C m−3) from upper 25 m | Field estimates Landry et al. (2021) |
| Large Prey Biomass (LPB) | Portion of LZB within larval-specific prey size range defined as a function of larval length SL (mm) | Field estimates Landry et al. (2021). SL from Shiroza et al. (2021) |
| Small Zooplankton Biomass (SZB) | Microzooplankton biomass (0.002–0.2 mm, μg C m−3) estimated from measured LZB multiplied by the ratio of SZB to LZB estimated from a biogeochemical model | Field estimates and Shropshire et al. (2021) |
| Small Prey Biomass (SPB) | Portion of SZB within larval specific prey size range defined as a function of larval length SL (mm) | Field estimates and Shropshire et al. (2021). SL from Shiroza et al. (2021) |
| Food Limitation Index (FLI) | Ratio of metabolic requirement to assimilated ingestion, values > 1 indicate prey-limited habitat | Field estimates and Shropshire et al. (2021) |
| Abiotic variables | ||
| Temperature | Temperature °C, 0–25 m | Gerard et al. (this issue) |
| Salinity | Salinity psu, 0–25 m | |
| Fluorescence | Fluorescence v, 0–25 m | |
| Cycle | 2 levels: C1, C5 | |
| Variable name | Variable definition | Data sources |
|---|---|---|
| Dietary metrics | ||
| Ingested prey | Sum of all prey in larval gut | Shiroza et al. (2021) |
| Ingested preferred prey | Sum of prey preferred by all larval sizes in gut (copepod nauplii and podonids) | |
| Ingested prey C | Sum of all prey carbon (C) weights in gut | |
| Ingested preferred prey C | Sum of preferred prey C weights in gut | |
| Prey habitat | ||
| Large Zooplankton Biomass (LZB) | Mesozooplankton biomass (0.2–1 mm, μg C m−3) from upper 25 m | Field estimates Landry et al. (2021) |
| Large Prey Biomass (LPB) | Portion of LZB within larval-specific prey size range defined as a function of larval length SL (mm) | Field estimates Landry et al. (2021). SL from Shiroza et al. (2021) |
| Small Zooplankton Biomass (SZB) | Microzooplankton biomass (0.002–0.2 mm, μg C m−3) estimated from measured LZB multiplied by the ratio of SZB to LZB estimated from a biogeochemical model | Field estimates and Shropshire et al. (2021) |
| Small Prey Biomass (SPB) | Portion of SZB within larval specific prey size range defined as a function of larval length SL (mm) | Field estimates and Shropshire et al. (2021). SL from Shiroza et al. (2021) |
| Food Limitation Index (FLI) | Ratio of metabolic requirement to assimilated ingestion, values > 1 indicate prey-limited habitat | Field estimates and Shropshire et al. (2021) |
| Abiotic variables | ||
| Temperature | Temperature °C, 0–25 m | Gerard et al. (this issue) |
| Salinity | Salinity psu, 0–25 m | |
| Fluorescence | Fluorescence v, 0–25 m | |
| Cycle | 2 levels: C1, C5 | |
The best model fit was achieved using the four variables: SPB, FLI, ingested preferred prey C and average temperature which explained 44.3% of the recent IW variance (Fig. 5, Supplementary Fig. S3, Table III). The most important explaining variable was estimated SPB, which surprisingly decreased with increasing recent IW, suggesting that other size categories could be more important. Indeed, although estimated LPB fell out of the model, that was mainly because it was positively correlated to FLI (r = −0.59), which did a better job of explaining the residual recent IW variance. The second most important variable, FLI, indicated larval growth rate was sensitive to food limitation. In addition, ingested preferred prey C also explained a significant amount of the residual variance, showing that the larvae grew faster with more of the preferred prey ingested. Interestingly, the abundances of preferred prey ingested or as well as total ingested prey C and abundances did poorly at explaining recent IW variance. Although in situ temperature was also significant, the positive relationship to recent IW was weak in comparison to the other variables. Both fluorescence and salinity were correlated with temperature and fell out of the model.
GAM statistics summary for recent IW
| Parametric coefficients | Estimate | Standard error | P value |
|---|---|---|---|
| Intercept | −0.03 | 0.06 | 0.634 |
| Smooth Function | ΔDE (%) | Edf | P value |
| Small prey biomass | 14.8 | 1.880 | <0.001* |
| Food limitation index | 13.4 | 2.390 | <0.001* |
| Ingested preferred prey C | 10.4 | 1.870 | <0.001* |
| Temperature | 4.3 | 1.658 | <0.020* |
| Parametric coefficients | Estimate | Standard error | P value |
|---|---|---|---|
| Intercept | −0.03 | 0.06 | 0.634 |
| Smooth Function | ΔDE (%) | Edf | P value |
| Small prey biomass | 14.8 | 1.880 | <0.001* |
| Food limitation index | 13.4 | 2.390 | <0.001* |
| Ingested preferred prey C | 10.4 | 1.870 | <0.001* |
| Temperature | 4.3 | 1.658 | <0.020* |
Top panel indicates the effect of the parametric coefficients on recent IW. Lower panel shows the estimated significance levels of the smooth functions; ΔDE is the loss in percent deviance explained caused by dropping the variable, “edf” is the estimated degrees of freedom for smooth terms. An asterisk (*) denotes statistical significance (α = 0.001).
GAM statistics summary for recent IW
| Parametric coefficients | Estimate | Standard error | P value |
|---|---|---|---|
| Intercept | −0.03 | 0.06 | 0.634 |
| Smooth Function | ΔDE (%) | Edf | P value |
| Small prey biomass | 14.8 | 1.880 | <0.001* |
| Food limitation index | 13.4 | 2.390 | <0.001* |
| Ingested preferred prey C | 10.4 | 1.870 | <0.001* |
| Temperature | 4.3 | 1.658 | <0.020* |
| Parametric coefficients | Estimate | Standard error | P value |
|---|---|---|---|
| Intercept | −0.03 | 0.06 | 0.634 |
| Smooth Function | ΔDE (%) | Edf | P value |
| Small prey biomass | 14.8 | 1.880 | <0.001* |
| Food limitation index | 13.4 | 2.390 | <0.001* |
| Ingested preferred prey C | 10.4 | 1.870 | <0.001* |
| Temperature | 4.3 | 1.658 | <0.020* |
Top panel indicates the effect of the parametric coefficients on recent IW. Lower panel shows the estimated significance levels of the smooth functions; ΔDE is the loss in percent deviance explained caused by dropping the variable, “edf” is the estimated degrees of freedom for smooth terms. An asterisk (*) denotes statistical significance (α = 0.001).
Environmental cycle differences
Since greater zooplankton or prey biomass could not be directly linked to faster growth in the GAM we indirectly compared in situ zooplankton availability between the two cycles. The mesozooplankton biomass (sizes 0.05–1.0 mm) for C1 and C5 were not significantly different (Table IV). Overall, C5 had more mesozooplankton biomass in the 0.2–1.0 mm range, however C1 had a greater biomass of smaller mesozooplankton in the 0.05–0.2 mm range. Since much of this total biomass is not part of ABT larval diet, we more closely examined the in situ availability of specific taxa that Shiroza et al. (2021) identified as the preferred ABT larval prey. C5 had consistently higher C biomass of the preferred taxa Podonids and Calanoid copepods (Table IV). Although most Appendicularia and Copepod nauplii may be too small to be efficiently caught by the 0.200 mm mesh sized net used, and did not differ significantly between cycles, C5 did contain more of them as well. In addition, the majority of C5 larvae (61%) experienced warmer temperatures (>25°C) which could also allow larvae to grow wider increments reflecting faster growth.
Daytime prey group biomass from cycles C1 and C5 collected using 0.055 and 0.2 mm mesh size plankton nets
| Cycle | Gear | Reference | t test | P | ||
|---|---|---|---|---|---|---|
| In situ bulk zooplankton sizes (mm) | C1 | C5 | ||||
| μg C m−3 ± SD | ||||||
| 0.05–0.2 | 643.77 ± 156.0 | 421.85 ± 137.26 | Bongo-20 (0.05 μm mesh), 0–25 m depth | In situ | 1.64 | 0.18 |
| 0.2–0.5 | 377.8 ± 116.8 | 763.5 ± 443.5 | Ring net (0.2 mm mesh), 0–100/135 m depth | Landry and Swalethorp, 2021 | −1.698 | 0.150 |
| 0.5–1.0 | 637.1 ± 363.4 | 729.3 ± 623.9 | −0.226 | 0.830 | ||
| In situ preferred taxa > 0.2 mm | μg C m−3 ± SD | |||||
| Copepod nauplii | 0.3 ± 0.3 | 0.8 ± 0.6 | −1.326 | 0.242 | ||
| CladoceraCalanoid | 1.3 ± 1.9127.7 ± 30.4 | 48.5 ± 19.5376.9 ± 92.8 | Bongo-20 (0.2 mm mesh), 0–25 m depth | Shiroza et al., (2021) | −4.083−4.386 | 0.01*0.007* |
| Appendicularia | 0.6 ± 0.4 | 0.9 ± 0.5 | −0.696 | 0.518 | ||
| Abiotic variables | Mean ± SD | |||||
| Temperature, °CSalinity, psu | 24.49 ± 0.1336.40 ± 0.01 | 25.23 ± 0.4435.90 ± 0.27 | CTD, 0-25 m | Gerard et al., (this issue) | 1**180** | <0.001*<0.001* |
| Fluorescence, volt | 0.07 ± 0.002 | 0.07 ± 0.002 | Selph, et al., (2021) | 2.921 | 0.995 | |
| Cycle | Gear | Reference | t test | P | ||
|---|---|---|---|---|---|---|
| In situ bulk zooplankton sizes (mm) | C1 | C5 | ||||
| μg C m−3 ± SD | ||||||
| 0.05–0.2 | 643.77 ± 156.0 | 421.85 ± 137.26 | Bongo-20 (0.05 μm mesh), 0–25 m depth | In situ | 1.64 | 0.18 |
| 0.2–0.5 | 377.8 ± 116.8 | 763.5 ± 443.5 | Ring net (0.2 mm mesh), 0–100/135 m depth | Landry and Swalethorp, 2021 | −1.698 | 0.150 |
| 0.5–1.0 | 637.1 ± 363.4 | 729.3 ± 623.9 | −0.226 | 0.830 | ||
| In situ preferred taxa > 0.2 mm | μg C m−3 ± SD | |||||
| Copepod nauplii | 0.3 ± 0.3 | 0.8 ± 0.6 | −1.326 | 0.242 | ||
| CladoceraCalanoid | 1.3 ± 1.9127.7 ± 30.4 | 48.5 ± 19.5376.9 ± 92.8 | Bongo-20 (0.2 mm mesh), 0–25 m depth | Shiroza et al., (2021) | −4.083−4.386 | 0.01*0.007* |
| Appendicularia | 0.6 ± 0.4 | 0.9 ± 0.5 | −0.696 | 0.518 | ||
| Abiotic variables | Mean ± SD | |||||
| Temperature, °CSalinity, psu | 24.49 ± 0.1336.40 ± 0.01 | 25.23 ± 0.4435.90 ± 0.27 | CTD, 0-25 m | Gerard et al., (this issue) | 1**180** | <0.001*<0.001* |
| Fluorescence, volt | 0.07 ± 0.002 | 0.07 ± 0.002 | Selph, et al., (2021) | 2.921 | 0.995 | |
Ring net sampled 0–100 m during C1 and from 0 to 133 m during C5. Bongo-20 sampled from 0 to 25 m. Of the taxa-specific groups, only those that were positively selected by ABT larvae were considered (Shiroza et al., 2021, Fig. 10). Please note that most copepod nauplii were small and not efficiently captured by the 0.2 mm mesh nets. Bottom panel indicates mean values for temperature, salinity and fluorescence for the upper 25 m determined from CTD casts. Fluorescence (volt) values are uncalibrated but close to mean values determined from extracted samples (mg m−3; Selph et al., 2021). Asterisk (*) indicates significant results at the 0.05 level for t-test; double asterisk (**) indicates Wilcoxon test.
Daytime prey group biomass from cycles C1 and C5 collected using 0.055 and 0.2 mm mesh size plankton nets
| Cycle | Gear | Reference | t test | P | ||
|---|---|---|---|---|---|---|
| In situ bulk zooplankton sizes (mm) | C1 | C5 | ||||
| μg C m−3 ± SD | ||||||
| 0.05–0.2 | 643.77 ± 156.0 | 421.85 ± 137.26 | Bongo-20 (0.05 μm mesh), 0–25 m depth | In situ | 1.64 | 0.18 |
| 0.2–0.5 | 377.8 ± 116.8 | 763.5 ± 443.5 | Ring net (0.2 mm mesh), 0–100/135 m depth | Landry and Swalethorp, 2021 | −1.698 | 0.150 |
| 0.5–1.0 | 637.1 ± 363.4 | 729.3 ± 623.9 | −0.226 | 0.830 | ||
| In situ preferred taxa > 0.2 mm | μg C m−3 ± SD | |||||
| Copepod nauplii | 0.3 ± 0.3 | 0.8 ± 0.6 | −1.326 | 0.242 | ||
| CladoceraCalanoid | 1.3 ± 1.9127.7 ± 30.4 | 48.5 ± 19.5376.9 ± 92.8 | Bongo-20 (0.2 mm mesh), 0–25 m depth | Shiroza et al., (2021) | −4.083−4.386 | 0.01*0.007* |
| Appendicularia | 0.6 ± 0.4 | 0.9 ± 0.5 | −0.696 | 0.518 | ||
| Abiotic variables | Mean ± SD | |||||
| Temperature, °CSalinity, psu | 24.49 ± 0.1336.40 ± 0.01 | 25.23 ± 0.4435.90 ± 0.27 | CTD, 0-25 m | Gerard et al., (this issue) | 1**180** | <0.001*<0.001* |
| Fluorescence, volt | 0.07 ± 0.002 | 0.07 ± 0.002 | Selph, et al., (2021) | 2.921 | 0.995 | |
| Cycle | Gear | Reference | t test | P | ||
|---|---|---|---|---|---|---|
| In situ bulk zooplankton sizes (mm) | C1 | C5 | ||||
| μg C m−3 ± SD | ||||||
| 0.05–0.2 | 643.77 ± 156.0 | 421.85 ± 137.26 | Bongo-20 (0.05 μm mesh), 0–25 m depth | In situ | 1.64 | 0.18 |
| 0.2–0.5 | 377.8 ± 116.8 | 763.5 ± 443.5 | Ring net (0.2 mm mesh), 0–100/135 m depth | Landry and Swalethorp, 2021 | −1.698 | 0.150 |
| 0.5–1.0 | 637.1 ± 363.4 | 729.3 ± 623.9 | −0.226 | 0.830 | ||
| In situ preferred taxa > 0.2 mm | μg C m−3 ± SD | |||||
| Copepod nauplii | 0.3 ± 0.3 | 0.8 ± 0.6 | −1.326 | 0.242 | ||
| CladoceraCalanoid | 1.3 ± 1.9127.7 ± 30.4 | 48.5 ± 19.5376.9 ± 92.8 | Bongo-20 (0.2 mm mesh), 0–25 m depth | Shiroza et al., (2021) | −4.083−4.386 | 0.01*0.007* |
| Appendicularia | 0.6 ± 0.4 | 0.9 ± 0.5 | −0.696 | 0.518 | ||
| Abiotic variables | Mean ± SD | |||||
| Temperature, °CSalinity, psu | 24.49 ± 0.1336.40 ± 0.01 | 25.23 ± 0.4435.90 ± 0.27 | CTD, 0-25 m | Gerard et al., (this issue) | 1**180** | <0.001*<0.001* |
| Fluorescence, volt | 0.07 ± 0.002 | 0.07 ± 0.002 | Selph, et al., (2021) | 2.921 | 0.995 | |
Ring net sampled 0–100 m during C1 and from 0 to 133 m during C5. Bongo-20 sampled from 0 to 25 m. Of the taxa-specific groups, only those that were positively selected by ABT larvae were considered (Shiroza et al., 2021, Fig. 10). Please note that most copepod nauplii were small and not efficiently captured by the 0.2 mm mesh nets. Bottom panel indicates mean values for temperature, salinity and fluorescence for the upper 25 m determined from CTD casts. Fluorescence (volt) values are uncalibrated but close to mean values determined from extracted samples (mg m−3; Selph et al., 2021). Asterisk (*) indicates significant results at the 0.05 level for t-test; double asterisk (**) indicates Wilcoxon test.
DISCUSSION
This study provides new information on the growth of ABT larvae in their GoM spawning grounds. By comparing two sites that differed in temperature and food availability we identified the conditions most important for regulating larval growth. Furthermore, we provide validation for a newly proposed ABT larval FLI. Below we discuss our findings further in relation to environmental conditions and compare with other published studies.
Least square regressions for aged ABT larvae examined for (a) body size (SL, mm) at age (dph); (b) otolith radius (μm) at age for C1 (○) and C5 (●). Linear regression results are shown for C1 (˙˙˙˙) and C5 (˙˙˙) for length at age yC1 = 0.39x + 1.43, R2 = 0.78; yC5 = 0.37x + 1.74, R2 = 0.86 and exponential regression results OR at age yC1 = 7.19e0.14x, R2 = 0.83, yC5 = 6.52e0.15x, R2 = 0.91, respectively.
(a) Mean increment width trajectories (μm) at increment number (dph) for ABT from cycles C1 (○) and C5 (●). Error bars indicate standard error and plotted when at least five larvae were included. The dashed vertical (—) line placed between 7 and 8 dph indicates a significant observed change in IW between C1 and C5 larvae. This time roughly corresponds with the onset of flexion.
Least square regressions for aged ABT larvae from cycles C1 (○) and C5 (●) examined for standardized residuals of recent otolith growth ~ age (y-axis) and standardized (a) small prey biomass (b) FLI (c) residual of ingested preferred prey C ~ age (d) temperature °C from 0 to 25 m.
ABT larval growth
All biometrics changed significantly as the larvae grew older, with larval length, body weight and OR all continuously increasing during the first 17 days of larval life. Body depth, also referred to as muscular height and, is reported here for the first time, increased exponentially during early life. These findings are consistent with other studies on the ontogenetic development of larval bluefin tuna (García et al., 2017; Hernández et al., 2021; Malca et al., 2017). Small preflexion stage larvae were not as well represented at C1 compared to C5 and could have affected the slope at the base end of the growth curves. Nonetheless, biometrics were not significantly different between the two nursery areas, indicating that surviving larvae from cycles C1 and C5 grew similarly on average for the size range examined. We do note, however, that growth curves for the two cycles intersect between 8 and 14 dph, when larvae are in the flexion stage. Since larvae undergo substantial transformation in morphology, foraging capabilities and diet during flexion (Morote et al., 2008; Shiroza et al., 2021), this could indicate that conditions essential for preflexion and postflexion larval growth differed between the two nursery areas.
Combining larvae from both cycles, our average growth rate of 0.37 mm SL d−1 is somewhat lower than previously reported estimates from the GoM, although the ages reported here include individual dph corrections calculated for the first time. The dph corrections account for observed variability among larvae associated with estimating the missing increments between 7 μm from the otolith core to the first observed increment. On average, 2.29 ± 0.94 day were added to increment counts. Malca et al. (2017) reported larval growth rates of 0.46 mm SL d−1 in larvae collected throughout the northern and eastern GoM in 2000–2012 based on raw increment counts alone. Despite this disparity, average growth rates (SL ~ dph) reported in this study are 18% lower than those in 2012 (0.55 vs 0.67 mm d−1). The average growth rates reported in Malca et al. (2017) were corrected for dph by adding 2 days to increment counts. Compared to other ABT larval studies, our rates are slightly lower, but similar to those reported for other tuna species (0.3–0.51 mm d−1; Jenkins and Davis, 1990; Lang et al., 1994; Tanaka et al., 2006; Zygas et al., 2015; Gleiber et al., 2020a, 2020b). The slower growth rates of GoM ABT larvae in 2017–2018 reflect a spatially restricted sampling effort spanning 3 days during daytime-collections. These short cycles likely included small subsets of all cohorts from the 2017 and 2018 spawning seasons. In Malca et al. (2017), larvae were aged from multiple water masses that included over 100 larvae aged from samples collected at temperatures >27°C. In contrast, the mixed-layer temperatures for C1 and C5 ranged from 24.1 to 25.9°C.
Residual analysis of recent otolith IW-at-age, body depth-at-age and length-at-age revealed a high correlation among all three growth metrics. While it is generally assumed that otolith growth tracks somatic growth, this may not be the case for some species during certain developmental periods or under specific environmental conditions (e.g. Morales-Nin, 2000; Swalethorp et al., 2016). Nevertheless, our results support the use of IW and recent IW for assessing the full history and recent growth differences in wild-caught ABT larvae. IW increased exponentially during the first 14 dph, but growth curves intersected between 7 and 8 dph for the two cycles, suggesting more favorable growth conditions for preflexion larvae during C1. Conversely, during the flexion and postflexion stage conditions were more favorable in the C5 nursery area. IWs were comparable to those reported by Malca et al. (2017) for the GoM in 2012 only during the first 2 days post hatch, after which they were much narrower. However, IWs were comparable to Balearic Sea ABT larvae during the first 11 days and then widened (García et al., 2013). Ontogenetic changes impact growth (Hare and Cowen, 1995), and should be considered vital during otolith microstructure studies. It is possible that maternal investment could have differed between C1 in 2017, C5 in 2018 and other studies. Maternal investment has the capacity to affect larval survival and growth later in life (Berkeley et al., 2004; Masuma, 2009). Although a maternal effect has been proposed for growth of ABT larvae (Laiz-Carrión et al., 2019; Uriarte et al., 2016), it has not yet been established (Medina, 2020).
Environmental drivers of growth
In addition to potential maternal effects, the availability of zooplankton prey is generally considered important for sustaining larval growth (Houde, 1987; Houde, 2008). For ABT larvae of a given size or developmental stage, prey need to be of a suitable size, catchability and nutritional quality (Cushing, 1990; Robert et al., 2013). Within the thermal range that ABT larvae are adapted, higher temperatures can support faster development, if food is sufficient to support the increased metabolic and growth demand. We tested the importance of these and other measured variables using data and data products from C1 and C5 (Table II). Our GAM approach revealed that greater food limitation negatively impacted recent larval IW. The multivariate FLI, which is indicative of sufficiency of zooplankton biomass of suitable larvae-specific size to support ingestion rates that satisfy metabolic requirements, (Shropshire et al., 2021) was determined for each individual larvae based on its size and location. Not surprisingly, it was among the best explanatory variables for larval growth, which supports its use in assessing ABT larval habitat quality (Shropshire et al., 2021).
Interestingly, however, we did not find a direct positive relationship between our estimates of prey biomass availability and larval growth. SPB, which only consider the fraction of 0.002–0.2 mm zooplankton that falls within the prey size spectra of larvae (Shropshire et al., 2021), was negatively correlated with recent IW. This is surprising since much of the prey ingested by aged larvae in this study falls within this size range. The LPB, which considers the 0.2–1.0 mm fraction within the larval prey size spectra, fell out of the GAM model, but only because it was strongly and negatively correlated to FLI. Considering this, availability of larger-sized prey, such as other larval fish, may be important in regulating growth. Some gut content studies either exclude size classes that may be piscivorous or remove larval fish from diet contributions due to degradation and uncertainty of body length needed in carbon-length conversion factors (Gleiber et al., 2020b). This practice can avoid the reality that fast-growing larvae such as ABT begin targeting other fish as soon as opportunity and larval development align. Our four piscivorous ABT larvae included in this study show that prey >1 mm are targeted. Although these precocious larvae had larval fish in their guts, these were likely first time captures as the corresponding IW were within range of other larvae of the same age and size. Piscivorous larvae could have had larger recent IW had the prey in their guts been fully digested at time of capture.
It can be difficult to detect significant relationships between fish larval growth and bulk estimates of zooplankton available as prey (Robert et al., 2009; Swalethorp et al., 2016). One of the main issues is that fish larvae are nonrandom predators, feeding selectively on specific prey taxa and sizes available to them (Robert et al., 2009; Robert et al., 2013; Swalethorp et al., 2016; Shiroza et al., 2021). We did not have robust in situ estimates of preferred prey taxa availability covering the entire prey size spectra of the larvae to include in the GAM. Younger smaller larvae feed extensively on prey not quantitatively sampled by a 0.2 mm mesh Bongo net, and 0.05–0.2 mm samples for taxonomic analysis were only collected on C5. However, we did have diet information on feeding success expressed here as the amount of prey carbon each larvae had ingested. Using this information, ingestion of Calanoid nauplii and Podonidae, the two most preferred prey taxa by developing larvae (Shiroza et al., 2021), was positively correlated to recent IW. Interestingly, only including ingested preferred taxa explained growth variability substantially better than when including all ingested prey. This is a significant finding as it highlights the importance of considering prey taxa and sizes that are positively selected by larval fish, not only bulk zooplankton, when assessing nursery habitat quality and its potential for supporting growth and survival.
Temperature was the weakest explanatory variable for larval growth, likely reflecting the narrow range of habitat temperatures sampled in this study, all well within the optimal range for ABT larvae (Domingues et al., 2016; Muhling et al., 2010). Even though temperature has been found to be a principal driver of larval growth (Satoh et al., 2013), if prey availability limits feeding success (i.e. ingested prey), it becomes more important (Swalethorp et al., 2016).
The intersection of all growth curves between the two nursery areas coincided with the important transition of larval flexion suggested that C1 had better conditions to support preflexion larval growth, while C5 had better conditions for flexion and postflexion growth. As previously mentioned, we did not have estimates of preferred prey taxa availability in the 0.05–0.2 mm fraction from C1. However, small zooplankton biomass was highest at C1 within this size range (Table IV), suggesting that more of the preferred Copepod nauplii and small Calanoid copepodites could have been available there. Preflexion larvae were also found to feed significantly on ciliates (Shiroza et al., 2021) and could have been feeding on other microplankton which are hard to detect in stomach contents, but can provide important feeding opportunities for first-feeding fish (Overton et al., 2010; Scura and Jerde, 1977). C1 had the highest concentrations of ciliates; and highest concentrations and production of large autotrophs, including dinoflagellates (Landry et al., 2021; Selph et al., 2021; Shiroza et al., 2021). Considering the more oligotrophic conditions at C1 compared to C5, which had a stronger connection to the nutrient and plankton rich northeastern shelf (Gerard et al., this issue), it is expected that plankton communities would be shifted toward the smaller median sizes preferred by first feeding larvae. However, C1 prey communities did not match the larval ontogenetic shift in prey preferences.
For flexion and postflexion larvae, C5 was significantly richer in the preferred larger Podonidae and Calanoids. This richness explains why larvae at this developmental stage grew significantly faster at C5. It also further emphasizes the importance of advection of shelf zooplankton communities into ABT larval nurseries (Kelly et al., 2021; Landry and Swalethorp, 2021), to support faster larval growth and likely survival, since faster growing individuals are more adept at avoiding predators and less likely to starve (Tanaka et al., 2006; Watai et al., 2017).
CONCLUSION
We investigated the growth of larval ABT in the GoM and its environmental drivers by contrasting different developmental stages in two nursery areas that differed in prey availability. All growth biometrics examined increased as the larvae grew older, but their trajectories differed between areas. The principal drivers of this growth were related to food limitation/availability and feeding success. As a significant finding, the ingestion of preferred prey (Copepoda nauplii and Podonidae) better explained growth variability than total ingestion. While mean growth rates were similar, one nursery habitat appeared better suited to faster growth of preflexion stage larvae, while the other had faster growth of flexion and postflexion larvae due to greater availability of the preferred prey sizes and taxa of more developed larvae. Our findings underline the importance in considering ontogenetic differences in preferred prey rather than bulk zooplankton biomass when assessing habitat quality, and that growth limitation can occur at different larval ages among nursery areas that vary in both the quantity and composition of food resources.
DATA ARCHIVING
Data presented here have been submitted to the NOAA’s National Centers for Environmental Information data repository and will also be archived at Biological and Chemical Oceanography Data (BCO-DMO) Management Office site https://www.bco-dmo.org/program/819631.
ACKNOWLEDGEMENTS
We thank the crew of the National Oceanic and Atmospheric Administration (NOAA) ship R/V Nancy Foster for their support in conducting the survey during both years. Unlike other systematically planned surveys, the navigating officers were required to adjust course and maneuver in situ during our search for the larval patches, and the deck crew were always standing by to support the deployment and the collection of the gear and samples. We also thank Sarah Privoznik for her assistance in organizing the cruises, Amelia Jugovich for assistance analyzing the ichthyoplankton samples, and Jason Mostowy and Barbara Muhling for providing valuable edits to an earlier draft of this manuscript. Finally, we acknowledge the ECOLATUN (CTM-2015-68473-R MINECO/FEDER) project for establishing the sampling protocol.
FUNDING
BLOOFINZ-GoM Program NOAA RESTORE Science Program (NOAA-NOS-NCCOS-2017-2004875 to J.T.L. and T.G.; NA15OAR4320 071 to M.R.L.; NA16NMF4320058 to M.R.S.); US National Science Foundation (OCE-1851558 to M.R.L.); Cooperative Institute for Marine and Atmospheric Studies (#NA20OAR4320472); and Nova Southeastern University Ph.D. Fellowship to E.M.
References
Fox, J. and Weisberg, S. (
