Exposing changing phenology of fish larvae by modeling climate effects on temporal early life-stage shifts

Changing environmental conditions are influencing the seasonal timing in life history events of organisms. Such shifts in phenology are often linked to increasing temperatures that stimulate faster developments or earlier arrivals. This phenomenon has been demonstrated in terrestrial and aquatic realms, but data and knowledge are limited on how early life stages of fish are affected over long-term and broad environmental scales. Here, we analyze 2 decades (1974−1996) of size class-specific Baltic herring Clupea harengus membras L. larval data along the whole coast of Finland to expose shifts in phenology linked to changes in environmental covariates. We use a novel Bayesian hierarchical spatio-temporal hurdle model that describes larval occurrence and abundance with separate processes. Abundances are modeled with the Ricker population growth model that enables us to predict size-specific larvae groups in relation to the environment while accounting for population density dependence. We quantify shifts in phenology at multiple life stages, based on first appearances of smallest larvae (<10 mm) and by detection of higher proportions of larger larvae (>15 mm) appearing earlier than they have done historically. Our results show a strong signal in shifting phenology of the larvae toward an earlier development of 7.7 d per decade. Increasing temperature had a positive effect on the earlier development of larger larvae. Additionally, we highlight that the survival of larvae becomes more density dependent as their size increases. Our modeling framework can reveal phenological shifts of early life stages in relation to environmental change for survey data that do not necessarily cover the onset of reproduction.


INTRODUCTION
Long-term changes in environmental conditions due to progressing climate change alter species abundances and their distribution (Last et al. 2011, Bellard et al. 2012. The seasonal timing of recurring biological events, i.e. phenology, plays an important role in linking species with their abiotic and biotic surroundings, ranging from individual physiological responses to interspecific relationships (Forrest & Miller-Rushing 2010, Burthe et al. 2012. Changes in the phenology of species are considered to be the fin-gerprint of climate change while most often being linked to increasing temperatures. Such changes favor earlier occurrences or faster developments (Parmesan & Yohe 2003), which may result in mismatches of resource availability and lower recruitment success (Durant et al. 2005, Thackeray et al. 2010, Asch et al. 2019, Régnier et al. 2019. Studies addressing changes in the phenology of species typically include shifting firsts of ecological events, such as first budding and flowering of plants or arrival and departure of migrating species (e.g. Hopkins 1918, Fitter et al. 1995, Sparks & Menzel 2002. Although there is a body of literature also investigating phenological changes in aquatic realms, it remains comparably small due to the suitability and accessibility of data to study such changes (Thomas et al. 2014), i.e. long-term and high-frequency observations. While much is known about the ecology and biology of commercially important fish stocks, including their spawning seasons, the first days of spawning or hatching of fish larvae are difficult to record as well as being subject to annual variation. Hence, long-term shifts in such timings are rarely documented (but see Asch 2015).
Here, we developed a modeling framework to study possible changes in fish larval phenology through advancements in egg-hatching and developing times and applied it to spring-spawning Baltic herring Clupea harengus membras L. Herring comprise ecologically and commercially important species in the Atlantic and Pacific oceans (Blaxter 1985, Levin et al. 2016 as well as in the Baltic Sea (Parmanne et al. 1994, Natural Resource Institute Finland [Luke] 2020). For some Baltic herring stocks, there is evidence of earlier spawning and higher larval survival in years of early ice breakup and early increases of sea surface temperature (Arula et al. 2014(Arula et al. , 2016, but there is currently no study that systematically addresses phenological changes in this species. With the semi-enclosed and shallow Baltic Sea being among the fastestwarming large marine ecosystems globally (Belkin 2009), we hypothesize that the temperature-dependent spawning and growth cues have led to a progressively earlier development of Baltic herring larvae.
We test our hypothesis with data from a long-term monitoring program that surveyed Baltic herring larvae over 22 yr (1974−1996) at 7 different areas along the coast of Finland (Parmanne 2001) (Fig. 1). Since no information on the exact first dates of larval hatching are available from the area, we developed phenology measures based on the probability of occurrences of the smallest larvae (<10 mm) and on the fractions of 3 different size classes (<10, 10−15 and >15 mm) throughout the sampling season from May to August. Earlier dates in high occurrence probabilities of the smallest larvae indicate relatively earlier spawning, and earlier dates in high fractions (over 10% of the total larval pool) of the largest larvae (>15 mm) indicate faster development. Our method is based on a spatio-temporal hurdle model (Cragg 1971, Potts & Elith 2006, which can distinguish the environmental and density-dependent effects on both occurrence and abundance patterns separately. The model also enables us to link phenology changes to environmental covariates and to predict size-specific Baltic herring larval groups. To increase the resolution of our environmental covariates, we constructed spatiotemporal maps through spatio-temporal interpolation from a large number of in situ observations. Our results expose shifts in herring larval phenology and show size class-specific responses in occurrence and abundance while highlighting the long-term environmental drivers and density dependence effects influencing Baltic herring larvae.

Study area and survey data
Our study area covers Finnish coastal regions in the northern Baltic Sea between 59° and 65°N latitude (Fig. 1). The included latitudinal range is characterized by strong environmental gradients in water temperature, salinity and primary production. Salinity values range from ca. 6 to 7 psu in the Archipelago Sea in the south to under 3 psu in the Bothnian Bay in the north and the eastern Gulf of Finland in the east. While the northern areas are less affected by nutrient loadings, the Archipelago Sea and the Gulf of Finland are eutrophied due to the highly populated and agriculture-dominated catchment areas in the vicinity. During 1974 to 1996, a Baltic herring larval survey was conducted at 7 areas along the Finnish coast, with samples taken approximately every tenth day from May to August ( Fig. 1) (Parmanne 2001). Herring is a pelagic species that spawns in shallow coastal areas, where its adhesive eggs attach to littoral substrates such as macrophytes. After hatching, the yolk-bearing larvae (<10 mm) remain concentrated in the vicinity of the spawning beds, whereas the larger larvae (14− 18 mm) are more homogeneously distributed in the coastal areas (Polte et al. 2017). Once the larvae have finished their metamorphosis to juvenile fish, they return to the open sea. Each area was sampled at 5 (Bothnian Bay and western Gulf of Finland) to 6 (all other areas) locations, where transects were conducted using a modified Gulf V larvae sampler with a 300 μm collection mesh and an opening diameter of 20 cm. The sampler was pulled by 7 to 11 m boats at a speed of 4 knots (2 m s −1 ) during daytime. Samples were taken from each 1 m layer, and maximum depths at the various sites ranged from 5 to 18 m (for more details see Parmanne & Lindström 2003). Larval abundances were the highest in the Archipelago Sea (Areas 3 and 4, Fig. 1) and by far the lowest in the Bothnian Bay (Area 7, Fig. 1). For detailed tem poral size class-specific abundances, see Fig. S1 in Supplement 1 at www. int-res. com/ articles/ suppl/ m666 p135 _ supp/

Environmental covariates and spatio-temporal covariate model
We selected surface water temperature (°C), salinity (psu) and chl a (μg l −1 ) (as a proxy for primary production and hence resource availability) as covari-ates to explain the environmental contribution in the variation of Baltic herring larvae. All environmental covariates were constructed so that they summarize the spatial average over the sampling area (including all sampling transects of respective areas) at the date of larval sampling. Hence, the covariates vary between sampling areas and sampling occasions, which increases the information on the effect of covariates upon the phenology.
To produce these spatial averages, we first combined 10 000 point observations for each covariate from or within the vicinity of the sampling areas.
These observations were derived from the open data service of the Finnish Environment Institute (VESLA: Finnish Environment Institute and the Centres for Economic Development, Transport and the Environment, http:// rajapinnat. ymparisto. fi/ api/ vesla/ 2.0/; Baltic Environmental Database: Baltic Nest Institute, http:// rajapinnat.ymparisto.fi/api/veslabnirajapinta/ 1.0/) and data of Parmanne (2001), who had measured the temperature in fixed locations close to the Baltic herring larval sampling transects (see Fig. 1). We then constructed a spatio-temporal statistical model (a Gaussian process, Banerjee et al. 2014) to interpolate the covariate values over each study area in a grid with a 10 × 10 km spatial resolution. The grids varied by area from 60 to 90 km 2 . We set the spatial resolution with respect to the estimated spatial smoothness of the covariate models and the frequency of the point observations. Next, we averaged the predicted values over each study area to derive a study area-specific covariate value. We created such predictions for each time point where there was herring sampling. The benefit derived from this approach is that it explicitly accounts for spatio-temporal correlations in measurements and filters out random variations in environmental data that arise from local phenomena (such as daily weather) or measurement errors. Supplement 1 gives a detailed summary of environmental data and the interpolation approach (Section S1). Model fitting and predictions were done with GPstuff (Vanhatalo et al. 2013) in Matlab 2018b. Supplement 2 provides all data and code to reproduce the spatio-temporal covariate model. Since temperature is a key driver of possible shifts in phenology, we highlight the signal in increasing temperature of each month (May− August) over time (Fig. 2).
We also included spawning stock biomass (SSB) data of the Baltic herring stocks referring to ICES subdivisions 25 to 29 and 32 (ICES 2018) and subdivisions 30 and 31 (ICES 2017). We assigned yearly SSB values to the respective study areas. SSB values for the northernmost subdivisions 30 and 31 ( Fig. 1) were only available from 1980 onwards; hence, the areas located within them were modeled only from 1980 onward. While the coarse resolution of SSB may not accurately reflect regional differences within areas in the subdivisions, it represents the best available data covering the spatio-temporal extent of this study and was hence considered to be the best proxy for regional stock sizes.

Population model for Baltic herring larvae
Since the processes underpinning larvae presence and abundance may depend differently on the environment and SSB, we developed a hurdle model which fits both components separately. While a hurdle model can be viewed as a single model, it is fitted with 2 separate components: one modeling presence−absence and the other modeling abundance conditional on presence. Here, we model the presence of larvae using a generalized linear probit model and the abundance, conditional on presence, using a Ricker-type population growth model (Hilborn 1985, Tang 1985, Subbey et al. 2014. Moreover, since the environment and competition may affect different-sized larvae differently, we constructed the model for 3 size classes: small (<10 mm), medium (10−15 mm) and large (>15 mm). We denote by R j(i),t,c the number of larvae, by N i,t the SSB and by x j(i),t the environmental covariates, with subindexes i for stock assessment region (see Fig. 1), j(i) for sampling area within the ith stock assessment region, t for sampling date and c for size class. The hurdle model is then: (1) where π j(i),t,c is the probability of presence of larvae and μ j(i),t,c is the expected log larval production rate (i.e. log of the number of larvae per m 3 water per unit SSB), with σ 2 c being the variance of the independent random error in log larval production rate.
As shown in Supplement 1 (Section S2), the expected log larval reproduction rate can be derived from the first principles of a Ricker model and represented as a linear mixed-effects model: ( 2) where α j(i),c is an area-specific intercept corresponding to the log-relative proportion of SSB in sampling area j, ƒ c (x j(i),c ) is a function describing the effect of environment on the rate of reproduction and larval  (5), June (6), July (7) and August (8). Dots are spatial average temperatures over the sampling area at dates of sampling. Trends are highlighted with a linear fitted regression growth and survival until size class c, β j(i),c < 0 is an area-specific regression coefficient corresponding to the density-dependent decrease in reproduction success of an individual (up to a given larval size class) and ε j(i),t,c is a spatio-temporally autocorrelated random effect. The random effect ε j(i),t,c explains both process stochasticity and residual error due to noisy data. Since we have our size-class specific parameters, the model allows different environmental and density-dependent processes for different size classes. For example, if β j(i),c was close to zero for the smallest size classes and significantly negative for larger size classes, there would be density-dependent competition only between larger larvae. We modeled the probability of presence of larvae with probit regression with a similar regression function as in the log larval production rate, that is π t,c) . However, the ecological interpretation of the parameters of this model component is different from the Ricker model: α j(i),c explains the mean occurrence probability, ƒ c (x j(i),t ) + β j(i),c N i,t describes the effect of environment and SSB to occurrence probability and ε j(i),t,c accounts for spatio-temporally correlated random processes in a sampling area over the spawning seasons. In the occurrence probability model, the parameter β j(i),c is not restricted to be negative, and it does not correspond to density-dependent competition.
We modeled the covariate effects in both model components with polynomials of the form ƒ c ( As a marine species in the brackish water conditions of the Baltic Sea, Baltic herring is at its physiological distribution limit with respect to salinity (Illing et al. 2016). Therefore, we assumed a linear relationship of salinity and Baltic herring larvae. The quadratic response functions for chl a and surface water temperature encode an assumption that the observed conditions may include optimal levels of these covariates (Peck et al. 2012, Abdellaoui et al. 2017). Prior to the analysis, we log transformed the chl a data.
When fitting the multiple size class model to data, we applied a hierarchical Bayesian approach. The fixed effects (β j(i),c , b c,sal , …, b c,chl a2 ) were given a hierarchical prior such that, for example: where m j(i) is the mean density-dependent mortality parameter over size classes in sampling area j(i). The other fixed effects had analogous priors. Analogous to Hartmann et al. (2017), we modeled the intercepts α j,c as study area-level multivariate random effects with between-size class correlations and the spatiotemporal random effects ε j(i),t,c as per sampling area independent multivariate temporal Gaussian processes with between-size class correlations. That is, t,3 ]. T denotes the matrix/vector transpose and t denotes time. The benefit from using hierarchical priors over size classes is that they allow information flow between them.

Statistical analysis
With the Ricker model not covering zero larval counts, the above hurdle model (Eq. 1) factorizes with respect to the occurrence and abundance models so that we can conduct statistical inference for these 2 model components independently. However, as we aim to model all size classes simultaneously, we need to implement them as multivariate generalized linear mixed models. For this, we used the hierarchical modeling of species communities (HMSC) (Ovaskainen et al. 2017, Tikhonov et al. 2020, from the HMSC R package (Tikhonov et al. 2019, which allows us to model the correlations between size class-specific responses to environmental covariates and SSB as well as the correlations between size class-specific random effects. We fitted the models using the Markov chain Monte Carlo method implemented in HMSC. We examined all simulations and ensured convergence with the Gelman-Rubin potential scale reduction factor. We assessed the explanatory and predictive powers of the 2 hurdle model components separately, the presence−absence model through size class-specific area under the curve (AUC) (Pearce & Ferrier 2000) and coefficient of discrimination (Tjur R 2 ) (Tjur 2009) values and the abundance model by R 2 . To compute the predictive power, we performed a 2-fold crossvalidation for both the presence−absence model and the abundance model. We further followed Ovaskainen et al. (2017) to calculate the proportional variance partitioning among the fixed and random effects of each larvae size class as well as overall for both the presence−absence model and the presenceconditioned abundance model. Herring analyses were performed in R version 3.6.3. Supplement 3 provides the code for fitting and post-processing of the HMSC models.

Phenology measures
Due to the lack of exact information when the first individuals of fish stocks spawn and subsequently the first eggs hatch, we took advantage of the detailed size class distribution of Baltic herring larvae. We studied (1) whether the high probability of small larvae (<10 mm) presence had shifted to earlier dates and (2) whether the higher expected proportion of large larvae (>15 mm) were appearing earlier in the year over time. For (1), we recorded the first day of the year where the probability of occurrence of small larvae was greater than 75 and 90% for each year from 1974 to 1996. We chose these probability thresholds to ensure that general patterns are consistent regardless of the threshold. The higher threshold at 90% comprises only high probabilities of occurrences and therefore corresponds to a conservative measure for the presence of larvae. The lower threshold at 75% probability represents a less stringent cutoff and is, thus, more sensitive to changes in phenology. Shifts toward earlier appearances in the high proportions of large larvae indicate faster development and/or earlier spawning dates, even if the actual spawning dates are outside of the sampling time frame. Hence, for this phenology measure (2), we calculated the posterior probability that the fraction of large larvae is higher than 10% of the total larvae pool for each year from 1974 to 1996. We set a >10% abundance threshold rather than a first observation, to ensure the robustness of our findings throughout the sampling season, as the springspawning Baltic herring spawns over several months (e.g. Aneer 1989). Herring larvae are approximately 7 to 9 mm in length at hatching (Geffen 2002), and daily larvae growth rates in Finnish waters have been reported to range between 0.18 and 0.52 mm per day (Hakala et al. 2003); hence, it takes approximately 10 to 40 d for them to grow to the large size class after hatching. If spawning rate remained constant and all larvae survived, 10% of all larvae would be in the large size class earliest at the moment when the first spawned larvae had grown from the small to the large size class. If the growth rate increases, the proportion of large larvae reaches 10% earlier. However, since larvae experience mortality, the 10% threshold can be seen as a conservative measure for the development. As the growth rates of fish larvae are heavily influenced by temperature (Fey 2001, Hakala et al. 2003, we additionally relate the probability of the fraction of large larvae to the realized temperature gradient to highlight the impact of surface water temperature on larvae growth.

Model evaluation
We assessed the performance of the occurrence and abundance components of the hurdle model individually ( Table 1). The presence−absence model component was assessed with AUC values that indicate how well the model was capable of distinguishing between classes (0 and 1). Mean AUC values of 0.93 among the size classes suggest a high model performance. The overall explanatory power of the presence−absence model, reflected in Tjur R 2 , was 0.51, and the 2-fold cross-validation-based predictive power was 0.45. The presence-conditioned abundance model had a mean explanatory power (R 2 ) of 0.56 and a predictive (cross-validation) power of 0.43.

Estimation of Baltic herring larval drivers
We partitioned the explained variation to the fixed and random effects (Fig. 3a,b) and summarized the weights of the covariate effects (Fig 3c,d) in the presence−absence and abundance models separately. The temporally explicit spatial random effect (combination of year, day and area; ε j(i),t,c ) explained most of the overall variation for presence− absence and abundance of Baltic herring larvae as well as similar amounts in both submodels. It showed decreasing importance for larval presence with increasing larval size class. No such change was detected for larval abundances. While the proportion of explained variation attributed to the area-level random effect (α j(i),c ) remained similar among the size classes in the presence− absence model, it changed considerably among size classes in the abundance model, with decreasing influence of area from small to large larvae. SSB at area showed consistent amounts of explained variation for occurrences of each size class. Moreover, SSB was the most important fixed effect for abundances, and its relevance increased with larval size (small = 10.1%, medium = 28.8%, large = 42.7% of explained variation), reflecting that the density dependence of larval survival in the Ricker model is more important for larger larvae. The importance of environmental covariates differed markedly between the 2 model components. Over one-quarter of the explained variation could be attributed to temperature, chl a and salinity together in the presence−absence model, while they accounted roughly for 10% in the abundance model. In the presence−absence model, temperature and chl a separately accounted for an overall similar proportion of explained variation, approximately 11%, but with different amounts for varying size classes. The importance of temperature increases with size reflected in the explained variation at each size class (small = 5.4%, medium = 7.6%, large = 19.4%), while the explained variation attributed to chl a was relevant particularly for the medium size class (small = 9.2%, medium = 18.1%, large = 5.8%). Temperature and chl a played lesser roles for the abundance of larvae, with a total average of explained variation of 2.7% for temperature and 3.1% for chl a. Salinity showed opposing trends for the size classes, being most relevant with a positive effect for presence− absence of large larvae and showing a negative effect for abundance of small larvae (Fig. 3).

Shifts in phenology
First appearances in which small larvae were present with > 90% (> 75%) posterior probability have on average decreased by 4.6 d (7.9 d) per decade (Fig. 4, all areas). With the exception of the eastern Gulf of Finland at the 90% probability threshold, there was a consistent trend in the earlier days of larval presence over time over all areas. The northernmost area, in the Bothnian Bay, did not show any trend, as all predicted occurrences had a lower probability than 90%.
However, we found strong patterns showing earlier days of the year when considering the predicted development of large larvae fractions over time (Fig. 5). When averaged over all areas, the proportion of large larvae being over 10% of the total  1975 1980 1985 1990 1995 Area 2: Western Gulf of Finland 1975198019901995 Year Area 6: Kvarken 1975 1980 1985 1990 1995 Area 3: Åland (islands) 1975198019901995Year Area 7: Bothnian Bay 1975198019901995 Area 4: Archipelago Sea 1975 1980 1985 1990 1995 All areas combined Upper regression line equation and correlation coefficient is the linear regression fit for the first day of > 90% posterior probability per year, and lower regression line is the same for 75% posterior probability. Area 7 does not show trends since the data did not meet the threshold criteria. Color code represents the area colors in Fig. 1 abundance was reached 7.7 d earlier per decade. Area-specific changes ranged between 1.4 and 13 d per decade. The trend was weakest in the northernmost area of the Bothnian Bay but steepest in the 2 other northern areas, Kvarken and the Both nian Sea (Fig. 5). Higher surface water tem-perature was positively associated, with an increasing probability for large larvae comprising more than 10% of the total larvae abundance. While this observation holds true for all areas, the 2 northernmost areas (Bothnian Bay, Kvarken) showed the weakest increase (Fig. 6).  1975 1980 1985 1990 1995 Area 2: Western Gulf of Finland 1975198019901995 Year Area 6: Kvarken 1975 1980 1985 1990 1995 Area 3: Åland (islands) 1975 1980 1985 1990 1995 Year Area 7: Bothnian Bay 1975 1980 1985 1990 1995 Area 4: Archipelago Sea 1975 1980 1985 1990 1995 All areas combined Color code represents the area colors in Fig. 1

DISCUSSION
Our results indicate that Baltic herring have shifted their phenology toward earlier larval development over time. This finding is centered on 2 complementing model-based results. First, the earlier detection of first Baltic herring larvae (<10 mm) occurrences at high probabilities (Fig. 3) suggests that northern Baltic herring stocks spawn and/or larvae hatch earlier. Second, fractions of large Baltic herring larvae (>15 mm), contributing to more than 10% of the total larvae abundance, are reached on average 7.7 d per decade earlier (Fig. 4), suggesting faster development. While both results essentially support the same conclusion, the latter especially provides novel solutions for modeling phenological change in larvae and early life stages of fish and other animals. Since seasonal fish larval surveys were designed to cover the main spawning times of a stock, it is likely that the sampling time frame does not match the first spawning events. However, by utilizing size class-specific information, we can analyze the temporal change in the fraction of large larvae that can reveal changing phenology, as those early detected large larvae hatched and grew already a few weeks before the sampling. Our results show that temperature was the strongest driver in terms of explained variation of occurrences of large larvae in our model (Fig. 3a) and that it had essentially a linear response, with a positive effect on all size classes in both models (Fig. 3c,d). With growth rates being a function of temperature (Oeberst et al. 2009), we further highlight that increasing temperatures are leading to higher probabilities of detecting fractions of large larvae earlier, pointing toward an earlier larvae development.

Drivers of Baltic herring larval occurrences and abundances
Using a hurdle model framework, we disentangle and quantify the drivers of Baltic herring larval occurrences as well as abundances separately (Fig. 3). While the climate-governed environmental drivers (surface water temperature, chl a and salinity) were important covariates in explaining occurrence patterns, SSB per area contributed most of the explained variation of the fixed effects in both models, for larval occurrences and abundances but particularly for the latter (Fig. 3a,b). With SSB corresponding to the density dependence effect in the abundance model (Ricker model), our results indi-cate that density dependence is stronger in larger than in smaller larvae, showing increasing importance with increasing size class (Fig. 3b,d).
The environmental structuring is mainly determining if Baltic herring larvae are present or not. One explanation for this could be the temperaturedependent prey availability during the transition from small yolk-sac larvae to medium-sized larvae, when autonomous feeding starts (ca. 8−12 mm) (Rannak & Simm 1979). The importance of chl a, a proxy for prey availability (Garrido et al. 2008, Druon et al. 2019, was particularly pronounced for the small and medium size classes determining occurrences. In creasing trophic mismatches have been linked to future warming scenarios and can result in lower recruitment success (Régnier et al. 2019). In cases where temperature determines prey supply for the smaller size classes, and hence their survival, the occurrence patterns of large larvae could be indirectly affected by the effects on smaller larvae. The observed increasing importance of temperature with size ( Fig. 6) is in line with this assumption and reflected in the positive weight of the temperature effect (Fig. 3c,d). The weaker effect of temperature on the probability of abundance of large larvae in the northernmost area, the Bothnian Bay (Fig. 6), suggests an overriding effect of low salinity. With lowest salinity values of < 3 psu, compared to 6−7 psu at the southernmost in cluded areas, salinity may be the limiting factor for successful reproduction, with recruitment being lower with decreasing salinity (psu < 3) (Illing et al. 2016). Hence, the low salinity may also account for the absence of occurrences with high probabilities we observed in the Bothnian Bay (Fig. 4).
Archipelago zones comprise fine-scale environmental gradients that are difficult to capture. With our modeling approach for the environmental covariates not being designed for capturing such fine-scale gradients, we more importantly focused on capturing the within-year and between-year variation in the covariates rather than the accurate difference between sampling site-specific values. The resolution in the measured salinity, temperature and chl a data did not allow accurate fine-scale (i.e. sampling site-specific) prediction of these within-year and between-year variations. For this reason, we summarized the environmental covariates as spatial averages over each of the 7 sampling regions. That is, the daily covariate values attached to each sampling site and year were the daily averages over the region where that sampling site was located -and all sampling sites within the same region had the same daily values. Hence, our covariate values capture the within-year progression and the variation between years in the average covariate values over the sampling sites. The effect to herring larvae from between-site variation in covariates within each sampling region was then captured by the random effects of the larval models.
While our analysis does not account for possible current-induced larval transportation, we cannot evaluate the proportion of larvae potentially not remaining within a region. However, since included SSB is at the ICES subdivision scale, movement of larvae by currents should not affect our conclusions on the SSB effect. Furthermore, if currents have moved larvae, we mitigate that covariates could have some mismatch with the actual spawning and early larval phase areas by the fact that we have used sampling region averages for the covariates and not covariate values at sampling sites, as described above. The risk that larvae from one sampling area are transported to another and thus are associated with different environmental covariate values can be neglected, since the sampling areas are usually >100 km apart, with complex islands mosaics separating them.

Phenological changes of Baltic herring larvae
We show that both developmental stages in Baltic herring larvae, small (<10 mm) and large (>15 mm), have shifted their timing to earlier dates (Figs. 4 & 5). First, our results suggest an earlier onset of small larvae. Depending on the probability threshold of detection, the decrease in the first days in occurrences of small larvae varied between 4.6 d per decade for very high probabilities of more than 90% and 7.9 d per decade for a still high but less strict threshold of more than 75% probability of occurrences. Second, our results show that the fractions of large larvae being higher than 10% of the total abundance has shifted by an average of 7.7 d per decade, varying between regions from 1.4 to 13 d per decade, and are associated with warmer conditions (Fig. 5). Our results are comparable to earlier reported magnitudes of phenological shifts in marine ecosystems (Genner et al. 2010, Poloczanska et al. 2013, Asch 2015. However, with an average earlier development of almost 8 d per decade, this advancement is faster than many observed shifts in phenological change in terrestrial systems (Menzel et al. 2006, Miller-Rushing & Primack 2008. Since the 4.6 d per decade advancement in the first detection of small larvae at a probability threshold of 90% is smaller than the 7.7 d per decade advancement for the large fraction of large larvae, our results indicate overall faster larval development. However, when considering the less strict 75% probability threshold for detecting the first small larvae (with a 7.9 d per decade shift), the overall change in phenology seems consistent and supports the robustness of our results and the magnitude of change.
Strong advancements in phenology can become critical for the survival of species when there is a mismatch in timing between prey and consumer (Edwards & Richardson 2004). Across realms, primary producers and primary consumers can advance seasonal timings at about twice the rate of secondary consumers (Thackeray et al. 2010 Fig. 6. Relationship of probabilities for large larvae >10% of total larval abundances as a function of realized temperature during the years 1974 to 1996. Regression lines are linear fits to area-specific data. Grey-shaded areas represent the 95% CIs of the linear regression lines. Color code represents the area colors in Fig. 1 year classes, i.e. lower SSB, via slower growth rates and starvation (Platt et al. 2003, Durant et al. 2007. Species that adapt to gradually changing environments by shifting their seasonal timings with sufficient pace to track their prey resources will likely not show strong negative effects of climate-induced phenology shifts (Cleland et al. 2012). In this light, strong shifts in phenology of secondary consumers such as Baltic herring larvae would likely be of benefit if they track the pace of their zooplankton prey resources and/or outpace potential predation pressure by higher consumers. However, more often it has been recorded that the diverging phenology of different trophic levels sets trajectories to mismatching resources and declines in consumers (Beaugrand et al. 2003, Mackas et al. 2007, Koeller et al. 2009). Despite the lack of adequate data to link changes in food resources for Baltic herring larvae at the spatio-temporal scale of this study, we know that the recruitment success of Baltic herring has decreased from the early 1980s to the early 2000s along with decreasing SSB in subdivisions 25 to 29 and 32 (ICES 2018). Our results that chl a is an important covariate in structuring the occurrence especially in the small and medium life stages of Baltic herring larvae (Fig. 3) highlights the importance of prey availability being a driver of year class strength (see also Rahikainen et al. 2017). Our analysis is based on the causal description for larval abundance (the Ricker model, Section S2 in Supplement 1), for which reason the conclusions of our results are more confirmatory than conclusions from a traditional approach that relies exclusively on correlations between the environmental covariates and SSB to explain changes in herring larval phenology. However, since our conclusions are only as good as our model, it is important to also consider alternative pathways that could have contributed to an earlier emerging fraction of large larvae which are not directly governed by the environment or do not affect adult stages of the Baltic herring. One such case could be size-specific predation pressure (e.g. Litvak & Leggett 1992) if altered over time. An increase in predation pressure on smaller larvae or a decrease in predation pressure on larger larvae could potentially contribute to shifting size class proportions. The sharp increase of sticklebacks in coastal areas of the Baltic Sea during recent decades (Berg ström et al. 2015, Eklöf et al. 2020 could indeed facilitate such changes in predation pressure during the early life stages of herring, despite seemingly no strong associations between stickleback abundance and white fish/vendace larvae (Vanhatalo et al. 2020). Yet, as we here in addition to proportional abundances of different-sized larvae also observed a matching shift in the onset of larvae development, i.e. first occurrences of small larvae per year, highlighting earlier dates of detection, the influence of predation pressure on our findings is likely not the main driver of the observed pattern. Another such case could be intensified fishing pressure, which is, for herring, most intense on pre-spawning aggregations in winter and spring in the Baltic. This could speculatively result in higher mortalities of late spawners and thus in higher proportions of early-spawned larvae. However, due to our coupled approach in modeling both earlier occurrences as well as earlier fractions of large larvae >10%, this assumption would also only influence the latter. While both ecological interaction as well as fishing pressure may contribute to the mechanisms of larval spawning and survival, our results provide strong support for environmental factors driving the shift in phenology.
The status of fish stocks affects the prospects of fisheries under climate change. With shifting timings of larval development of commercially important fish species, fisheries should consider an adaptive management approach, where stock assessment and management plans take the phenological advances of the species into account, to ultimately mitigate climate change effects on fisheries.

Conclusions
Our work highlights the shifting phenology of Baltic herring larvae in the northern Baltic Sea and provides a modeling framework to detect such changes based on the different life stages of larvae and early life stages of animals in general. We show that Baltic herring larvae have advanced their first occurrences as well as their growth-dependent abundance over a 22 yr time period consistently over a large spatial scale. Surface water temperature and chl a (serving as proxies for prey availability) were particularly strong drivers in explaining the occurrence patterns, while larval abundances showed strong density dependence effects explained by the specific SSB. With increasing temperatures having positive effects on occurrence and abundance, indicating high probabilities of earlier development and detection of larvae, a continued advancement of Baltic herring larval phenology can be expected under future climate projections. Long-term data series, such as the one utilized in this study, need to be maintained and broadened to also include lower trophic levels, as this will be fundamental to under-standing how the re sulting predator−prey dynamics form altered species phenology in the Baltic Sea and elsewhere.