Land use change and coastal water darkening drive synchronous dynamics in phytoplankton and fish phenology on centennial timescales

At high latitudes, the suitable window for timing reproductive events is particularly narrow, promoting tight synchrony between trophic levels. Climate change may disrupt this synchrony due to diverging responses to temperature between, for example, the early life stages of higher trophic levels and their food resources. Evidence for this is equivocal, and the role of compensatory mechanisms is poorly understood. Here, we show how a combination of ocean warming and coastal water darkening drive long‐term changes in phytoplankton spring bloom timing in Lofoten Norway, and how spawning time of Northeast Arctic cod responds in synchrony. Spring bloom timing was derived from hydrographical observations dating back to 1936, while cod spawning time was estimated from weekly fisheries catch and roe landing data since 1877. Our results suggest that land use change and freshwater run‐off causing coastal water darkening has gradually delayed the spring bloom up to the late 1980s after which ocean warming has caused it to advance. The cod appear to track phytoplankton dynamics by timing gonadal development and spawning to maximize overlap between offspring hatch date and predicted resource availability. This finding emphasises the importance of land–ocean coupling for coastal ecosystem functioning, and the potential for fish to adapt through phenotypic plasticity.


| INTRODUC TI ON
At high latitudes, phenological events tend to be tightly synchronized across trophic levels (Durant et al., 2007;Stenseth et al., 2002;Sundby et al., 2016).This is due to large seasonal differences, which in marine ecosystems cause a pulsed planktonic production in spring (Cushing, 1990;Hjort, 1914;Hughes, 2000).During winter, nutrient concentrations are high, but phytoplankton growth and concentration is constrained by short days and a deep mixed layer, constantly moving phytoplankton out of the shallow euphotic zone.Towards spring, as day length and temperature increase and the water column stabilize, phytoplankton is retained in the euphotic zone and the spring bloom starts (Behrenfeld & Boss, 2014;Doney, 2006;Lindemann & St. John, 2014;Sverdrup, 1953).Around the same time, zooplankton ascend from their overwintering depths to graze on phytoplankton and initiate spawning (Broms & Melle, 2007).In turn, zooplankton eggs and nauplii are an important food resource for newly hatched fish larvae (Beaugrand et al., 2003;Cushing, 1990).
The increased turbidity associated with the phytoplankton bloom also provides fish larvae with some protection from visual predators without impairing their own foraging ability (Fiksen et al., 2002).This synchronized propagation from physics to fish is well-established (Hjort, 1914) and is often inferred when linking variation in bloom timing and fish spawning to recruitment success in fish (Malick et al., 2015;Platt et al., 2003;Schweigert et al., 2013).
While the theory of phenological match-mismatch was initially formulated for marine ecology to explain how climatic and environmental variability drive interannual fluctuations in fish stocks (Cushing, 1974(Cushing, , 1990;;Hjort, 1914), the concept of phenological asynchrony driven by long-term climate change stems largely from terrestrial ecology.This background is likely due to the comparative ease of collecting phenological data in terrestrial systems (Samplonius et al., 2021), eliciting long timeseries on phenological events such as plant flowering time, insect bursts and bird egg laying (Both et al., 2008;Thackeray et al., 2016;Walther et al., 2002).This contrasts that of marine ecology, where phenological timeseries are generally scattered and fragmented, and inference of climate change on trophic synchrony has largely been pieced together with information from different ecosystems and time periods (Asch et al., 2019;Cooley et al., 2022;Poloczanska et al., 2013Poloczanska et al., , 2016)).
In marine ecology, the most common explanation for the occurrence of climate driven trophic asynchrony is the discrepancy between the direct temperature-effect on fish physiology and hence vitellogenic rates and spawning phenology (Kjesbu et al., 2010;Neuheimer & MacKenzie, 2014;Pankhurst & Munday, 2011), and the indirect temperature-effect on phytoplankton phenology mainly though environmental drivers such as stratification (Asch et al., 2019;Doney, 2006).Asch et al. (2019) estimated that increasing temperatures will shift fish spawning time at a rate twice that of the phytoplankton bloom onset, and that ocean warming will lead to a gradual decoupling between the two.A recent IPCC meta-analysis (Cooley et al., 2022) came to an opposite conclusion and found that observed phenology change was more than two times faster for phytoplankton (−7.8 days decade −1 ) compared to fish (−3.2 days decade −1 ).The same study also found that changing phenology in lower trophic levels such as phytoplankton and zooplankton were more consistent (72-81%) with climate change expectations than that of fish and seabirds (42-65%).Another meta-analysis of 51 observations of spring phenology across five functional groups (phytoplankton, benthic invertebrates, zooplankton, fish and seabirds) found no clear link between phenological change and that expected from changing temperatures (Poloczanska et al., 2013).
An additional, but largely unexplored effect of climate change, relates to the land-ocean coupling.Several studies have suggested that centennial land use change, reduced sulphur emissions and a warmer climate in boreal areas has caused an increase in coloured dissolved organic matter (CDOM) in freshwater lakes and rivers, known as browning (de Wit et al., 2016;Kritzberg, 2017;Monteith et al., 2007).In addition, hydrology is a key driver of CDOM-export from catchments (de Wit et al., 2016), and a long-term increase in annual precipitation in Northern Europe since the 1950s (Ballinger et al., 2023;Hurrell, 1995) has caused an overall increase in freshwater run-off to the coast (Saetre et al., 2003).In turn, the CDOM content in this freshwater causes reduced light penetration in coastal waters (Aksnes et al., 2009;Opdal et al., 2023).
A general challenge for marine studies is that phenology changes at different trophic levels are derived from separate ecosystems and/or time periods, and with large variations between species and regions (Poloczanska et al., 2016).While several studies have detected significant correlations between fish phenology and temperature (McQueen & Marshall, 2017;Rogers & Dougherty, 2019;Wieland et al., 2000), a recent review of 109 papers covering 129 taxa (Samplonius et al., 2021) found that only a few studies could document that climate warming led to trophic decoupling (Edwards & Richardson, 2004;Philippart et al., 2003).Simultaneous long-term observations of ocean temperature, timing of resource availability and higher trophic level phenology in a single region or ecosystem are rare (Samplonius et al., 2021).
Apart from the lack of good data, one explanation for this apparent discrepancy between hypothesised climate driven trophic asynchrony and observations in the field could be the effect of adaptive phenotypic plasticity, which may mitigate and counteract the expected trophic mismatch (Charmantier et al., 2008;Merilä & Hendry, 2014).For fish, several studies have shown that adaptive phenotypic plasticity is a common response to increasing temperatures (Crozier & Hutchings, 2014), but has to our knowledge not been found in direct response to phenological shifts at lower trophic levels.In a study comparing 21 different populations of Atlantic cod, Neuheimer and MacKenzie (2014) found that cod accurately match their spawning time to the variability in the local spring bloom onset, and hence offspring food availability, despite a wide temperature range (ca.10°C).Similarly, Opdal, Wright, et al. (2024) found that the Northeast Arctic cod appeared to spawn in tight synchrony with the spring phytoplankton bloom, largely independent of temperature variations.This could suggest that some sort of environmental cues rather than simple temperature-driven rates may be important to tune spawning behaviour to match local environmental conditions (Neuheimer et al., 2018;Opdal, Wright, et al., 2024).A similar local adaptation to match egg hatching time with the spring bloom appear to have evolved in shrimp populations (Koeller et al., 2009).
Thus, to what degree adaptive phenotypic plasticity in fish matches long-term changes in the phenology of lower trophic levels remains largely unresolved.
In this study, we use a weekly resolved biological dataset for spawning time in Northeast Arctic (NEA) cod dating back to 1877, alongside bi-weekly hydrographical depth profiles since 1936, all at the same location in Lofoten, Norway.We also make use of a previously published timeseries of non-phytoplankton light attenuation for the freshwater endmember of the Norwegian Coastal Water (NCW, Opdal et al., 2023), that is, a background attenuation that account for light attenuation due to CDOM but not phytoplankton.This quite unique combination of high-resolution and long-term timeseries allow us to (1) establish a 145-year timeseries of spawning phenology in cod, (2) create a hindcast of the phytoplankton spring bloom timing and (3) disentangle the effect of biological and physical drivers on the phenological dynamics of the highly ecologically and economically important northeast Arctic cod population.

| MATERIAL S AND ME THODS
All below-mentioned variables, including abbreviations, units and values are listed in Table 1.

| The phytoplankton spring bloom
In Lofoten, the main spawning area for NEA cod (Figure 1), we estimated the phytoplankton spring bloom onset using two approaches.
For the years 1998-2021, the bloom onset was derived from remote sensing, while for the longer time-period , we used a hydrographical approach based on long-term observations at the coastal monitoring station Skrova (68.12 N, 14.65 E) situated in Lofoten (Figure 1).The latter approach was a first-order approximation, where bloom onset was defined as the day when the critical depth (Z CR , m) becomes deeper than the mixed layer depth (MLD, m;Sverdrup, 1953).

| Remote sensing estimate
For the years 1998-2021, the spring bloom onset in Lofoten was derived from satellite images of surface chlorophyll a concentration estimates (mg m −3 ) available through the Copernicus Marine Service (European Union).We utilized the chlorophyll-a concentration product (CHL) from the Global Ocean Colour level-3 multi sensor product with daily gap-free estimates at a spatial grid resolution of 4 × 4 km (https:// doi.org/ 10. 48670/ moi-00280 ).
This product is known to perform well in comparison with in situ measurements (Garnesson & Bretagnon, 2022) and has previously been used successfully in Norwegian coastal waters (Opdal, Wright, et al., 2024).
For comparison and calibration with the hydrographical bloom estimate (see below), the spring bloom onset was calculated for 25 (5 × 5) grid-cells located around the Skrova monitoring station.The spring bloom onset was initially defined (B R1 ) as the first day-of-year when the chlorophyll a concentration exceeded the annual median with 5% or more over 3 or more consecutive days (Henson et al., 2009).In addition, we tested an alternative definition (B R2 ) (Brody et al., 2013), corresponding to the first day-of-year where the chlorophyll a concentration doubled the next day.Due to little light in winter, satellite images were available at the earliest on February 23th (Day 54).To avoid a confounding effect between estimated bloom timing and image availability, estimates that occurred at the same day as the first available image were discarded.
The critical depth (Z CR ) is defined as the depth, where the vertically integrated light-limited gross primary production equals vertically integrated phytoplankton losses (Marra, 2004;Nelson & Smith, 1991;Sverdrup, 1953), and can be described by where K (m −1 ) is the non-phytoplankton attenuation coefficient of downwelling irradiance, I 0 (mol photons m −2 day −1 ) is the daily integrated photosynthetically active radiation (PAR, 400-700 nm) just below the surface, and I c (mol photons m −2 day −1 ) is the irradiance at the (photo-) compensation depth.
The non-phytoplankton light attenuation, K, of Norwegian Coastal Water (NCW), where NCW is defined as the water masses with salinity <34.5 psu (Figure 1, Saetre, 2007), reflects the amount of CDOM of freshwater origin and can be approximated according to a conservative mixing model with freshwater (FW, zero salinity) and North Atlantic Water (NAW, salinity 35.2 psu) as end members (Aksnes, 2015;Opdal et al., 2023): where g = S/35.2and S is the NCW salinity that defines the mixture between FW and NAW.K NAW is the non-phytoplankton light attenuation of the NAW and is set to 0.03 m −1 (Aksnes, 2015).This value was estimated for the wavelength of 440 nm and corresponds to oceanic water type IA/B (Jerlov, 1978).K FW is the non-phytoplankton light attenuation of the freshwater endmember, which has been gradually increasing throughout the 20th century (Dupont & Aksnes, 2013;Kahru et al., 2022;Opdal et al., 2019Opdal et al., , 2023)).Here, we use the annual estimates of K FW   To get the daily integrated PAR just below the surface (I 0 , μmol photons m −2 day −1 ), we started by obtaining hourly values of (downward) shortwave sea surface radiation (I surf , W m −2 ) from the ERA-20CM product (Hersbach et al., 2015), provided by the European Center for Medium-range Weather forecast (ECMWF), which in ice-free regions are found to match airborne observations (Müller et al., 2023).For each day with hydrographic observations in the period 1936-2021, I 0 was calculated according to where 86,400 is the number of seconds in a day, cf w = 4.57 is a conversion factor from W m −2 to μmol photons m −2 s −1 (Thimijan & Heins, 1983), cf surf = 0.6, accounts for light reduction at the surface due to albedo and absorption, and cf PAR is the fraction of I surf that is PAR and is set to 0.47, which is the mean of published correlations between the PAR photon flux density and daily broadband solar irradiance (Jacovides et al., 2004).The total range of values reported by Jacovides et al., 2004 is between 0.39 and 0.54.Finally, the irradiance at the photo-compensation depth (I c ) was set to 1.0 mol photons m −2 day −1 , as estimated by Siegel et al. (2002) for the North Atlantic at latitudes between 65 and 70° N.
The MLD was estimated from the density stratification derived from temperature and salinity (TS) bi-weekly observations at the Skrova monitoring station (Figure 1).We defined the MLD as the first depth where the density exceeds the density at 10 m depth by a given difference (e.g., de Boyer Montegut et al., 2004)-referred to here as the density difference criterion (DC).The DC was calculated as the average density difference (across years) between 10 m depth and the critical depth (Z CR ) at the day of the bloom onset as derived from remote sensing. (1)

| Spawning time for Northeast Arctic cod
The spawning time for NEA cod in the Lofoten area was calculated based on the relationship between weekly resolved commercial landings of NEA cod and NEA cod-roe between 1877 and 2022.These data were drawn from the Official Fisheries Statistics of Norway, available as scanned books and reports for the period 1877 to 2014 (Anon, 1877(Anon, -2022)), and in digital format for the years 1990-2022 (Hopland & Aasheim, 2023).The complete data set is available for download in Opdal, Lindemann, et al. (2024).A subset of these data  has previously been published by Pedersen (1984), who showed that the time of spawning, estimated from egg sampling, was strongly correlated to the day-of-year (t) with the highest average female gonadosomatic index (GSI).The weekly (w) mean female GSI in any given year (yr) was defined as the total weight of landed roe (R(w,yr), kg) relative to the total weight of landed fish (W (w, yr), kg) multiplied by the male-female ratio, assumed to be 0.5.
Because data are weekly resolved, a piecewise polynomial spline interpolation was fitted to the roe (R, kg) and landing (W, kg) data so that for each year, we obtain a daily (t) mean female GSI estimate (GSI(yr, t)) used to find the day-of-year of peak GSI (t GSImax , Figure S1) For the years 1980 to 2019, the t GSImax estimation was compared with an independent estimate of the spawning time of NA cod (dayof-year when >50% are spawning) based on the gonad development stage of >126.000sampled individuals along the Norwegian coast, including the Lofoten area (Opdal, Wright, et al., 2024).

| Surface temperature in Lofoten, 1868-2022
While the bi-weekly TS-profiles at Skrova were used to calculate the  .To check for consistency between the two timeseries, they were compared using the linear model T S = α T A + β, fitted using reduced major axis regression (Sokal & Rohlf, 2012).

| Statistical analyses
All timeseries datapoints are presented as means with error bars denoting the standard error of the mean.Datapoints without error bars are thus single values.For graphical purposes, we also calculate and show smoothed splines with 95% confidence intervals [CI].For this, we use a locally weighted scatter plot smooth with linear interpolations and span (f) set to 0.7 (Cleveland, 1979).
Temporal trends were analysed for all timeseries using simple linear regressions.Due to their segmented form, possible changepoints were also investigated.We used the findchangepts() function in Matlab R2023a (Killick et al., 2012;Lavielle, 2005), and changepoints defined as a significant change in mean and slope between 1936 and 2022 (statistic = 'linear').

| RE SULTS
The results of the hydrographical estimation of the phytoplankton spring bloom onset, which is based on Sverdrup's critical depth hypothesis (Equation 1), is illustrated in Figure 2. We note that during late winter (Day 50) incoming PAR just below the surface (I 0 ) was lower than 7 mol photons m −2 day −1 with an average of 3.8 mol photons m −2 day −1 .Towards mid-spring (Day 110), this increases to between 5 and 23 mol photons m −2 day −1 , with an average of 16 mol photons m −2 day −1 (Figure 2a).
In late winter, the non-phytoplankton light attenuation (K) varies from 0.044 to 0.15 m −1 , depending on year.Thus, the combination of I 0 and K results in critical depths (Z CR ) ranging between 0 and 40 m in late winter (Figure 2b).As such, observed variation in solar irradiance at this time will shift Z CR by 40 to 140 meters depending on K, while observed variation in K will shift Z CR by up to 150 meters depending on I 0 (Figure S2).As I 0 increases towards the spring, the Z CR becomes gradually deeper (Figure 2b).(4) GSI(w, yr) = R(w, yr) W(w, yr) • 0.5 The sea surface temperature anomaly in the Lofoten area in spring (1867-2022, Figure 3a) was estimated by combining temperature recordings from the hydrographical station Skrova  and the Andenes lighthouse (1867-1971; Figure 1 and Figure S5).The changepoint analysis suggest a shift in temperature trend in 1988, after which surface temperature starts to increase by an average of 0.24°C decade −1 (Table 2).Prior to this, temperature has remained relatively stable, although a slight linear increase of 0.048 [0.013 0.82] °C decade −1 is evident when analysing the entire period from 1867 to 1988 (N = 111).For salinity (Figure 3b), shown as the depth-integrated mean anomaly (0-30 m, March-May), the changepoint analysis detects a shift in mean and linear slope in 1988 (Table 2).While the linear slopes before and after 1988 are not statistically different from zero, the average salinity after 1988 is statistically lower than before 1988.Considering the entire period (1936-2022, N = 86), an overall decline of −0.055 [−0.079-0.030]psu decade −1 is observed.For the seasonal formation of the mixed layer (Figure 3c and Figure S4), which depends on both temperature and salinity, there is a shift towards earlier formation of the mixed layer after 1978 (Table 2).Before 1978, we see no statistically significant trend, and on average, the MLD is 30 m around Day 152.
After 1978, the time the MLD reaches 30 m advances by 11 days per decade, and occurs on average 2 weeks earlier, on Day 135.
The combination of increasing light attenuation of the freshwater endmember (K FW ) of the NCW (Opdal et al., 2023) and decreasing salinity (Figure 3b) drives an overall increase in non-phytoplankton light attenuation (K) in Lofoten (Figure 3d).The changepoint analysis suggests an increase in K of 0.004 m −1 decade 1 up to 1988, before gradually levelling out (Table 2).This causes a delay, of 3.0 days decade −1 between 1936 and the changepoint in 1988, in the estimated hydrographical spring bloom onset.After 1988, there is a stabilisation.
For the estimated spawning time of NEA cod (Figure 2e), which was based on the time of peak GSI (Figure S1), we found an overall delay between 1877 and 2022 of 2.1 [1.6-2.5]days decade −1 (N = 143).A changepoint was detected in 1994, after which there is a levelling off (Table 2).Between 1936 and1994, cod delayed spawning time by 4.1 days decade −1 .For the years 1980-2019, we could compare the timing of peak GSI (t GSI max ) with the time when 50% of individuals have gonads in spawning stage (t P50 ).A ranged major axis regression, t GSImax = a• t P50 + β suggests reasonable coherence between the two; R 2 = 0.42 and a = 1.41 [0.9-2.4],N = 34 (Figure S7).Here, we also note that over the period 1980-2019, both estimates indicate that cod spawned gradually earlier.A linear regression analysis gives an advancement of −6.7 [9.1-4.4]days decade −1 for the timing of peak GSI (t GSImax ), and − 7.   1, thin yellow lines).Thick blue and yellow lines with shading indicate the 10th and 90th percentiles of the MLD and Z CR , respectively.Circles denote the day of year where Z CR ≥ MLD, at which point, a light-limited phytoplankton growth rate will be equal to the loss rate.This point is the estimated the phytoplankton bloom onset.For illustrative purposes, circles denoting bloom onset for the first (1936)(1937)(1938)(1939)(1940)(1941)(1942)(1943)(1944)(1945) and last (2012-2021) decades have been coloured pink and purple, respectively.Panel (c) shows the day of year (colour scale) where Z CR ≥ MLD relative to the mean mixed layer depth in winter (Jan-Feb) and the corresponding mean non-phytoplankton light attenuation (K).Dots denote the same estimates as in panel (b).

| DISCUSS ION
In the coastal waters of Lofoten, the main spawning area for the NEA cod, our observations provided detectable changepoints between 1978 and 1994.These include the three independently observed variables; temperature, salinity and cod spawning (1988)(1989)(1990)(1991)(1992)(1993)(1994), as well as the three estimated variables; light attenuation, shallow mixed layer formation and phytoplankton spring bloom onset (1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988).For the years prior to the changepoints, we found that a combination of freshwater browning and coastal water freshening has gradually increased the light attenuation, which in turn drive long-term delay in the phytoplankton spring bloom onset.The NEA cod, in turn, appear to gradually adjust their spawning time so that synchrony with the changing phytoplankton bloom timing is maintained.In the years after the changepoints, light attenuation, the timing of the phytoplankton bloom, and the spawning time of NEA cod stabilize.As such, it appears that atmospheric, terrestrial and hydrographical factors (see conceptual Figure 4) have similar effect on the timing of the spring phytoplankton bloom as well as on the timing of the cod spawning.
While water column density stratification (as derived from the observed temperature and salinity), determines the mixed layer depth in our study, the incoming solar irradiance (which also affects temperature) and the underwater non-phytoplankton light attenuation define the critical depth.Together, these factors drive the variation in the phytoplankton spring bloom onset in this study (Figure 4).The initial delay  in phytoplankton spring bloom onset occurred during a period of quite stable temperatures but increasing light attenuation.The latter essentially shoals the critical depth, thereby requiring a relative shallower mixed layer depth to initiate the spring bloom.A delay in bloom timing associated with elevated light attenuation, is mechanistically wellsupported through general theory of phytoplankton growth-loss dynamics (Opdal et al., 2019;Urtizberea et al., 2013).For the period 1988-2022, the mean temperature increased by ca 1°C, and a clear advancement is seen in the formation of a shallow mixed layer.After 1988 there is a tendency towards earlier bloom timing (as well as cod spawning) as would be expected from, solely (i.e., unaltered critical depth), an advancement in the formation of a shallow mixed layer.
For the period 1998-2021, our hydrographical estimates of the onset of the phytoplankton spring bloom are largely in agreement with the estimates from remote sensing, giving confidence to the methodology.Despite questions concerning the universality of the critical depth concept (e.g., for the Pacific Ocean, Behrenfeld, 2010) originally suggested by Gran and Braarud (1935) and later formalised by Sverdrup (1953), our results suggests that it can be a useful proxy for the spring bloom onset in the NCW.
The effect of freshwater browning on light attenuation (i.e., increased K FW and associated darkening of the water column), however, cannot be validated using remote sensing products since the most prominent change in K FW occurred before these products existed.
We note that the long-term average surface temperature in Lofoten has been stable up to the late 1980s, after which it gradually increases.A similar pattern of gradual warming after the 1980s has also been shown for the North Atlantic in general (Wang & Dong, 2010).The long-term decline in salinity in Lofoten reflects  2).
the general freshening of the NCW, as also seen at the hydrographical station in Sognesjøen (61.02°N, 4. 84° E, 1936-2005) (Aksnes et al., 2009), as well as for the entire southern Norwegian coastal waters (<61° N, 1903-2021) (Opdal et al., 2023).This long-term coastal water freshening, and an associated increase in the NCW non-phytoplankton light attenuation (K), has been termed 'coastal water darkening' (Aksnes et al., 2009) thought to be driven by a long-term increase in precipitation in Northern Europe (Ballinger et al., 2023;Hurrell, 1995) and a subsequent increase in freshwater run-off to the NCW (Saetre et al., 2003).In addition to this freshening effect, additional NCW darkening is caused by a centennial increase in the light attenuation of the freshwater itself (i.e., increased K FW ), primarily due to land-use change and an expansion of the forest cover around the Baltic Sea (Opdal et al., 2023).Since freshwater and its light attenuation, K FW , is endmember in our mixing model for the NCW (Equation 2), a temporal increase in K FW causes, according to Equation 2, increased non-phytoplankton light attenuation, K (Figure 4).
With future ocean warming and freshening, we expect a continued advancement in the formation of a shallow mixed layer, and thus an advancement of the spring phytoplankton bloom (Doney, 2006).
On the other hand, further increase in rainfall in Northern Europe (Ballinger et al., 2023) and associated freshening of the Norwegian Coastal Water may also increase K, possibly exacerbated by terrestrial greening that further increases K FW (Crapart et al., 2023;Kritzberg et al., 2020;Opdal et al., 2023).In turn, this may push bloom timing in the opposite direction of the temperature effect (Opdal et al., 2019).Ambient temperature is expected to dictate physiological rates in ectotherms (e.g., Dillon et al., 2010;Gillooly et al., 2001), with direct effects on phenology (Kjesbu et al., 2010;Parmesan, 2007) although some studies question such a simplified relationship between temperature and metabolic rates (Clarke, 2004;Schulte et al., 2011).A TA B L E 2 Changepoint analysis for all timeseries variables between 1936 and 2022.common mechanism to explain variation in spawning time is variation in the vitellogenic temperature, i.e., the temperature during oocyte development in autumn and winter.High or low temperatures will reduce or prolong oocyte development time, and in turn advance or delay spawning time, respectively (Kjesbu et al., 2010;Neuheimer & MacKenzie, 2014).This mechanism has also been proposed to de-couple the spring bloom and the fish spawning time, because temperature is expected to change the spring bloom timing (indirectly through hydrography) at a slower rate than the direct effect of temperature on oocyte development time (Asch et al., 2019).
In our study, we do not detect such a proposed decoupling.The gradual delay in spawning time between 1877 and 1994 of around 40 days, cannot be explained by temperature driven vitellogenic rates.This would have required declining temperatures during this period, which was not the case.Rather, for the period where estimates of phytoplankton bloom onset are available , the rates of change in spring bloom onset and in cod spawning are remarkably similar up to the 1988 and 1994 changepoints, respectively.After this, there was a significant increase in temperature, but no significant trend towards earlier bloom or spawning.However, considering the period 1980 to 2021, where we had two independent estimates of cod spawning time, a clear trend towards earlier spawning was evident, and could thus be related to increasing temperatures.Even though temperature does influence physiological rates and gonadal development in NEA cod (Kjesbu et al., 2010), a recent study analysing both experimental and spatiotemporal variation in NEA cod spawning found that cod could advance or delay spawning by several weeks, largely independent of ambient temperatures (Opdal, Wright, et al., 2024).The latter corroborates the findings in our study, where cod spawning appears to synchronise with the phytoplankton bloom onset.
In another study, based on a subset of the NEA cod roe data , Pedersen (1984) found a two-weeks delay in spawning time, and hypothesised that this could be caused by a concurrent 3-year decrease in the mean age of the spawning stock.The reasoning was that, historically, older fish appear at the spawning grounds earlier in the season than the younger ones (e.g., Rollefsen, 1938Rollefsen, , 1939;;Sund, 1938Sund, , 1939)), and that a truncation in the spawner demography towards younger individuals would inevitably delay the time of peak spawning.Earlier spawning for older cod has also been found in the Baltic Sea (Wieland et al., 2000), while they were found to spawn later in the season in the Northwest Atlantic (Hutchings & Myers, 1993) These differences, however, may reflect methodological differences in defining spawning time (Morgan et al., 2013).
Unfortunately, apart from some modelling estimates dating back to 1913 (Hylen, 2002) reliable data on the mean age of the NEA cod spawning stock are unavailable prior to the 1930s (Jørgensen, 1990;Opdal & Jørgensen, 2015).These studies all suggest that the spawner age started to decline around the 1950s, during the time with the highest recorded catches from the Barents Sea trawl fishery (Godø, 2003).Hence, the high fishing pressure has likely reduced the mean age of spawners, directly through demographic effects (Jørgensen, 1990), but likely also through fishing induced evolution (Heino et al., 2002;Jørgensen et al., 2007).While our findings do not contradict the hypothesis of Pedersen (1984), we note that NEA cod spawning time has occurred gradually later since the 1880s, more than four decades before the onset of the Barents Sea fishery in the 1920s (Godø, 2003).As such, changes in spawning time prior to this will likely have other drivers than the age of the spawning stock.
The reproductive success of the NEA cod has been suggested (Ellertsen et al., 1989;Sundby, 2000) to be highly dependent on spawning time in relation to prey abundance for their first-feeding larvae, typically zooplankton eggs and nauplii (Ellertsen et al., 1989;Sundby, 2000).Similar relationships between fish recruitment success and match between spawning and the phytoplankton bloom are found for several other systems (Malick et al., 2015;Platt et al., 2003;Schweigert et al., 2013).In the NCW, abundance and egg production in zooplankton (Calanus finmarchicus) are strongly coupled to the phytoplankton spring bloom (Broms & Melle, 2007).
Hence, if cod detect the presence-absence of an early bloom signal or associated hydrographical characteristics they may delay or advance spawning to maximize first feeding success and survival of their offspring.Since temperature influences the spring bloom and vitellogenic maturation differently (Asch et al., 2019), and offspring survival relies primarily on the spring bloom dynamics, we should expect a strong selection towards a behavioural plasticity that allows individuals to de-couple or offset spawning time from interannual fluctuations in temperature (Opdal, Wright, et al., 2024).Local adaptations to match egg hatching time with spring bloom timing seem to have evolved in different populations of shrimp (Koeller et al., 2009) and Atlantic cod (Neuheimer et al., 2018).Our analysis strengthens the hypothesis that such plasticity may not only exist between populations, but also within a single population of NEA cod (Opdal, Wright, et al., 2024).
Our results suggest a spawning plasticity that agrees with laboratory experiments (Kjesbu et al., 2010) and field observations (Neuheimer et al., 2018;Opdal, Wright, et al., 2024) that have un- Individual plasticity aside, our study reveals a thought-provoking link between atmospheric and terrestrial processes in Scandinavia on one hand and the spawning time of a large oceanic fish stock on the other.A warmer climate and increased precipitation combined with less pasture and grazing have led to more forests and therefore a higher run-off of freshwater as well as higher concentrations of coloured dissolved organic matter in rivers, lakes and coastal waterswhich, in turn, affect the timing of the spring bloom (Figure 4).Our results suggest that cod have adjusted to these long-term changes, likely through environmental cues enabling them to control their spawning intensity.
presented inOpdal et al. (2023) which they found valid for the NCW, including the monitoring station Skrova in the Lofoten area.This means that our estimates of the nonphytoplankton light attenuation at Skrova vary according to the variation in the measured salinity as well as with the long-term change in K FW .Given the low phytoplankton concentration prior to the bloom onset our approach does not consider phytoplankton-driven attenuation as a driver of bloom onset variability.

1
Geographical overview of the Northeast Arctic cod feeding and spawning areas.The Northeast Arctic cod has its main feeding and nursery areas in the Barents Sea.In late autumn, mature cod migrate south to spawning grounds along the Norwegian coast.Spawning takes place between January and May.Eggs and later hatched larvae are carried by the northbound Norwegian Coastal Current (NCC) back to the Barents Sea.The insert shows the main spawning grounds (hatched) around the Lofoten archipelago based onBrander (1994).Red colour denote the Norwegian Coastal Water (<34.5 psu), based on the depthintegrated salinity in the upper 30 m between January and March in the years 2015-2022(Lien et al., 2013).
Figure1).Here, surface temperature has been regularly recorded between 1868 and 1963 by the keeper of Andenes lighthouse and made available as monthly means byFrogner (1948) for the years1868-1945 and by the Norwegian Meteorological Institute (1946-  1963).The Andenes data sets are available for download inOpdal,   Lindemann, et al. (2024).A complete record of surface temperature in Lofoten (1868-2022) was achieved by combining the Andenes timeseries (T A , with the temperature at 1 m depth at the Skrova station (T S , 1936-2022) by their monthly anomalies in spring (March, April and May) relative to the mean temperature in the overlapping years.To check for consistency be-

For
the MLD, the density criterion (DC) was estimated to 0.32 [0.24 0.41] kg m −3 (N = 23) based on the B R1 bloom estimate, and to 0.27 [0.21 0.32] kg m −3 (N = 21) based on the B R2 bloom estimate.See Figure S3 for annual plots of daily depth-resolved densities and DC-depth.Below, all estimates are based on DC = 0.32 kg m −3 , unless otherwise mentioned.From Figure 2b, we see that the MLD becomes shallower over the season, but with large interannual variations (Figure S4).During late winter, a variation in MLD of more than 50 m is observed.The estimated phytoplankton bloom onset (Equation 1) occurs between late February (Day 51) and late March (Day 86), depending on year.The respective effects of MLD and K on phytoplankton bloom onset day are shown in Figure 2c.Overall, the earliest bloom onsets are characterized by shallow MLD and low K, while the latest bloom onsets generally occur at deeper MLD and at higher K.

F
I G U R E 2 Estimation of the phytoplankton spring bloom onset in Lofoten for the years 1936-2021.Panel (a) shows the solar irradiance just below the surface (I 0 ) between mid-February to early May.Each line denotes a single year.Panel (b) shows, for each year, the seasonal development of the mixed layer depth (MLD, thin blue lines) and the critical depth (Z CR , Equation

F
I G U R E 3 Hydrographical and phenological timeseries and changepoints in Lofoten, Norway, 1868-2022.Panel (a) estimated annual mean spring surface temperature based on measurements taken at the Andenes lighthouse (squares, 1867-1935) and at the hydrographical station Skrova (circles, 1936-2022).Panel (b) the annual mean spring salinity anomaly between 0 and 30 m. Panel (c) the day-of-year anomaly when the mixed layer (ML) reaches 30 m. Panel (d) mean non-phytoplankton light attenuation between 0 and 30 m. Panel (e) the estimated onset of the phytoplankton spring bloom, shown as day-of-year anomaly.Panel (f) the estimated spawning time anomaly for Northeast Arctic cod, based on the seasonal development of gonad size.For panels a, b, d and e error bars denote the standard error of the mean.For all panels, lines and shade show a fitted locally weighted scatter plot smooth and its 95% CI.Vertical lines mark the estimated trend changepoints (Table Note: Changepoints (CP) were determined using the findchangepts() function in Matlab R2023a evaluating change in both mean and linear trend.See methods for further details.Based on the estimated changepoint, the statistical parameters (Stat.)mean (x) and change per decade (α) are calculated for each variable before (Pre CP, 1936-CP) and after (Post CP, CP-2022) the estimated changepoint.Change per decade (α) is estimated based on a simple linear regression, variable value = β + α•year.Statistically significant differences in means (before and after change point) and linear trends are denoted in bold.F I G U R E 4 A conceptual sketch of the most important atmospheric, terrestrial and hydrographic drivers and processes behind phytoplankton bloom dynamics and fish spawning time.Boxes with black frames and black whole-line arrows show drivers and processes explicitly addressed in this study.These boxes have colours referring to the same colour-schemes used in Figures 2 and 3. Boxes with grey frames and arrows are known drivers and processes but are not explicitly addressed in this study.Dashed black arrows indicate potential direct physiological effects on phytoplankton growth rate and fish oocyte development rate.I surf = surface irradiance (mol photons m −2 day −1 ), K FW = non-phytoplankton light attenuation coefficient of the freshwater endmember (m −1 ), Z CR = Sverdrup's critical depth (m), K = non-phytoplankton light attenuation in NCW (m −1 ), MLD = mixed layer depth (m).
covered a high degree of plasticity in the spawning time of Atlantic cod.This plasticity is partly independent of temperature, but tightly coupled to the availability of food for the offspring.Still, the underlying mechanisms for this plasticity remain unresolved.What cues are being utilized, and what physiological processes are at play in delaying or advancing spawning time?To make reliable predictions of fish phenology under climate change, we should know the physiological mechanisms that drives spawning time variability, and which may outweigh temperature-driven enzyme kinetics.