Spawning phenology of a rapidly shifting marine ﬁsh species throughout its range

P


Introduction
In temperate marine ecosystems, the rate and timing of seasonal warming and cooling sets the rhythm for biological activity, making these systems vulnerable to global climate change as it rapidly alters ocean thermal properties.Generally, the seasonality of midto high-latitude ocean temperatures is changing with spring warming advancing earlier and fall cooling occurring later.These phenological shifts are occurring more rapidly in the ocean than in terrestrial systems (Burrows et al., 2011).Large-scale shifts in the timing of seasons in the physical ocean can affect the phenology, the timing of biological events, and distribution of marine organisms (Harley et al., 2006;Richardson, 2008).For marine fish, responses to changing ocean seasonality can include earlier spawning migrations (e.g.Dufour et al., 2010) and advanced spawning phenology (e.g.Rogers and Dougherty, 2019), whereas responses to ocean warming can include poleward range shifts (Kleisner et al., 2016) and changes in population productivity (Free et al., 2019).
Fish spawning phenology should theoretically be under strong selective pressure to promote successful reproduction and larval recruitment (De-Camino-Beck and Lewis, 2008).The timing of spawning balances reproducing in favourable conditions for adults and for larval hatching and growth.For adults, optimal breeding conditions can be influenced by temperature, habitat quality, and local population density (Hendry et al., 2001), and the presence or absence of such conditions can influence annual reproductive output for iteroparous adults (McBride et al., 2015).For larvae, variation in survival is associated with the degree to V C International Council for the Exploration of the Sea 2021.This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/ licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.which ontogenetic timing (e.g.hatching and first feeding) matches the timing of seasonal increases in their prey, a phenomenon demonstrated in Atlantic cod (Gadus morhua; Kristiansen et al., 2011) and haddock (Melanogrammus aeglefinus; Head et al., 2005).When larval hatching does not overlap with high prey abundances, a mismatch can occur at a critical life stage and lead to low recruitment success for that year (Cushing, 1990).Spawning phenology is also influenced by the degree of seasonality experienced from low-to-high latitudes.Shorter growing seasons at higher latitudes necessitate accurate timing of spawning to allow substantial larval growth before the onset of winter, increasing the probability of early life survival (Conover, 1992).Thus, temperature typically has an immediate physiological effect and, therefore, is considered a proximal cue that sets the timing of reproductive development and spawning (Clark et al., 2005).For fish that experience seasonality, daylength can also be a relevant cue (Bromage et al., 2001) and can synergistically interact with temperature to set photothermal cue regimes for reproductive development (Pankhurst and Porter, 2003).Understanding the influences of seasonality on fish spawning is important because, as fish distributions move to follow the direction and rate of local climate shifts (Pinsky et al., 2013), they increasingly experience conditions where the relationships between temperature, daylength, and other seasonal cues differ from that in their historic range.Therefore, there is a need to investigate the interaction between shifting spawning phenology and distributional range shifts in marine species.
Marine fish with wide latitudinal spawning distributions will experience different environmental conditions depending on their spawning location, which can lead to intraspecific differences in spawning timing, duration, and reproductive output.In lower latitudes where growing seasons are longer, spawning duration is usually extended but fecundity per spawning event can be lower (Vila-Gispert et al., 2002).Conversely, in higher latitudes where growing seasons are shorter, spawning duration is truncated but fecundity per spawning event is higher and time between spawning events is shorter to garner total annual reproductive output comparable to the lower latitudes (Conover, 1992).For example, in a common reef fish (Pomacentrus coelestis), fish spawning at a higher latitude had both a shorter interspawning interval and higher clutch weight compared to their lower latitude counterparts (Kokita, 2004).Similarly, for largemouth bass (Micropterus salmoides), in a theoretical reciprocal transplant study, Garvey and Marschall (2003) found that the spawning strategies in the lower latitude fish were unsuccessful at higher latitudes, leading to significantly lower lifetime fitness.Thus, fish that are adapted to spawn across a wide latitudinal range will exhibit different spawning strategies to best optimize recruitment success throughout their entire range.However, for fish undergoing current poleward range shifts as a response to ocean warming (Fields et al., 1993;Pinsky et al., 2013), it is unknown if fish spawning at these new higher latitude regions will exhibit spawning strategies adapted for the shorter seasonal timeframe of acceptable spawning conditions at higher latitudes.
Lessons from ocean warming focused studies can help us understand how temperature differences across latitude may impact partitioning of energetic resources towards reproduction.This is particularly relevant because, although adult fish may have wide thermal windows making them resilient to climate change, actively spawning fish have a narrower thermal window and are significantly more sensitive to ocean warming (Dahlke et al., 2020).
As temperatures warm, fish energy demand increases to maintain higher metabolic rates (Po ¨rtner and Farrell, 2008).Fish in warmer regions, if unable to compensate through higher consumption, must compensate physiologically (Perry et al., 2005) by allocating relatively more energy to metabolic maintenance than to reproduction and growth (Holt and Jørgensen, 2015).Energetic condition of a fish can be linked with reproductive output (Wuenschel et al., 2013), and estimation of condition pre-, during, and postspawning can illuminate potential trade-offs between metabolic demand and reproductive output.
The US Northeast Shelf (US NES) is a broad continental shelf that extends from Cape Hatteras, NC, to Georges Bank and provides habitat for a plethora of economically important fish species (Hare et al., 2016).Climate change is affecting the shelf through rapid ocean warming (Chen et al., 2020), including in the Gulf of Maine where temperatures have risen more rapidly than in 99% of the world's oceans (Pershing et al., 2015), and shifting phenology to earlier spring warming and later fall cooling extending the duration of summer (Friedland and Hare, 2007;Henderson et al., 2017).Ocean warming and shifting phenology of the US NES has already led to distributional range shifts for many important fish stocks (Bell et al., 2015;Kleisner et al., 2017) and higher fish abundance correlated to the lengthening of the summer season (Henderson et al., 2017), respectively.
The US Northern stock of black sea bass (Centropristis striata), a demersal reef fish, inhabits the US NES from Cape Hatteras, NC, to the Gulf of Maine (Musick and Mercer, 1977;Roy et al., 2012).Black sea bass are protogynous hermaphrodites and show female-biased sex ratios except at larger sizes and older ages (Mercer, 1978).Spawning occurs inshore throughout their entire distributional range during the late spring to early fall months (Moser and Shepherd, 2009).Recently, there has been a dramatic northward shift in the black sea bass population centre of biomass (Bell et al., 2015;Kleisner et al., 2017) and successful spawning and recruitment in the Gulf of Maine, a region where they were once infrequent inhabitants (McBride et al., 2018).These changes are attributed to recent ocean warming (Slesinger et al., 2019).Together with a wide spawning distribution and range expansion, spawning conditions differ for black sea bass based upon latitude and temperature (Friedland and Hare, 2007;Richaud et al., 2016).For example, during the beginning of the spawning season ($May), there can be a 6-10 C difference between bottom temperatures in the northern vs. southern inshore spawning locations (Narva ´ez et al., 2015;Richaud et al., 2016).
The recent black sea bass range expansion towards higher latitudes is important to study through the lens of phenology.Black sea bass, at multiple life stages, are voracious generalist feeders preying upon fish and epibenthic invertebrates (Drohan et al., 2007); here, the timing of spawning and recruitment success in the newly expanded range could significantly alter existing ecosystems.For example, in the Gulf of Maine, juvenile lobsters are vulnerable to black sea bass predation, threatening the most valuable fishery in the United States (McMahan and Grabowski, 2019).Also, black sea bass are an important recreational and commercial fisheries species along the US NES [NEFSC (Northeast Fisheries Science Center), 2017], and proactive fisheries management could benefit from increased knowledge of black sea bass spawning timing and output as their biomass shifts northward.To assess the implications of range expansion and shifting phenology on spawning phenology and output, we collected black sea bass throughout their entire spawning season to encompass 1011 Spawning phenology of marine fish species developing, spawning capable and post-spawning fish as well as throughout their distributional range from locations spanning from Virginia to Massachusetts.Metrics related to black sea bass spawning timing, reproductive output, and indirect energy indices were obtained for individual fish.We explored (i) the drivers of spawning season phenology and tested whether spawning season duration decreased with increasing latitude; (ii) whether reproductive output was higher in the northern locations compared to the southern locations; and (iii) differences in indirect energy indices throughout the spawning season and distributional range of black sea bass.

Fish collection
Fish collections targeted four inshore locations, representing almost the entire distributional range of the US Northern stock of black sea bass, that included (from north to south) the US states of Massachusetts (MA), New Jersey (NJ), Delaware (DE), and Virginia (VA).While the northern edge of the black sea bass range is now in the Gulf of Maine, collecting fish from this region was not attempted because of sampling challenges resulting from localized and patchy distributions (McMahan, 2017).In addition, fish collection occurred $110 km north of the southern edge of the US Northern stock, Cape Hatteras, NC, to eliminate the chance of unknowingly collecting fish from the US Southern stock due to stock mixing in this region (Roy et al., 2012).Time of collections specifically targeted black sea bass when they were first arriving inshore for spawning through the end of spawning when black sea bass are known to start offshore migrations for the winter.Collections in some regions were also constrained in time by the opening and closing of state mandated fishing seasons.Because the focus of this study was on spawning capable fish, only adults were kept and all fish <190-mm TL were returned to the water.In MA, fish were collected both from the Massachusetts Department of Marine Fisheries ventless trap survey and via hook-and-line sampling aboard a recreational charter boat.Sampling occurred in 2018 from the end of May to the end of August.In NJ, fish were collected both through a joint Rutgers University and New Jersey Department of Environmental Protection artificial reef trap survey and via hook-and-line sampling in conjunction with another ongoing black sea bass research project at Rutgers University at the time.Sampling occurred in 2018 from late April to the early November, but a gap in sampling during June occurred as an artefact of the sampling plan for the artificial reef trap survey.In DE, fish were collected through the Delaware Natural Resources and Environmental Control protected area trap survey in 2019 (end of May to the end of October).And in VA, fish were collected through hook-and-line sampling aboard a recreational charter boat in 2019 (end of May to mid-September).For all hook-and-line sampling, collections and euthanasia were performed by Slesinger and in accordance with Rutgers University IACUC Protocol (# PROTO201800054).Sampling duration, frequency, collection methods, and state or federal issued fishing permits are provided in Supplementary Table S1.
Bottom temperature was measured during each sampling event.For collections using fish traps, each trap was fit with a bottom temperature logger (HOBO TidbiT v2 Water Temperature Data Logger, Onset Computer Corporation) to provide a single averaged temperature value for each sampling day.During hook-and-line sampling, a handheld profiling CTD device (SonTek CastAway-CTD, Xylem) was lowered to the bottom.

Fish processing
Black sea bass were kept on ice and dissected in the laboratory no longer than 24 h post-capture.During dissections, length (mm TL) and weight (g) were measured, and the gonads and liver were removed and weighed to the nearest 0.01 g.Black sea bass were sexed (male, female, transitional) and macroscopically assigned a maturity stage of either developing (DEV), spawning capable (SPC), or post-spawning (POST) following Brown-Peterson et al. (2011).An ovarian mid-lobe cross-section was also removed, placed in a histology cassette, and fixed in 10% neutral-buffered formalin.After a minimum of 1 month in 10% neutral-buffered formalin, the histology cassettes were rinsed with deionized water and transferred to 70% ethanol for preservation.Ovarian crosssections were processed for histology using a haematoxylin and eosin stain (Histology Consultants, Inc., Everson, WA, USA).Histology was used to QA/QC the maturity stages assigned macroscopically to ensure proper classifications, especially during portions of the spawning season when black sea bass can transition from female to male and where macroscopic staging can become difficult due to the small size of the gonads (Klibansky and Scharf, 2015).Transitional status was assigned when male tissue was overgrowing the female tissue, there was significant atresia in residual oocytes, and male tissue comprised >5% of the gonadal cross-section (Klibansky and Scharf, 2015).Maturity stages were also used to estimate spawning season characteristics.For each region, a threshold of above 0.25 proportion SPC was set to estimate the active spawning season length, and a threshold of above 0.75 proportion SPC was used to estimate a range of peak spawning dates.Several other threshold values were tested, and threshold choice did not significantly alter the estimates of spawning season length (Supplementary Figure S1).Gonadosomatic index (GSI), calculated as ([gonad weight]/ [body weight]) Â 100, can be used as a quick to measure proxy for reproductive output if maturity stages are known (Lowerre-Barbieri et al., 2011) and was used in our study because estimating total annual or batch fecundities of black sea bass would have required more frequent sampling (Klibansky and Scharf, 2018).To eliminate biases related to oocyte size (Ganias et al., 2014), only GSI measurements from SPC fish, where hydrated oocytes were present, were used.
Two indirect energy indices, relative condition factor (Kn) and hepatosomatic index (HSI) were calculated.These indices were chosen because Kn can fluctuate throughout the spawning season indicating energy usage for spawning or recovery in energy stores postspawning (Slotte, 1999), and HSI can estimate liver energy reserves (Lambert and Dutil, 1997) used to determine pre-spawning energy as many fish recruit proteins and lipids from the liver for gonadal development (Brown and Murphy, 2004;Rosa et al., 2020).Kn was calculated as the proportion of the measured fish weight to a predicted fish weight derived from a log-log length-weight relationship from all fish collected (Blackwell et al., 2000).To remove the influence of reproduction on Kn, gonad weight was subtracted from both the final measured and predicted fish weights.HSI was calculated as ([liver weight]/[body weight]) Â 100.

Statistical analyses
Generalized linear models (GLMs) were used to independently predict measures of spawning season phenology, reproductive output, and indirect energy indices.For each model, predictor variables were scaled to the mean of that predictor if they were included in an interaction, and the variance inflation factor (VIF) was used to assess the impact of correlations among predictor variables where a VIF score of >4 indicated collinearity between predictor variables (O'Brien, 2007).For each model, the best fit models of spawning phenology, reproductive output, and indirect energy indices were determined using the second-order Akaike's Information Criterion (AICc; Burnham and Anderson, 2002).For the reproductive output and spawning phenology models, collection method (trap vs. hook-and-line) was included as a factor to account for differing size and sex selectivity of the sampling methods (Provost, 2013).If collection method was a significant predictor, results were scaled to a standardized collection method (hook-and-line).For predictor variables included in the best fit models, a dominance analysis, using the "dominanceanalysis" package in R (Navarrete, 2020), was used to compare the relative importance of each predictor in the model.The dominance analysis compares R 2 values (McFadden R 2 for GLMs) for each predictor across all subsets of the model to determine the relative contribution of one predictor over the other (Azen and Budescu, 2006).This analysis can be useful when all or most predictors are significant in the best fit model to further explore the ability of a predictor variable to explain the variance of the response variable.All analyses were performed using R (Version 4.0.1;R Core Team, 2019).GLMs were performed using the "glm" function, maps were made using the package "ggmap" (Kahle and Wickham, 2013), and GLM partial effects were visualized using the package "Jtools" (Long, 2020).
Spawning season phenology was investigated by fitting two separate candidate models with a binomial distribution that described the beginning and end of spawning separately.The beginning of spawning model included SPC and DEV fish, and the end of spawning model included SPC and POST fish.To explore the effects of Julian day and daylength, two separate candidate models with either Julian day (1) or daylength (2) were fit to the beginning and the end of spawning models.Separate models were used due to collinearity between the two predictors, but both were assessed as they provided slightly different information.Julian day was constant across latitude and linearly increased throughout the spawning season, whereas daylength changed across latitude and responded non-linearly and nonmonotonically with time as the spawning season overlapped with the summer solstice. (1) In both ( 1) and ( 2), m is the probability of being SPC.Predictor variables were latitude (continuous), temperature (continuous), and collection method (categorical with two levels).The interaction term was latitude Â temperature.In (1), Julian day was included as a predictor variable and Julian day Â latitude as an interaction term.In (2), daylength was included as a predictor variable and daylength Â latitude as an interaction term.The best fit models that included Julian day and daylength separately were competed against each other using AICc to determine which metric, Julian day or daylength best explained the variation in the data.
For reproductive output, a candidate model (3) was fit to GSI from SPC fish to allow the analysis of reproductive output only during the time when fish were actively spawning.
The response variable (m ¼ GSI) was a non-negative ratio and modelled with a gamma distribution.Predictor variables were length (continuous), sex (categorical with two levels), location (categorical with four levels), and collection method (categorical with two levels).The interaction term was sex Â length to account for sex-specific differences in size.
The change in indirect energy indices of black sea bass was investigated in reference to maturity stage (DEV, SPC, and POST) as well as collection location.Both HSI and Kn were modelled separately: Here, m are HSI or Kn, to assess the indirect energy indices separately, but were both were non-negative ratios and modelled with a gamma distribution.For (4), the predictor variables were length (continuous), sex (categorical with two levels), location (categorical with four levels), and maturity stage (categorical with three levels).The interaction term was sex Â length to include potential sex-specific differences between length and the indirect energy indices.

Results
Throughout the spawning seasons of 2018 and 2019, a total of 898 black sea bass were collected in locations relevant to the current state of the black sea bass fishery and population, and spanning almost the entirety of the US Northern stock distributional range (Figure 1).The total number of fish collected in each location can be found in Supplementary Table S1.Because collections targeted black sea bass located inshore for spawning, most sampling occurred in relatively shallow inshore waters, except for in the southern locations (DE, VA).Both the DE state survey, which focused on protected reef areas, and the VA black sea bass fishery operates further offshore.Due to the apparent difficulty in obtaining adult black sea bass inshore, collections were not attempted inshore.These further offshore sampling sites led to deeper collections throughout the southern locations.Because of increasing depth at lower latitude collection sites, measured bottom temperature counterintuitively increased with higher latitudes.All locations experienced warming temperatures during late summer to early fall indicating the beginning of fall overturn as a result of increasing storm frequency (Rasmussen et al., 2005; Supplementary Figure S2).We collected a broad range of sizes of adult black sea bass in each location (Supplementary Table S2).On average, black sea bass in MA and VA were larger than in NJ and DE.While these size differences were likely attributed to the higher probability of collecting larger fish via hook-and-line sampling than in traps (Provost, 2013), they provide a representative size range of black sea bass throughout their entire range.Similarly, there were different proportions of male-to-female fish in each location, which was likely an artefact of sampling and should not be taken as a population-wide sex proportion for black sea bass.Histology analysis agreed with the macroscopic evaluation of ovarian maturity stage in 97% of the samples and was able to determine misidentified fish that were either female fish early transitioning to male or female fish on the cusp of being spawning capable or spent.While we attempted to only collect sexually mature fish, there was one fish that was immature and was not used in analysis.

General seasonal and regional patterns
Black sea bass spawning phenology differed both temporally and spatially (Figure 2).For each region, black sea bass reproduction followed representative patterns of asynchronous batch spawning fish whereby it was rare to collect 100% SPC fish on a single day during the active spawning season.Except for VA, the duration of sampling was sufficient to collect black sea bass pre-and postspawning as seen by sampling days with 100% DEV or 100% POST fish, respectively (Figure 2).GSI closely followed female and male gonadal development with a peak when the proportion of SPC fish first reached a threshold of >0.75 proportion SPC (Figures 2 and 3) and then decreased as fish began to enter post-spawning phases.This indicates higher reproductive output in the beginning of the spawning season, which is common for multiple batch spawners.The variation in GSI measurements from VA fish was likely attributed to the proportion SPC fluctuating throughout the spawning season (Figure 2).Altogether, this suggests that GSI was a good approximator of the spawning season and reproductive output proxy for black sea bass.By using the 0.25 proportion SPC threshold, the duration of the spawning season showed an increasing trend with decreasing latitude (Figure 3).In addition, the regions where the spawning season was shorter, the duration of when black sea bass were >0.75 SPC was also shorter (Figure 3).

Spawning phenology
Both the beginning and end of spawning phenology were analysed separately and models using Julian day or daylength were competed against each other.For predicting both the beginning and end of spawning phenology best fit models, Julian day performed better when included in the model than when daylength was included (beginning: Julian day AICc ¼ 492.94, daylength AICc ¼ 669.53; end: Julian day AICc ¼ 339.18, daylength AICc ¼ 358.87).The model that explained the most variation in the probability SPC for the beginning of spawning included all main effects of latitude, Julian day, temperature, collection method, and a latitude by temperature and latitude by Julian day interaction.One other model had a DAICc < 6 (Richards, 2007; Table 1); however, this model had a lower AICc weight and was not chosen (Burnham and Anderson, 2004).In the best fit model, all main effects and interactions were significant (p-value <0.01; latitudejday p-value <0.05) except for latitude (p-value >0.05), which was still included because of the significant interaction between latitude and temperature and latitude and Julian day.Model residuals are found in Supplementary Figure S3.Overall, Julian day was the dominant predictor, followed by temperature, and a latitude by temperature interaction (Supplementary Figure S4a).Specifically, the probability SPC during the beginning of spawning was higher in fish collected with hook-and-line than with traps (Figure 4a).The probability SPC increased with Julian day and was similar across all latitudes, with a slightly earlier start in the most southern regions of the sampling locations (Figure 4b).The effect of temperature on probability SPC differed across latitude showing higher probability of SPC with lower temperatures in the southern latitudes and higher temperatures in the northern latitudes (Figure 4c).This is likely an effect of the sampling depths and associated temperatures within those regions.
The end of spawning probability SPC was best explained in a model with latitude, Julian day, temperature, collection method, and a latitude by Julian day interaction.All main effects and interactions were significant (p-value <0.01; temperature <0.05).There were two other candidate models with a DAICc < 6 (Richards, 2007), but the best fit model was chosen based on AICc weights (Burnham and Anderson, 2004).One model that included all terms within the best fit model and a latitude by temperature interaction received only slightly less support (DAICc ¼ 0.953 and AIC weight of 0.331 vs. 0.533 for the model without this interaction).All models with DAICc < 6 are in Table 1 and model residuals are found in Supplementary Figure S5.The dominance analysis showed Julian day as the dominant predictor and followed by temperature (Supplementary Figure S4b).For collection method, there was a higher probability of collecting SPC fish via traps than hook-and-line sampling (Figure 4d).The probability SPC decreased with increasing temperature (Figure 4e), indicating that the end of spawning was occurring when the fall temperatures were starting to increase associated with fall overturn in the region.The interaction between latitude and Julian day exhibited differing phenology across latitudes for the end of spawning with spawning ending on an earlier date in the northern regions and later in the southern regions (Figure 4f).

Reproductive output
The reproductive output proxy, estimated from the GSI of SPC fish, showed regional and sex-specific differences (Table 2 and Figure 5).The best fit model to explain variation in reproductive output included sex and location and both terms were significant (p-value <0.01).Six alternative candidate models had a DAICc < 2 that included additional terms and marginally improved the log likelihood.Therefore, the most parsimonious model was chosen (Arnold, 2010).All models with a DAICc < 6 are reported in Table 2, and model residuals are found in Supplementary Figure S6.The dominance analysis showed that location was the dominant predictor following the sex of the fish (Supplementary Figure S7).From the best fit model, the sex of the fish had a negative coefficient value, indicating higher reproductive output in females than in males.For example, a male fish from DE was predicted to have an average GSI of 5.46 compared to a female fish from DE with a predicted average GSI of 6.71.Location for NJ, DE, and VA all had positive coefficient values with respect to MA indicating higher reproductive output when compared to MA (Figure 5).Here, a female fish from MA was predicted to have a GSI of 3.61 (compared to a female fish from DE with predicted average GSI of 6.71).

Indirect energy indices
The indirect energy index, HSI, varied throughout the spawning season and collection locations.The model that best explained variation in HSI included location, maturity stage, length, sex, and a length by sex interaction (Table 3), and all main effects and interactions were significant (p-value <0.01), except for length, lengthjsex, and location between MA and NJ (p-value >0.05).One other model, without the interaction term, had a DAICc < 6 but did not improve the log likelihood and had a lower AICc weight (Burnham and Anderson, 2004; Table 3).Model residuals are found in Supplementary Figure S8.The dominance analysis indicated that maturity stage was the dominant predictor (Supplementary Figure S9).Specifically, for each region, as maturity stage advanced from DEV to SPC to POST, HSI decreased.HSI was higher in female than in male fish and increased in southern locations (DE and VA; Figure 6).The length by sex interaction showed that HSI slightly increased with length for male fish but was relatively similar across sizes of females.
The relative condition factor (Kn) was largely uninformative for black sea bass and the values were relatively constant, remaining close to 1, throughout the spawning season and locations.As such, results for Kn are found in the Supplementary Material.

Discussion
Our study demonstrates important intraspecific differences in spawning timing and output of the US Northern stock of black sea bass, whereby spawning duration was shortest in the northern locations and longest in the southern locations.The beginning and end of spawning was mostly driven by Julian day and temperature.Reproductive output was lower in the northern locations and higher in the southern locations, which contrasts with previous studies that suggest reproductive output is highest in the higher latitudes to compensate for the shorter spawning duration.This result could be related to lower pre-spawning energy stores (e.g.HSI) in the northern fish when compared to the southern fish, potentially related to migration dynamics.These results suggest potential decreased recruitment in the northern regions and further investigation is warranted as black sea bass continue to expand their range poleward.

Julian day and temperature effects on spawning phenology
Black sea bass spawning duration decreased with increasing latitude, which corroborates well documented fish spawning patterns (Conover, 1992).The timing of the start and end of spawning was mostly driven by Julian day and temperature with interactive effects of latitude.In many fish spawning phenology studies, daylength and temperature are dominant drivers of the initiation and/or termination of spawning (e.g.Holt and Riley, 2001).However, in our models, Julian day outperformed daylength, which is intriguing because daylength has been shown to affect the timing of reproduction in black sea bass (Howell et al., 2003).This result could be an artefact of our statistical methodology and/or reproductive biology of black sea bass.The summer solstice (i.e.peak day length) occurred during a portion of the spawning season with high proportions of SPC fish; thus, there was relatively little contrast in day length within our data.Along these lines, the daylength cue may have also been mediated through Julian day (with which it was highly correlated, r 2 ¼ 0.68, p-value <0.001) as some fish can use daylength as a cue to record time (Bromage et al., 2001).Biologically, many studies that investigate the effects of daylength on spawning phenology are focused on the timing of the initiation of vitellogenesis (e.g.Clark et al., 2005).Our classification of the beginning of the active spawning season focused on when fish became SPC, or in other words, had initiated the oocyte hydration process which means these fish had already undergone vitellogenesis.For some fish, temperature can be used as a cue to synchronize the final stages of maturation and initiate the transition into postspawning phases (Van Der Kraak and Pankhurst, 1997), which is more aligned to the reproductive processes measured in our study.Moreover, daylength varied with latitude while Julian day was fixed across collection locations.The better performance of Julian day could indicate that population-wide timing of spawning may be more driven by internal biological clocks and/or proximal effects of temperature, a cue that did vary with latitude.b) the duration of the spawning season as designated by 25% of all fish SPC (black dot and lines), duration of when 75% of all fish were SPC (orange dot and lines), and the day when the mean GSI was highest (blue dot).For reference, Julian day 140 is 20 May and Julian day 260 is 17 September.
Table 1.Candidate models evaluated for their ability to explain the probability of spawning capable fish for the beginning and the end of spawning.Temperature has also been shown to significantly influence spawning phenology for temperate fish species including mackerel (Scomber scombrus; Jansen and Gislason, 2011), walleye pollock (Gadus chalcogrammus; Rogers and Dougherty, 2019), and striped bass (Morone saxatilis; Clark et al., 2005).While black sea bass spawning phenology was influenced by temperature, our results differed compared to these other species.It should be noted, however, that the VIF in the model predicting the start of spawning was moderately high ($4)-a reflection in the correlation (r 2 ¼ 0.51) between latitude and temperature.The collinearity suggests that we cannot completely separate the effects of latitude and temperature on the start of spawning.Spawning began in cooler temperatures in the lower latitudes and warmer temperatures in the higher latitudes whereby the initiation of spawning at our northernmost site (MA) occurred at a bottom temperature of $15 C, while the equivalent temperature at our southernmost site (VA) was $10 C.This would indicate that across latitude the start of spawning was not as temperature sensitive.A lack of a strong relationship between temperature and spawning contrasts with other fish that will spawn earlier in warmer water (Lyons et al., 2015).While apparent temperature insensitivity could make black sea bass resilient to climate change, it may disadvantage them if matching spawning to earlier warming would be more adaptive.Spawning ended across each Table 2. Candidate models evaluated for their ability to explain the variation in GSI from SPC fish.

Response variable
All models are fit with a gamma distribution and an identity link function.From all candidate models, only those with DAICc < 6 are shown.The relative likelihood of a model is represented by the AICc where an AICc ¼ 1 represents the most likely.
sampling location as fall bottom temperatures warmed such that across all sites there was an average warming of 7.25 6 4.19 C from the beginning to end of the active spawning period.This fall warming is largely a destratification event-driven process and occurs as warm surface water mixes with the Cold Pool, a cold water mass below the mixed layer, due to the onset of fall storms (Rasmussen et al., 2005).Fall mixing events have been shown to affect the black sea bass behaviour, whereby some fish began to migrate offshore (Secor et al., 2019).Overall, the relationship between temperature and black sea bass spawning phenology is complex, and further research is needed to elucidate the different responses of the beginning and end of spawning to temperature.

Reproductive output across spawning locations
For black sea bass in this study, the reproductive output proxy, GSI, was consistently lower for male fish than for female fish (Figure 6).Sexual dimorphism in reproductive output has been found in other fishes [e.g.blackeye gobies (Rhinogobiops nicholsii); Yong and Grober, 2014] and showed higher female GSI than male GSI.Lower male reproductive output in this study is of interest for black sea bass because sperm limitation has been raised as a possible concern for hermaphroditic species that develop highly female skewed sex ratios (Provost and Jensen, 2015).Across our sampling locations, there was lower GSI in the most northern location (MA; maximum GSI ¼ 7.53) and higher GSI in the more southern collection locations (e.g.DE; maximum GSI ¼ 14.53).GSI from black sea bass in the southern locations was also higher than in locations comprising their newly expanded range from Southern New England to Maine (McMahan et al., 2020).This contrasts with other studies where fish spawning at higher latitudes will compensate for a shorter spawning season through higher reproductive output (Conover, 1992;Kokita, 2004).The lower reproductive output of individuals at the black sea bass northern range margin suggests that these fish may not be fully adapted to spawning at capacity at higher latitudes compared to individuals from their core historical range.While no previous studies on black sea bass examined GSI across the entire distributional range of the Northern stock as we have done here, patterns within sub-regions or across studies generally contrast with our findings and suggest higher GSI in the more northern regions compared to the southern regions.For example, in Wuenschel et al. (2011), black sea bass collected from Massachusetts and Rhode Island (i.e., the northern part of their range) from May to June had maximum GSI values $10, while Wilk et al. (1990) found black sea bass collected south of Rhode Island along the New York Bight (near the centre of their range) to have maximum GSI values $8 during July.In Rosa et al. (2020), black sea bass from New Jersey in June had higher GSI than fish collected $250km south off of the coast of Maryland during July and August.In all of these studies, the maximum GSI values from the different locations were measured during different months, thus confounding spatial comparisons of GSI.For example, comparing GSI from black sea bass caught in New Jersey in June to black sea bass caught in Maryland in August may lead to a comparison of GSI during peak spawning vs. towards the end of the active spawning season, respectively.Therefore, comparison across studies or at different time periods does not allow for a full assessment of regional differences in GSI.Our GSI comparisons were made across the entire duration of Table 3. Candidate models evaluating the variation in the indirect energy index, HSI.   the spawning season which allows for a more accurate comparison of estimated reproductive output amongst regions.
While the results for black sea bass indicate a reduction in reproductive output at higher latitudes, these results should also be considered alongside the use of GSI as a proxy for reproductive output.A female black sea bass with a higher GSI should have more oocytes (Mercer, 1978), but the oocyte stage can nonlinearly affect the weight of the gonad (i.e.hydrated oocytes are disproportionately heavy; Ganias et al., 2014).To account for this, we only analysed GSI from a standardized maturity stage, where hydrated oocytes were present.There could be variation in the weight of the gonad if a female fish was in-between spawning batches and recently released hydrated eggs.However, because black sea bass are asynchronous spawners and should not all be within or between spawning batches at the same time, assessing the GSI from multiple fish caught on a single day should provide enough variation to use this proxy across collection dates.

Energy allocation throughout spawning
Lower reproductive output at the northern regions of our sampling locations was unexpected but may be explained through the differences in energetic condition of black sea bass throughout spawning.HSI was lower in the northern region (MA) in developing fish when compared to the other sampling locations.However, post-spawning fish in MA had higher HSI than in the more southern regions, indicating a recovery in liver energy stores in this northern region.Fluctuations in HSI can be linked to reproductive output because the liver is a main energy storage depot (Brown and Murphy, 2004) and provides hepatic vitellogenin to be used in developing oocytes at the beginning of ovarian development (Hiramatsu et al., 2002).Thus, lower HSI in developing fish could lead to lower reproductive output.Although MA was our northernmost site, due to regional oceanographic patterns, this region warmed earlier and was the warmest during spawning compared to the other regions (Supplementary Figure S2), potentially decreasing energy stores due to higher energy demand.In addition, the patterns in HSI of MA black sea bass could be a product of migration distances.Black sea bass overwinter near the continental shelf edge off the coasts of Virginia to Maryland and migrate inshore for the summer to spawn leading to longer migration distances for the higher latitude fish ($100-200 km for southern fish vs. $400-500 km for northern fish; Moser and Shepherd, 2009).For MA fish, lower HSI before spawning could be related to lower energy reserves following migration, and relatively higher HSI in post-spawning fish (compared to other collection locations) could indicate increased feeding to prepare for a long migration to return to overwintering sites.While out of the scope of this study, future research is needed to focus on energy and lipid partitioning pre-, during, and post-spawning in black sea bass throughout their distributional range to disentangle the effects of migration and local oceanography on energetics and subsequent impacts on reproductive output.

Trailing edge of US Northern black sea bass range expansion
Understanding the potential changes at the trailing edge of the US Northern stock of black sea bass distributional range is important as this region is subject to ecosystem restructuring and thermal habitat loss (Kleisner et al., 2017).While the spawning season in VA started earlier in the year and had the second longest spawning duration after DE (Figure 3), GSI was high but variable due to fluctuations in the proportion of SPC fish towards the end of the spawning season (Figure 2d).Collecting adult black sea bass inshore was not attempted due to the fishery operating further offshore, which lead to the nearest inshore sampling site to be 35 km from the coast.Black sea bass off of VA have historically been collected during the spawning season farther offshore ($35-50 km;Mercer, 1978;Moser and Shepherd, 2009), but in recent years, spawning black sea bass abundance has been low from inshore to mid-shelf (NorthEast Area Monitoring and Assessment Program survey data; http://fluke.vims.edu/fishgis/faovims/index.htm; J. Gartland, pers.comm.).In our study, black sea bass collections during the middle of the spawning season were made where there was a higher probability of catching fish, which led to sampling sites close to the shelf break ($100-140 km from the coast).Altogether, the US Northern stock of black sea bass distribution in the southern portion of their range has existed further offshore in the past, but increasing difficulty in obtaining black sea bass inshore may suggest a range contraction.This could be a result of ocean warming, where this region has experienced significant warming from 1977 to 2016 (Wallace et al., 2018), but further research is required to understand the effects of warming at the trailing edge of US Northern stock of black sea bass population.
It should be noted that our conclusions are not extended to the entire population of black sea bass.The US Southern stock of black sea bass, which inhabits waters south of Cape Hatteras, NC, to southern Florida, are considered genetically distinct and mostly reproductively isolated from each other (Drohan et al., 2007;Roy et al., 2012) and may have higher thermal tolerances (Atwood et al., 2001).Spawning phenology and ocean warming effects on the Southern stock of these fish and throughout their range are out of the scope of this study but would be valuable for further interpretation of the results for the US Northern stock of black sea bass.

Study limitations
Due to logistical constraints, black sea bass were collected across two separate years (2018 and 2019).While the dominant predictor of both the start and end of spawning was Julian day, suggesting that these results are comparable across the 2 years, temperature did have a significant, yet more minor, effect on the timing of reproduction.Interannual temperature variability occurs across the US NES and has been shown to potentially impact interannual recruitment success of black sea bass (Miller et al., 2016).Therefore, further investigation of spawning duration and reproductive output results across all regions was conducted.We assessed the difference in average monthly bottom temperatures for each spawning location between 2018 and 2019 using modelled bottom temperature from the Doppio Regional Ocean Model System model (Wilkin et al., 2018) and found similar ($D2 À 4 C) bottom temperatures between 2018 and 2019 in each collection location (Supplementary Figure S10).Also, there was agreement on the timing of peak spawning between our study and previous studies.For example, Wilk et al. (1990) found highest GSI during July of black sea bass collected from the New York Bight, which was similar timing for black sea bass in this study caught in a similar location (NJ).GSI data from fish collected off of the NJ coast from 2011 to 2013 also suggest peak spawning times around July (Supplementary Figure S11; Provost et al., 2017).While agreement in estimated peak spawning times for black sea bass is likely, reproductive output may differ interannually.For the US Southern stock of black sea bass, significant interannual differences in proportions of spawning capable females and annual fecundity have been found (Klibansky and Scharf, 2018).While annual fecundity and GSI are linked, the degree to which interannual variability in fecundity is reflected in varability in GSI is unknown.Our GSI measurements likely provide a more conservative estimate of reproductive output, but future studies would benefit from investigating interannual changes in US Northern black sea bass GSI across their entire spawning range.

Conclusions
With ongoing climate change in the US NES, the intersection between marine fish poleward range shifts and spawning phenology will be important to study for many economically and ecologically important fish species.Our results showed that black sea bass do exhibit different spawning phenology across their latitudinal spawning gradient, yet their reproductive output did not compensate for a shorter spawning season in the northern regions.Black sea bass in state waters are managed with statespecific quotas that are based on regional biomass.As abundance increases in the north and re-evaluations of state fishing quotas are suggested, there should be caution towards setting annual total allowable catch limits due to the potential variability in recruitment in these new regions and future impacts from ocean warming.This caution is also extended to the trailing edge of the US Northern black sea bass distributional range where more data are required to investigate the potential impacts of ocean warming on adult energetics and productivity.The outcomes of this study, and potential impacts on recruitment and fishing quotas, are applicable globally to other range shifting species.Many midto high-latitude fish species are currently undergoing range shifts as a response to ocean warming and have migratory behaviours related to spawning.As we found, even if a fish may be able to follow cues to shift the timing and duration of spawning, this does not ensure reproductive output will increase to compensate for new spawning phenology, leading to potentially lower recruitment in newly expanded ranges.While range shifts may be seen as adaptive for fish negatively impacted by ocean warming, these benefits may not be fully realized if there is reduced spawning output in these new regions.

Figure 1 .
Figure 1.Collection locations for black sea bass throughout the entire sampling period for Virginia (purple), Delaware (red), New Jersey (blue), and Massachusetts (green).The grey lines from thinnest to thickest widths indicate the 25, 50, and 100 m isobaths.The box within the inset map of the United States denotes the full study region.Note: where the 100 m isobath occurs also approximates the shelf edge for the US NES.

Figure 2 .
Figure 2.For (a) Massachusetts, (b) New Jersey, (c) Delaware, and (d) Virginia, the percent of fish in each maturity staging category is shown in the top panel, where the percent of fish for each sampling day is cumulative, and the GSI (mean 6 standard error) is shown in the bottom panel.For the maturity stages, DEV, SPC, and POST.

Figure 3 .
Figure 3.The spawning season characteristics are shown for each location in (a) the spawning duration (the Julian day difference between the beginning and end) and (b) the duration of the spawning season as designated by 25% of all fish SPC (black dot and lines), duration of when 75% of all fish were SPC (orange dot and lines), and the day when the mean GSI was highest (blue dot).For reference, Julian day 140 is 20 May and Julian day 260 is 17 September.

Figure 4 .
Figure 4. Plotted partial fixed effects from the best fit model for the beginning of spawning (a-c) and end of spawning (d-f) GLM are shown.The beginning of spawning model best explains the probability of a fish being spawning capable [p(SPC)] over being developing based on (a) the single fixed partial effect for collection method, (b) the interaction between latitude and Julian day, and (c) the interaction between latitude and temperature.The end of spawning model best explains the probability of a fish being spawning capable [p(SPC)] over being post-spawning based on the single fixed partial effects for (d) collection method, (e) temperature and (f) the interaction between latitude and Julian day.In (a) and (d), HL ¼ hook-and-line and T ¼ trap.For (b), (c), and (f), the predictor variables are plotted on the y-and x-axes and the estimated probability SPC is shown in the heat map with darker colours indicating lower probability and lighter colours indicating higher probability SPC.For each partial effect plot, a higher probability SPC indicates that for (a)-(c) a fish will have started spawning, and for (d)-(f), a fish has not yet entered a post-spawning stage.
with a gamma distribution and an identity link function.Only models with DAICc < 6.The relative likelihood of a model is represented by the AICc where an AICc ¼ 1 represents the most likely.Stage: maturity stage.

Figure 5 .
Figure5.The reproductive output estimated from GSI measurements from SPC fish, for each location, is shown with female and male fish separated (female ¼ teal, male ¼ yellow).The violin plots show both the distribution and structure of the raw GSI data.The points and associated error bars are estimated GSI and standard errors, respectively, from the best fit GLM.

Figure 6 .
Figure 6.HSI for each region is shown with (a) female and (b) male fish plotted separately.The violin plots show both the distribution and structure of the raw HSI data.The points and associated error bars are the estimated HSI and standard errors, respectively, from the best fit GLM with the length of the fish kept constant to the mean length of all fish.For each region, HSI is separated by DEV, SPC, and POST.