Environment-driven trends in larval abundance predict fishery recruitment in two saltwater basses

Abstract

Environmental and biological factors influencing fish larvae can drive fishery cohort strength, yet larval abundance is typically a better indicator of spawning biomass. Under a changing ocean, studies that explore the relationships between environmental variables, larval abundance, and fishery recruitment remain valuable areas for ongoing research. We focus on a popular, recreational-only, multispecies saltwater bass fishery (genus Paralabrax) whose population status and recovery potential are uncertain. We resolved Paralabrax spp. larval data to species over a 54-year period (1963–2016) and used species distribution models to (i) generate and test species-specific standardized indices of larval abundance as indicators of adult stock status and fishery recruitment and (ii) evaluate long-term spatiotemporal trends in their population dynamics relative to environmental variables and climate forcing. Contrary to initial hypotheses, species-specific larval abundance predicted future catches, with higher recent larval abundances suggesting potential for fishery recovery. Temperature, zooplankton biomass, and isothermal depth were important predictors of bass larval abundance, indicating these variables could also be valuable for predicting fishery recruitment and anticipating population change. Our findings paint a path forward towards a more ecosystem-based fishery management approach for this important fishery and may serve as a template for data- or assessment-limited fisheries.

Introduction

An ecosystem approach to fishery management (Dolan et al. 2016) considers the many biological and environmental factors that have influenced and continue to shape fished populations. This, of course, also requires having data and resources to adequately evaluate the status of fish populations. Unfortunately, many small-scale fisheries lack such a robust approach, partly because the data and resources for managing these fisheries are typically limited. Given that the economic and ecological impacts of some marine recreational fisheries have rivalled those of commercial fisheries (Cooke and Cowx 2004, Lewin et al. 2019), and that overexploitation dampens the ecosystem resilience of fish populations (Perry et al. 2010, Ziegler et al. 2023), there is a clear need for research that leads to more robust population assessments for recreational-only and small-scale commercial fisheries.

Southern California is uniquely data-rich in terms of long-term coastal fisheries oceanography data and thus represents an ideal testbed for the development of ecosystem-assessment methods. The California Cooperative Oceanic Fisheries Investigations (CalCOFI) is the world’s largest and longest-standing fisheries oceanography survey (McClatchie 2014, Gallo et al. 2019), and since 1951, it has conducted quarterly surveys to collect biological samples (e.g. fish larvae abundance, zooplankton biomass) and hydrographic data at discrete sampling locations throughout the Southern California Bight (SCB). Due to its long temporal coverage and high spatial resolution, the survey has enabled researchers to detect a secular warming trend in southern California (McClatchie et al. 2018), ecosystem responses to climate forcing (Hsieh et al. 2009, Asch 2015, Thompson et al. 2017, 2022), and has assisted resource managers in fishery assessments (Moser et al. 2001b, McClatchie 2014, Gallo et al. 2019). Combined with other oceanographic and fishery data streams, the CalCOFI fish larval density time series provides a unique opportunity to assess historical and recent spatiotemporal trends in important fished populations.

One of the most popular and economically important recreational-only marine fisheries in California, USA, and perhaps even the world, is the southern California saltwater bass Paralabrax spp. fishery, which consists of Barred Sand Bass P. nebulifer (BSB), Kelp Bass P. clathratus (KB), and Spotted Sand Bass P. maculatofasciatus (SSB). This multi-species fishery began in the early 20th century, with significant declines in BSB and KB catches occurring in the mid-2000s, following a period of rising fishing mortality and recruitment failure (Jarvis et al. 2014a). Catches have remained at historic lows since 2013, largely reflecting a depressed BSB population (Jarvis et al. 2010, 2014a, Erisman et al. 2011), as BSB are particularly vulnerable to harvest due to the predictability of their large summer spawning aggregations (Jarvis et al. 2010). In contrast, KB spawning aggregations are smaller and more widely distributed (Erisman and Allen 2006). After the fishery decline, regulations on saltwater bass fishing became stricter in 2013 (Calif Code of Regs 2023; Jarvis et al. 2014a). Despite these changes, the multi-species fishery has not recovered, adding further uncertainty to the status and recovery potential of BSB and KB populations, in particular (SSB are primarily catch-and-release and presumed to have a healthier population). Adding to the uncertainty is the unknown impact of a rapidly changing ocean on these populations.

The long-term CalCOFI Paralabrax larval time series offers an invaluable opportunity to unravel uncertainty concerning the BSB and KB populations. Successful species discrimination across Paralabrax larval stages was recently made possible with the development of a robust, species-specific larval taxonomic key (Jarvis Mason et al. 2022), enabling us to reconstruct the CalCOFI Paralabrax larval time series for the first time. Larval abundance has long been considered a proxy for spawning stock biomass (Hilborn and Walters 1992, Cowan and Shaw 2002). Indeed, CalCOFI larval abundance data have been incorporated in west coast fishery stock assessments or used as an index of spawning stock biomass to monitor status and trends of adult populations (Gallo et al. 2019). Although CalCOFI is a uniquely intensive fish larval survey, regular surveys for larval fishes, even with minimal field equipment and small boat operations, hold potential to capture fisheries time series. On one hand, if larval abundance directly reflects the biomass of reproductive females, fishery managers can easily track trends in adult stocks. On the other hand, larval abundance may also predict fishery recruitment in some species (Murphy et al. 2018, Tripp et al. 2023), as environmental factors heavily influence larval dynamics, with growth and mortality processes shaping future fishery year class strength (Hjort 1914, Lasker 1981, 1987, Cushing 1990, Houde 2001). However, due to continued abiotic and biotic pressures on subsequent life stages and differential ontogenetic responses to environmental conditions, the juvenile or subadult stage, rather than the larval stage, is generally considered to have higher predictive ability for forecasting fishery recruitment (Bradford 1992, Cowan and Shaw 2002). Thus, by resolving the CalCOFI Paralabrax spp. larval data to species, we can explore long-term, species-specific bass larval trends and how they relate to population status and fishery recruitment.

BSB and KB grow at similar rates and become susceptible to hook-and-line fishing gear at 2 to 3 years of age, reaching size at fishery recruitment between 5 and 7 years for several decades (1959–2012, Love et al. 1996, Jarvis et al. 2014a). Following 2012, an increase in the minimum size limit from 12 to 14 inches corresponded with an average fishery recruitment age of 8 years (Walker et al. 2020). If bass larval abundance can predict fishery recruitment, rather than (or in addition to) being a proxy for female spawner biomass, then a cross-correlation analysis should provide evidence of larval abundance trends leading fishery catch trends by one or more yearly lags. These lags would correspond with fishery recruitment age and older, as some individuals of a strong cohort may eventually be caught at an older age. This would be in contrast to an expected lag of zero, if larval abundance is a proxy for female spawner biomass (and if saltwater bass catches are representative of the spawning stock).

Reconstructed larval time-series for BSB and KB would also provide an opportunity to explore each of their responses to changing environmental conditions, which may help with anticipating future change and improve management to aid in fishery recovery. For example, differences in the geographic ranges and reproductive strategies across saltwater bass species suggest distinct sensitivities to fishing pressure and environmental conditions, and thus, differences in their resilience to a changing ocean. The saltwater basses are typically found in southern California, south of Point Conception (Love and Passarelli 2020). However, KB have been known to occur as far north as the Columbia River of Washington, USA, while BSB have been recorded north only to Monterey Bay, CA, making the latter a less common vagrant to more temperate systems. The southernmost occurrence for KB is southern Baja California, Mexico, and for BSB, it is the tropical waters of Acapulco, Mexico (Heemstra 1995). As a group, abundance of saltwater basses appears to be positively correlated with warmer ocean temperature regimes (Moser et al. 2001b, Hsieh et al. 2005, Jarvis et al. 2014a), and at one time, the basses were common north of Pt. Conception (Hubbs 1948). Thus, it is possible that climate-driven increases in sea surface temperature (SST) might influence future distribution shifts or expansions (Auth et al. 2018). At the same time, events on shorter temporal scales, such as seasonal upwelling and El Niño Southern Oscillation (ENSO), may also influence population dynamics and distributions. Therefore, understanding environmental influence over saltwater bass populations, within the context of their vulnerability to harvest impacts, should help elucidate how species may respond to oceanographic climate change and which species may be more resilient.

Here, we unlock the CalCOFI Paralabrax spp. larval archive for the first time to

• generate and test standardized indices of saltwater bass larvae abundance as indicators of adult stock status and fishery recruitment, and

• evaluate spatiotemporal trends in bass larval abundance relative to environmental variables and climate forcing (e.g. environment–species relationships and latitudinal shifts in abundance).

In the face of a rapidly changing ocean, the ability to identify potential environmental indicators of population status or future fishery recruitment, i.e. potential management action triggers, will be ever more important for guiding sustainable fishery management and may be key in determining the fate of this popular fishery.

Methods

Survey samples and species identifications

BSB and KB are batch spawners, capable of daily spawning (Oda et al. 1993, Jarvis et al. 2014b), with a summer spawning peak in July (Walker et al. 1987); their larvae occur primarily in nearshore habitats off southern California and Baja California (McGowen 1993, Moser et al. 2001a). We obtained Paralabrax spp. larvae collected during July CalCOFI oceanographic cruises from 1963 to 2016 along the continental shelf off southern California, USA, and Baja California, Mexico (Fig. 1). Two changes to the CalCOFI zooplankton sampling occurred during this period: in 1969, the depth from which nets were obliquely towed increased from 140 to 210 m (or from within 5 m of the seafloor in shallower waters), and in 1978, the sampling gear changed from ring nets to paired bongo nets (Thompson et al. 2017). The period chosen for our study consisted of the most up-to-date taxonomically reviewed records available at the time of the study, i.e. we assume no misreporting.

July survey effort varied spatially and temporally, in which the most consistent coverage occurred in the SCB (Fig. S1); with the exception of 1983 (only two stations surveyed), each cruise in southern California consisted of at least 6 transects. July CalCOFI cruises off Baja California were spatiotemporally consistent from 1963 to 1978, but afterward, spatial and temporal coverage was limited (Fig. S1).

Larval species identifications and data processing

Despite limited July survey efforts off Baja California after 1978 (Fig. S1), we nevertheless chose to taxonomically identify all samples from both regions (Fig. 1) to provide a characterization of the spatial distribution and relative abundance of KB and BSB larvae collected in both regions overall and during the period of consistent spatial coverage. Subsequently, we focused the remainder of our analyses on southern California only.

Two taxonomists independently determined the larval stage of each larva based on notochord development and then assigned species using pigmentation patterns described in Jarvis Mason et al. (2022). We excluded yolk sac larvae as these could not be reliably identified. After performing quality assurance procedures and data processing (Supplementary Material S1), we adjusted the Paralabrax spp. counts in each sample by the percent of sample sorted and the standard haul factor to obtain the number of larvae ‘under 10 m² of sea surface’ (Thompson et al. 2017). We refer to this adjusted CalCOFI raw count as larval density.

General modelling approach

We focused our analysis on BSB and KB in southern California, the region with the most consistent sampling coverage. We built species-specific species distribution models (SDMs) of July larval density with two aims: (i) develop and test standardized indices of bass larval abundance as proxies for spawning stock biomass and fishery recruitment and (ii) evaluate the effect of environmental factors on bass larval densities. For our first aim, we built species-specific index of abundance models that accounted for variables that may affect ‘catchability’ but not abundance (e.g. survey design). By not accounting for environmental variability in the index models, the standardized indices of abundance could be more fairly tested against catch data (our proxy for female spawner biomass), which also contain environmental variability. The larval abundance indices were also used to explore latitudinal shifts in abundance (see the ‘Latitudinal shifts in larval abundance’ section). For our second aim, we constructed another set of models that related species-specific larval density as a function of environmental covariates while also accounting for factors influencing catchability (see the ‘Environmental influence on bass larval densities’ section).

Standardized indices of abundance

We used a geostatistical generalized mixed effects model framework that incorporated spatial and spatiotemporal random fields (effects) to account for latent variables that may contribute to observations being correlated in space and time. The spatial random field captures latent spatial variation that is unchanging through time (e.g. depth), while the spatiotemporal random field captures latent time-varying spatial patterns (e.g. dynamic biological and oceanographic processes not captured by covariates in the model). This type of SDM can therefore account for unbalanced sampling effort (e.g. a few missing years/stations) and is more precise than other SDM methods at estimating abundance (Thorson et al. 2015, Brodie et al. 2020). We fit SDMs with the R package sdmTMB (Anderson et al. 2022), which estimates parameters of the Gaussian random fields using the Integrated Nested Laplace Approximation and implements maximum marginal likelihood estimation with Template Model Builder (Anderson et al. 2022).

We modelled bass larval density Y at location s and time t using a Tweedie observation error family [positive continuous density values that also contain zeros, Shono (2008)] and a log link:

μ_s,t represents mean density at location s and time t,

X_s,t is a vector of main effects, time-varying effects, or spatially varying effects corresponding to location s and time t,

β is a vector of estimated parameters, and

ω_s and ε_s,t are the spatial and spatiotemporal random effects, respectively, assumed to be drawn from a multivariate normal distribution with covariance matrices |$\sum_{\rm{\omega }_s}$| and |${\sum _{{\varepsilon _{s,t}}}}$|⁠. A user-specified triangulated mesh approximated the spatial component of the Gaussian random fields (Fig. S2, Supplementary Material S2, Anderson et al. 2022).

For both species, candidate index models included various combinations of year, moon phase, day of year (Julian day), hour of day, net type, and distance to the mainland coast. We did not include a factor for the increase in tow depth because southern California Paralabrax spp. larvae are most abundant inshore of the 36 m depth contour (Lavenberg et al. 1986). We modelled year effects and all other covariate effects with penalized smooth functions, i.e., penalized regression ‘P-splines’ (Eilers and Marx 1996), similar to GAMs (generalized additive models; Wood 2017). To account for potential daytime net avoidance related to the gear change from ring to bongo nets (Thompson et al. 2017), we specified a day/night by net type interaction where a separate smooth was created for each combination of net type and day/night. We allowed the model to select the basis dimension, i.e. k, or ‘wiggliness’, for each variable, except for the hour and continuous numeric moon phase covariates, in which we specified a cyclic smooth over 24-h and 4-w cycles, respectively. In the case of BSB, in which the inclusion of a smooth year effect was not supported, we incorporated year effects as a time-varying random walk (analogous to a dynamic linear model), which allowed for each year to function as the intercept. For prediction of larval abundance in years without data, we specified missing years for model estimation of larval abundance. In cases of non-convergence, we tried simplifying the catchability covariates to their linear form, i.e. no smooth function.

For both species, we identified the best-fit random effects model structure to include both an independent spatiotemporal random field across years and a spatial field (we also tested models with the spatial field off). We used Akaike (AIC) values (model specified with Restricted Maximum Likelihood [REML] = off) to evaluate model parsimony across candidate models with different main effects. We report results only for models that passed all model diagnostics; this was determined using the function ‘sanity’ in sdmTMB, which checks for convergence, a positive definite Hessian matrix, large standard errors associated with model parameter estimates, and more (Anderson et al. 2022). On this smaller subset of candidate models, we conducted 8-fold cross-validation, in which we used eight folds containing random data subsets in space and time to train the model. For example, in each iteration, one data fold was excluded and the remaining data was used as the training set. The model’s predictive performance was then evaluated on the excluded fold. We used the sum of the negative log likelihood across all folds for each candidate model to determine the model with the highest predictive performance for each species (maximizing the likelihood for held-out data).

Trends in larval abundance

To derive an area-weighted standardized index of July larval abundance for southern California, we used our index model to predict abundance across a 2 km × 2 km fine-scale grid covering the survey region for all years (Anderson et al. 2022). For KB, in which day was a predictor, we predicted on the data-scaled Julian value for July 15th. We used the ‘get_index’ function within sdmTMB to derive relative larval abundance per year based on the summed spatial predictions for each year and plotted these temporal trends for each species in normal and natural log space. To explore changes in the July distribution of larvae through time, we plotted spatiotemporal trends in larval abundance of both species, including both fixed and random effects. We also separately plotted trends in the spatiotemporal random effects (⁠|${\varepsilon _{s,t}}$||$)\,\,$|to visually explore patterns driven by latent (unaccounted for) spatiotemporal processes.

Latitudinal shifts in larval abundance

We used larval abundance predictions from the standardized index of abundance models and the function ‘get.cog’ within sdmTMB to derive a mean predicted July centre of gravity (COG) across available years for both species, with the exception of 1983, which included only two stations (Thorson et al. 2016, Anderson et al. 2022). Given that temperatures within the CalCOFI region have shown a secular warming trend over much of the time series (Gallo et al. 2019), and given the oceanography of the system (cooler temperate conditions to the north), we predicted that any associated latitudinal shifts in bass larval distribution through time would be to the north. We visually and statistically evaluated the data for a northward progression of predicted COG estimates by mapping the COG by year, colour coding the predictions on a graduated scale from dark (early decades) to light (later decades), and testing for a difference in the median predicted COG northings (latitudinal coordinates in km) between the 1960s and 2010s using a Wilcoxon test. Additionally, we used a simple general linear model to test for a linear relationship between mean July SSTs and mean annual predicted COG northings. We obtained mean July SSTs from the corresponding CalCOFI cruise stations, with the exception of 1981 (no SST data), in which we used the mean July SST off Pt. Dume, CA, which is approximately centrally located along the southern California coast (Carter et al. 2022). Finally, we calculated 50, 75, and 90th kernel quantile densities for both species’ COGs using the ‘kde2d’ function within the R package MASS (Venables and Ripley 2002) to provide a comparison of the size and location of the two species’ overall COG in southern California.

Relationship of saltwater bass larval abundance to fishery catch data

We used Spearman cross-correlation to test for a positive monotonic correlation between predicted July larval abundance and annual adult catches and considered only unidirectional lags, from zero to ten years. We assumed that if bass larval abundance is a function of the biomass of reproductive age females in the water, then that relationship should be detected with the strongest positive cross-correlation at a lag of 0 between our standardized index of abundance and our proxy for spawning stock biomass, i.e. adult catches. If instead, bass larval abundance can better predict fishery recruitment, then we should see evidence of larval abundance leading catches, having the highest correlations occurring at lags corresponding with fishery recruitment age (5–8 years) and older, as some individuals of a strong cohort may eventually be caught at an older age. We standardized all data (predicted larval abundance and catch) to a mean of zero. We performed cross-correlation analysis using the R package funtimes (Lyubchich et al. 2023), in which a 95% confidence band provided guidance for interpreting correlations that may be influenced by the presence of autocorrelation in either or both datasets, i.e. coefficients falling within the band may be an artifact of non-stationarity.

We used catches to represent female spawner biomass given that BSB catch and adult densities are correlated (Jarvis Mason et al. 2024) and because KB catch directly corresponds to fishery catch-per-unit-effort (Jarvis et al. 2014a). Additionally, the catch data represented the only long-term index of saltwater bass abundance that matched the geographic extent of the CalCOFI larval abundance data. We used two catch data sets as proxies: (i) harvested catch data (=fish kept) from California Department of Fish and Wildlife Commercial Passenger Fishing Vessel logbooks (Fig. S3; these sportfishing vessels take paying customers fishing and are the predominant fishing mode for BSB and KB), and (ii) estimates of total catch across all fishing modes (Jarvis et al. 2014a). The total estimates may represent a more complete picture of adult abundance because, in addition to kept fish, it also includes catch-and-release fishing, i.e. released fish, including fish below the minimum size limit (‘shorts’) or smaller adults. We obtained the total catch estimates from the Marine Fisheries Statistical Survey (1980–2003) and the California Recreational Fisheries Survey (2005–2016; RecFIN 2023); see Supplementary Material S3 for a description of the surveys and associated data processing for cross-correlation analysis. For both species, the bulk of the catch occurs during the summer spawning period (June through August, CDFW 2021, 2023).

Environmental influence on saltwater bass larval densities

We used the R package sdmTMB to model the influence of site-specific prey availability (e.g. zooplankton biomass), SST, and larger-scale oceanographic indices on BSB and KB CalCOFI July larval densities. As with the standardized index of abundance models, we tested covariate effects as smooth functions, and in cases of non-convergence, we tried simplifying the environmental covariates to their linear form, i.e. no smooth function. Not all covariates were available from 1963, so we first modelled a smaller set of variables across the entire time series (1963 and 2016, Time Period 1) and then incorporated additional variables over the shorter period (1984–2016, Time Period 2). This also provided the opportunity to visually inspect for non-stationary environment–species relationships between the two time periods, i.e. are the magnitude and direction of the covariate effects consistent? For these environmental models, we identified the best fit random effects model structure to include an independent spatiotemporal random field across years and no spatial field.

In Time Period 1, the entire time series, we included station-specific CalCOFI SST (°C, averaged across the upper 10 m) and square root transformed zooplankton biomass measured as zooplankton displacement volume (cm³/1000 m³ strained). This metric can be a crude index of zooplankton community structure as it is highly influenced by pelagic tunicates (e.g. salps) and other large gelatinous organisms (Lavaniegos and Ohman 2007). To better elucidate the potential impact of pelagic invertebrate prey, we also explored the relationship between bass larval density and biomass of the calanoid copepod, Calanus pacificus, from spring cruises (data was not available for July cruises). We reasoned that increased copepod biomass available to adult bass prey species (e.g. coastal forage fishes) in the spring might translate into positive maternal conditioning of larvae in summer (Walsh et al. 2024).

In Time Period 1, other environmental covariates included the North Pacific Gyre Oscillation (NPGO) index for the month of July, the Oceanic Niño Index (ONI) for June/July, and mean Bakun upwelling indices (m³/s) for June/July off southern California (latitude 33°N) and northern Baja California (latitude 30°N; Fig. S4). The NPGO index is based on sea surface height variability throughout the North Pacific and tends to track regional nutrient fluctuations; mechanisms driving plankton ecosystem dynamics (Di Lorenzo et al. 2008). Upwelling was included as a seasonal measure of potential nutrient flux (versus the larger-scale NPGO), where positive values indicate wind-driven offshore transport (Bakun 1973). We considered upwelling off Baja California because older BSB recruits (∼ age 2) showed a negative correlation with upwelling south of the SCB at a 2-year lag (Jarvis et al. 2004). The ONI is based on surface temperature anomalies in the east-central equatorial Pacific and defines ENSO states (NOAA 2023).

In Time Period 2, the latter part of the time series, we compared models that included all covariates from Time Period 1 and the isothermal layer depth (ILD, Fig. S4), which we included because spawning BSB orient to the thermocline during spawning season (McKinzie et al. 2014). In the study region, the ILD roughly corresponds to the mixed layer depth (MLD), which is the depth where temperature differs by 0.5°C from the sea surface. Additional upwelling indices we considered in Time Period 2 were the BEUTI (Biologically Enhanced Upwelling Index) and CUTI (Coastal Uptake Transport Index) specific to the southern California region (Fig. S4). The BEUTI and CUTI provide a measure of nitrate flux through the mixed layer and coastal vertical transport, respectively (Jacox et al. 2018a, data from Hunsicker et al. 2022). Positive BEUTI values indicate nitrate flux into the mixed layer, and positive CUTI values indicate coastal upwelling. For KB, we also considered kelp canopy extent (Fig. S4), because giant kelp Macrocystis pyrifera is important adult habitat (White and Caselle 2008). We obtained Bight-wide areal kelp canopy data by year from kelpwatch.org.

All species-specific environmental models also accounted for the catchability factors we identified in the index of abundance models (e.g. distance to mainland coast, day of year), except for year effects. As with other mixed models, our spatiotemporal model in sdmTMB did not permit missing values, so prior to analysis, we used the R package gstat (Pebesma 2004) to fill in within-cruise CalCOFI temperature data with the inverse distance-weighted interpolation method. On average, within each cruise, there were 11 station temperatures interpolated from 48 stations. There were too few stations surveyed in 1983 to interpolate data, and there was no temperature data available for 1981 (Fig. S4). We checked for multicollinearity among model covariates using the R packages corrplot (Wei and Simko 2021) and performance (Lüdecke et al. 2021). Most of the 66 correlations among environmental covariates were weak (<0.20, Fig. S5), and all model checks for multicollinearity produced ‘Low Correlation’ results, with most covariate variance inflation factor values close to 1 and none greater than 4. We selected models and evaluated model performance in the same way as stated for the index models. We visually explored the conditional effects of important explanatory variables using the R package visreg (Breheny and Burchett 2019).

Results

We processed a total of 1 716 Paralabrax spp. larvae sorted from 143 original CalCOFI station vials from July surveys off southern California and Baja California between 1963 and 2016 (Figs S6 and S7; Table S1); however, the results of our analysis below are focused on southern California, the region with consistent survey coverage over the entire study period. For a characterization of the geographic extent and relative abundance of CalCOFI Paralabrax spp. larvae identified by species and region, see Supplementary Material S4.

Trends in saltwater bass larval abundance

The most parsimonious index standardization models with the highest performance for BSB included time-varying year effects and distance to mainland coast, and for KB included year and day of year (Table 1). The 54-year period showed four strong sporadic larval pulses in 1981, 1989, 2012, and 2014 against relatively low background levels; the most recent larval peaks in the mid-2010s were near all-time highs (Fig. 2a and c). Overall, predicted July mean natural log abundance (±SD) was higher and less variable for KB (10.90 ± 0.856) than BSB (9.01 ± 1.28; ANOVA: F = 44.82, P < 9.5e-9, Fig. S8). Back transformed, mean July KB larval abundance was 6.6x higher than BSB. In terms of model coherence, the predicted larval abundance trends aligned with expectations, following with both empirical trends observed in this study (adult catches; Fig. S3) and young-of-the-year trends in the literature (Miller and Erisman 2014, Jarvis Mason et al. 2024). Both species showed dramatic increases in larval abundance between the 1960s–70s and the 1980s (Fig. 2b and d). BSB larval abundance steadily declined from the mid-1990s through the early 2000s before increasing again to higher levels through the mid-2010s (Fig. 2b). In contrast, a steady decline in KB larval abundance began earlier, from the 1980s, through the early 2000s, before increasing again into the 2010s (Fig. 2d).

Predicted July spatial distributions also matched species-specific expectations, with BSB larval abundance being higher nearshore and along the central and southern coast (Fig. 3), while KB larval abundance was relatively higher throughout the SCB, with hotspots occurring primarily at the northern Channel Islands (Fig. 4, see Fig. S2 for location reference). Spatiotemporal random effects plots for BSB and KB also showed different patterns driven by latent forces through time, in which KB showed more spatiotemporal structuring than BSB (Figs. S9 and S10).

Latitudinal shifts in bass larval abundance

The July BSB larval COG in southern California was smaller and located more nearshore and central along the coast than the KB larval COG (Fig. 5). Over the study period, the estimated July COG for both species showed no apparent northward trajectory (Fig. 5a and c), and there was no difference between the median predicted COG northing in the 1960s and 2010s for BSB (Wilcoxon: W = 10, P = .69, 1960s: 3725 km, 2010s: 3726 km) or KB (Wilcoxon: W = 11, P = .84, 1960s: 3714 km, 2010s: 3712 km). However, we found that mean July SSTs explained 26% of the variability in July COG northings for KB, where higher SSTs generally resulted in a more northward larval extent (glm: β = 12.061 ± 3.878 SE, R² = 0.26, P = .004; (Fig. 5b and d).

Relationship of saltwater bass larval abundance to fishery catch data

We found no relationship at a lag of zero between July bass larval abundance estimates and annual fishery catch data (our proxy for the relationship between larvae and spawners; Fig. 6, see Fig. S11 for the index estimated for all years). In contrast, for both species, we found moderate to strong lagged relationships between larval abundance and both sources of catch data (reported Commercial Passenger Fishing Vessel harvest and total catch estimates across all fishing modes). For example, the highest correlation coefficients between BSB larval abundance and total catch estimates occurred at 6- and 7-year lags, corresponding to the age of most fishery recruits (Fig. 6a). We found a similar trend with the BSB harvest data, in which the highest correlations indicated larval abundance led catch by 7 to 10 years. For KB, the species with greater catch-and-release fishing, there were moderate to strong correlations between larval abundance and total catch estimates from 3- to 10-year lags; with harvest, KB larval abundance showed weak, marginally significant correlations at 4- to 6-year lags and strong correlations at fishery recruitment age and beyond (7- to 10-year lags, Fig. 6b).

Environmental influence on saltwater bass larval densities

Given that we standardized predictor variables, the relative effect of each was directly comparable. Over the entire study period, 1963–2016 (Time Period 1), CalCOFI station-specific zooplankton biomass and SST were important positive predictors of July bass larval density for both species, but the effect of zooplankton was estimated to be smaller (Fig. 7). Although the distribution of zooplankton biomass was positively skewed (mostly low biomass values), bass larval density was estimated to be higher at higher zooplankton biomass extremes (Figs. S12 and S13). Inclusion of C. pacificus did not improve model performance (Table 2). For both species, there was also a strong negative relationship with distance to mainland coast (higher larval density nearshore; Figs. 7, S12 and S13). The effect of temperature and distance to mainland coast was stronger for BSB than KB, while the positive effect of day of year was stronger for KB than for BSB. The NPGO, ONI, and Bakun upwelling indices were not important covariates for either species’ larval abundance (Fig. 7).

Over the latter period, 1984–2016 (Time Period 2), the ILD was strongly negatively related to bass larval density (higher bass larval density with shallower ILD; Figs. 7, S14, and S15). Among the different upwelling indices, none had a strong effect on bass larval density (Table 2). Bight-wide kelp canopy extent did not improve model fit for KB (Table 2).

There was some consistency in the estimated covariate effects across the two time periods examined (both in magnitude and precision of the effect). The positive relationship with zooplankton biomass and negative relationship with distance to mainland coast remained important for both species during Time Period 2. For BSB, SST also remained important, while the effect of the NPGO changed from no effect to a positive effect, and the previously weak relationship with day of year changed to no effect (Fig. 7).

Discussion

Our reconstructed larval indices of abundance span over half a century, providing the longest species-specific, fishery-independent time series for the saltwater basses in the region. We demonstrate that bass larval abundance has value in predicting fishery recruitment, as both BSB and KB larval abundance had moderate to strong predictive power in forecasting fish catches (Fig. 6). This predictive ability was further characterized by strong links between larval density and environmental factors, paving the way for integrating larval abundance indices and environmental data into ecosystem-based fishery assessments.

Importantly, our study also highlights key differences in the population dynamics of BSB and KB, despite their shared management under the same sportfishing regulations. First, although both species exhibited sporadic larval pulses, mean KB larval recruitment was less variable and ∼7x higher compared to BSB. This finding is consistent with a southern California KB population that has been more persistent and more resilient to harvest impacts (Jarvis et al. 2014a). Second, relative to KB, we found BSB larval recruitment was more closely tied to temperature, which is consistent with a history of sporadic, warm water BSB juvenile recruitment events in southern California (Jarvis Mason et al. 2024).

Since environment-driven recruitment variability influences both a stock’s resilience to overfishing and its recovery potential (Kuparinen et al. 2014), ongoing warming and the rise in marine heatwaves (Oliver 2019) may lead to more stable and substantial local larval recruitment in the future (Jarvis Mason et al. 2024). Third, the offshore islands in the northern SCB are a notable hotspot for KB larval recruitment, which has implications associated with habitat protection (e.g. Marine Protected Areas). The Northern Channel Islands comprise more than a quarter of the area within southern California’s network of Marine Protected Areas, which, through protection of adult KB and other harvested species such as California Sheephead Bodianus pulcher, an herbivorous predator, may serve to enhance larval settlement and KB recruits through the maintenance of healthy kelp habitat (Hopf et al. 2022). Here, we discuss potential factors influencing the predictive power of bass larval abundance, implications for both species’ fishery recovery and future resilience relating to important species-specific differences in their population dynamics and environmental influence, and the importance of continued monitoring and assessment.

Saltwater bass larvae predict future catches

Historically, the challenge of predicting fishery recruitment has been complex (Haltuch et al. 2019), as predictive ability is influenced by factors affecting both pre- and post-recruitment stages, including maternal effects (Houde 2008). In fact, very early life stages like eggs and preflexion larvae are more often a suitable proxy for spawning stock biomass (Hilborn and Walters 1992, Gunderson 1993) and typically show weak relationships with subsequent juvenile and fishery recruitment (Cury et al. 2014, Szuwalski et al. 2015). Thus, it is notable we found that young bass larval abundance predicted fishery recruitment rather than spawning stock biomass (Fig. 6). Although rare, there are other instances where fishery recruitment is determined in the first year of life (Cushing 1990, Stige et al. 2013, White et al. 2019, Schilling et al. 2022), in which rapid initial larval growth and a brief preflexion larval stage may provide survival benefits into adulthood (Fontes et al. 2011, Robert et al. 2023). Assuming there is no major interannual variation in batch fecundity for the basses (Oda et al. 1993, Jarvis et al. 2014b), a strong relationship between bass larval abundance and spawning stock biomass in the same year would have to rely on (i) appreciable spawning occurring every year, (ii) spatiotemporally consistent mortality prior to larvae being surveyed, and (iii) appreciable local spawning and retention of larvae. In the case of the saltwater basses, the apparent breakdown in this relationship may be partly driven by variable larval mortality or spawning behaviour associated with environmental conditions and harvest impacts. For BSB, there also appears to be a lack of consistent, appreciable, locally sourced recruitment (Fig. 2; Jarvis Mason et al. 2024). Additionally, populations at their geographic margins tend to experience higher larval recruitment variability (Myers 1991, Neill et al. 1994, Levin et al. 1997), which is typically reflected in future fish catch (Armsworth 2002). Schilling et al. (2020), e.g. found a positive correlation between predicted Bluefish (Pomatomus saltatrix) larval ‘settlement’ (transformation to the juvenile pelagic phase) and fishery catch-per-unit-effort at the southern end of its distribution.

Regardless of the mechanism, our findings are consistent with other saltwater bass studies showing positive relationships between early life history indices of bass abundance (larvae, young-of-the-year juveniles) and future saltwater bass fishery recruitment, despite different collection methods and study time frames (Miller and Erisman 2014, Jarvis et al. 2014a, Jarvis Mason et al. 2024). Although the catch data in our analysis comprised multiple size classes, we were able to detect moderate to strong correlations at lags corresponding to the age of fishery recruits and older (Fig. 6). Future analyses using data focusing solely on fishery recruits may strengthen predictive ability (Jarvis et al. 2014a, Schilling et al. 2022). Interestingly, by including all size classes, we were also able to detect a species-specific difference in correlations that we attribute to the greater popularity of catch-and-release fishing for KB (e.g. correlations between KB larval abundance and total catch also occurred at lags corresponding to ages younger than size at fishery recruitment; Fig. 6a and c).

Using catch data as a proxy for spawning stock biomass assumes a proportional relationship between catch and abundance. For other aggregation-based fisheries like BSB, there is potential that hyperstability in catch data, i.e. consistent, high catches on spawning aggregations, may result in a nonlinear relationship between catch and abundance and an overestimate of abundance at low population sizes (Sadovy de Mitcheson 2016). However, a strong positive correlation was found between BSB harvest and fishery-independent adult density measures in southern California (Jarvis Mason et al. 2024), suggesting that our results were unlikely to be significantly affected by hyperstability.

Saltwater bass larval recruitment is sporadic

Consistent with BSB and KB young-of-the-year recruitment patterns derived from coastal power plant entrapment data (Miller and Erisman 2014) and dive surveys (Love et al. 1996, Jarvis Mason et al. 2024), our indices of larval abundance showed evidence of sporadic recruitment. Both basses showed the same few sporadic larval pulses, with the most recent in the mid-2010s (Fig. 2). Despite this similarity, overall, KB larval abundance was higher and less variable than BSB (Fig. S8), a pattern consistent with locally sourced larval recruitment (Selkoe et al. 2007) and higher densities (Cowen 1985). In contrast, BSB has high genetic diversity among individuals from southern California to Baja California (Paterson et al. 2015) and may be more influenced by El Niño-driven sporadic northward advection of larvae from Baja California (Jarvis Mason et al. 2024). This advection has been suggested for other fishery species in the region (Smith and Moser 1988, Allen and Franklin 1992, Ben-Aderet et al. 2020). Larval dispersal models demonstrated transboundary larval connectivity of up to 500 km between northern Baja California and the SCB in the fall for other species (Arafeh-Dalmau et al. 2023) with a similar pelagic larval duration to BSB (∼lunar month, Findlay and Allen 2002). Off Baja California, BSB spawning peaks in summer and fall, with higher larval abundance during El Niño events (Avendaño-Ibarra et al. 2009). In the future, it is predicted that connectivity between the SCB and Mexico will be significantly reduced, thus highlighting the importance of protecting habitats important for larval connectivity and binational conservation efforts to manage the collective resource (Arafeh-Dalmau et al. 2023).

SST, zooplankton, and ILD predict saltwater bass larval densities

We found strong relationships between saltwater bass larval density and SST, zooplankton biomass, and ILD. In general, our model results indicated CalCOFI bass larval density was higher when SSTs were warmer, zooplankton biomass was higher, and the ILD was shallower (Fig. 7). These relationships were relatively stronger for BSB than they were for KB, suggesting BSB is more sensitive to cooler temperatures and declines in zooplankton biomass. Indeed, BSB was historically considered a southern, subtropical/tropical species with an inconsistent presence in southern California (Young 1963, 1969, Frey 1971, Feder et al. 1974). These environment–species relationships are not surprising, as peak spawning for KB and BSB occurs during the warm summer months (Erisman and Allen 2005, Jarvis et al. 2014b) following the spring transition in the SCB (McClatchie 2014). We further hypothesized that increased spring copepod biomass available to adult bass prey species (e.g. coastal forage fishes) might translate into positive maternal conditioning of bass larvae in summer (Walsh et al. 2024). However, incorporating spring biomass of C. pacificus in our models did not improve fit.

In a prior analysis of environmental covariates on the biological response of larval fish communities in the California Current Ecosystem (CCE), the ILD was most predictive of community composition and ecosystem response within the southern California region (Hunsicker et al. 2022). In our study, ILD had the strongest environmental relationship with bass larval abundance. This relationship might be explained by ‘seasonal trophic amplification’ in nutrient-rich regions like the CCE (Xue et al. 2022). In these regions, phytoplankton and zooplankton become more concentrated within shallower MLDs, thereby increasing prey encounter rates and grazing rates of zooplankton, such that higher feeding efficiency in these regions is driven more by seasonal changes in the MLD than by the amount of food present (Xue et al. 2022).

Additional abiotic processes, including the timing, direction, and strength of ocean currents, can impact larval dispersal, connectivity, and settlement, particularly in the physically dynamic SCB (Cowen 1985). In the summer months, the counterclockwise southern California eddy forms within the SCB (Hickey 1993). The primary KB larval hotspot was the northern Channel Islands, which is located in the ‘transition zone’ where several cool and warm ocean currents converge (Hamilton et al. 2010). White and Caselle (2008) reported that in southern California, the density of giant kelp had a positive effect on KB larval settlement, but the relationship was conditioned on larval supply, whereby adult densities were largely a function of whether there was a match between areas of high larval supply and areas of higher kelp stipe density. In that study, the authors explain how eddies and temperature fronts may drive a match-mismatch between KB larval supply and kelp spore supply. This fine-scale spatial dependency may explain why Bight-wide kelp canopy extent was not an important driver of our spatially explicit KB larval densities.

Larval dispersal and retention are likely equally important for BSB. Nearshore flow in the SCB is mostly alongshore with a minor cross-shelf component (Winant and Bratkovich 1981, Barnett et al. 1984). Whereas we found BSB larvae were almost exclusively distributed nearshore along the mainland coast, KB larvae also occurred at the islands (Figs. 3 and 4). Indeed, we observed differences in the two species’ models’ spatiotemporal fields (ϵ_s,t), which indicated different latent forcing of dynamic biotic and abiotic processes (Figs. S9 and  S10). Latent variables not captured in our models (e.g. currents, eddies) may either help retain larvae close to shore (BSB; Barnett et al. 1984) or provide longshore and offshore dispersal to increase the probability of reaching areas with giant kelp, including at offshore islands (KB; see asymmetric directionality of KB larval transport in Watson et al. 2010). Differences in spawning aggregation formation/migration may also contribute to the observed spatiotemporal variability, as BSB form large spawning aggregations at a handful of predictable primary locations, while KB spawning aggregations are smaller and more broadly dispersed in both time and space (Erisman and Allen 2006, Jarvis et al. 2010).

No northward latitudinal shift in saltwater bass larval abundance

Given that historical bass populations extended much farther north than southern California during the region’s tropicalization of the mid-1800s (Hubbs 1948), we presumed secular ocean warming and increased heat content in surface waters within the study region (Gallo et al. 2019) may have resulted in a northward latitudinal shift in either or both species’ larval COG within the region through time. However, we found that temperature only partially explained the variation (26%) in annual COG northings of KB (but not BSB), and we did not find evidence of a northward latitudinal shift in either species’ COG northings over the 54-year period (Fig. 5). Northward larval dispersal of nearshore species beyond the geographic break at Pt. Conception in the SCB, at least during this study period, is thought to be limited more so by hydrographic features than temperature (Cowen 1985). However, it is important to note that our predicted annual COG northings were informed by larval occurrence only within southern California (survey coverage off Baja California post-1978 was inconsistent), and thus, temperature-related latitudinal shifts in their overall geographic larval distributions could not be explored.

Implications for fishery recovery and resilience

Our indices of larval abundance can inform saltwater bass stock assessments and provide insights into fishery recovery potential. For example, despite the importance of the basses to southern California marine anglers, currently no formal stock assessment or management strategy evaluation exists for either species. Our larval abundance indices could be used as a young pre-recruit index (Le Pape et al. 2020) or to ground truth an assessment's predicted recruitment deviations. We observed near all-time highs in both species’ larval abundance during the mid-2010s, one of the warmest periods on record (Jacox et al. 2018b). Even without a stock assessment, this pulse in larval abundance (Fig. 2) and the ability of bass larval abundance to predict future catches (Fig. 6) suggest that near-term catches for both species could rebound. Indeed, harvest increased in 2020 (KB), 2021 (BSB), and 2023 (BSB; CDFW 2024), and BSB harvest in 2023 was the highest in over a decade (though still low compared to average historical harvest), surpassing KB harvest for the first time during the same period (CDFW 2024). Fishery-independent measures of adult BSB and KB densities also show an upward trend in recent years (CDFW 2021, 2023).

Predicting future resilience of the saltwater basses under a changing ocean remains challenging. Continued global increases in MLDs (Somavilla et al. 2017, Sallée et al. 2021) could counteract surface warming benefits to BSB larvae by decreasing the concentration of prey items. MLD deepening also has the potential for shallow, nearshore habitats such as kelp forests becoming nutrient-poor (Parnell et al. 2010), resulting in negative impacts on KB larval settlement (White and Caselle 2008). Additionally, some environment–species relationships are not always stationary (White et al. 2019). Here, most environment–species relationships over the 54-year study period held up when only using the latter part of the time series that spanned 32 years (Fig. 7); however, the influence of temperature and zooplankton biomass on KB larval abundance was somewhat weakened in the latter time series. This time series began just before the 1988–89 North Pacific climate shift (Hare and Mantua 2000) after which, some species apparent responses to temperature changed (Litzow et al. 2020) and previously established physical (e.g. surface temperature) and ecological relationships with the NPGO (Di Lorenzo et al. 2008) became weaker (Litzow et al. 2020). This potential for nonstationary environment–species relationships highlights the importance of ongoing monitoring and reassessment. Furthermore, long-term data streams, such as those provided by CalCOFI, are invaluable tools for assessment of smaller-scale fisheries, which typically have limited resources for full-scale, advanced quantitative assessments. Future research into the role that maternal effects, prey type and condition, and externally sourced recruitment play in bass larval dynamics will likely provide additional, important ecosystem considerations for managing this culturally and economically important recreational-only fishery.

Acknowledgements

We are grateful to CalCOFI cruise participants, past and present, M. Human and S. Charter for help with accessing archival specimens, L. Bulkeley, K. Farno, L. Martz, R. Quaal, A.C. Salazar Sawkins, and V. Tang for help with sorting samples, and L. Bulkeley for help with identifying the formalin preserved Paralabrax spp. larvae. This manuscript was greatly improved by constructive comments and suggestions from three anonymous reviewers. Fish illustrations by Amadeo Bachar, Studio ABachar.

Author contributions

E.T.J.M.: conceptualization, funding acquisition, investigation, data curation, formal analyses, methodology, project administration, visualization, writing—original draft; W.W.: data curation, resources, validation, writing—review & editing; E.J.W.: methodology, validation, writing—review & editing; A.R.T.: conceptualization, funding acquisition, data curation, resources, validation, writing—review & editing; B.X.S.: conceptualization, funding acquisition, methodology, resources, supervision, validation, writing—review & editing

Conflict of interest

None declared.

Funding

This research was supported by a grant awarded by the California Ocean Protection Council (Proposition 84 Competitive Grant Program, Project R/OPCSFAQ-09) and administered by the California Sea Grant College Program, as well as by a NOAA Cooperative Institute for Marine, Earth, and Atmospheric Systems (CIMEAS) grant awarded to B.X. Semmens.

Data availability

Data and code pertaining to the sdmTMB models and cross-correlation analysis are available online in a GitHub repository: https://github.com/ETJarvisMason/Basses. We performed all analyses in R 4.3.2 (R-Core-Team 2023).

References