Long‐term shifts in phenology, thermal niche, population size, and their interactions in marine pelagic copepods

Under climatic warming many species shift their seasonal timing of life cycle events (phenology) and seasonal abundance distribution, but whether they maintain the same thermal niche is still poorly understood. Here, we studied multidecadal trends in abundance and phenology of seven major copepod species across three stations (Stonehaven (SH), Helgoland Roads (HR), and Plymouth L4) on the North–West European shelf, spanning ~ 6.5° of latitude. All seven species consistently occupied colder temperatures at the northern station compared to the southerly station, but they maintained the same realized thermal niche over years. Expected phenological shifts (i.e., earlier when warmer) in some stations were obscured possibly by the long‐term drop of copepod density in spring–summer, which may be due to a variation in the food/predators abundance. The ongoing spring–summer declines in abundance (~ 50%) of many North Atlantic pelagic species over the last five decades, as found in recent studies, may have also influenced the metrics of seasonal timing. To separate the seasonal timing of life events from that of seasonal abundance distribution, we used a time series of egg production rate (EPR) of Calanus helgolandicus at L4, and found that this shifted later into the summer–autumn over the last 30 yr of warming, coincident with declining spring–summer food and increasing predator abundance. Overall, direct temperature effects do appear to influence the seasonal timing of the copepods, but to explain impacts at individual stations or long‐term trends in population size or phenology, understanding the changing balance of food and predators appears to be critical.

Climatic warming is affecting many aquatic and terrestrial species and has led to varying degrees of ecological disruption (Walther et al. 2002;Burrows et al. 2011;Poloczanska et al. 2013).Three "universal" responses to climate change have been described in relation to planktonic species: body size response to warming (Atkinson 1994, Daufresne et al. 2009;Horne et al. 2017), phenological response (Parmesan and Yohe 2003;Edwards and Richardson 2004), and geographical redistribution (Beaugrand et al. 2002;Chivers et al. 2017).A key concept behind the distributional responses in both space and time can be found in the conservation of the thermal niche (or bioclimatic envelope, Ter Braak and Barendregt 1986).Whereby a shift in phenology, or geographical distribution, can be seen as a response that conserves the same thermal niche by changing timing in the year, or through a shift in latitude (Socolar et al. 2017, Beaugrand andKirby 2018).However, the distribution of abundance of a species (both in space and time) can also be affected by other physical and biological factors besides temperature, such as currents, food availability and predation (see for example Peterson et al. 2011, Fanjul et al. 2019, Uriarte et al. 2021).The thermal niche measured in situ, which is the focus of this study, can be defined as the "realized" thermal niche (Grüner et al. 2011).This differs from the more theoretical (at times laboratory derived) fundamental thermal niche, which is not the focus of our study.Following common convention (see Thackeray 2012), we use the term "phenology" to denote seasonal timing of metrics of copepod abundance, even though this is a balance between birth, death rates, and migration (at a given location), and thus not strictly comparable to classic phenological indices such as dates of first egg laying, bud burst, or flowering (Ji et al. 2010).
Phenological adjustment is one potential mechanism to conserve the same realized thermal niche in a changing environment.In fact, temperature can often explain up to $ 25-30% of the variability in timing of peak abundance from yearto-year in copepods (Mackas et al. 2012;Atkinson et al. 2015), with different intensities and directions of shifts observed in different species and at different latitudes (Edwards and Richardson 2004;Mackas et al. 2012, Usov et al. 2013;Reygondeau et al. 2015;Uriarte et al. 2021).One of the possible implications of shifts in the timing of seasonal copepod abundance is the trophic mismatch between the timing of consumer and resource abundance (Cushing, 1990).However, evidence for penalties incurred from a trophic mismatch remains scarce, and mainly apply to specialized taxa that show an abrupt increase or pulse of abundance (Thackeray, 2012;Atkinson et al. 2015, Samplonius et al. 2020).
Copepods make up $ 80% of the mesozooplankton abundance and have a crucial role in the marine ecosystems as the primary trophic link between unicellular primary producers and upper trophic levels, including commercially important fish species (Sundby 2000, Orlova et al. 2005).Copepods are also excellent candidates to explore responses to climate change since they are ectothermic and generally (in temperate and tropical regions) multivoltine, such that their abundance can rapidly respond to thermal changes.Multidecadal time series are essential to explore long-term patterns and year-to-year changes, while the high-frequency (weekly) sampling is valuable to detect shifts in timing.Several long time series (> 10 yr) are now available, mainly from the North Atlantic and European shelf seas (Mackas and Beaugrand 2010).Copepods are good exemplars of all three "universal" climate warming responses: shifts in range (poleward movement and expansion, Beaugrand et al. 2002;Chivers et al. 2017), in phenology (earlier peak of abundance, Mackas et al. 2012;Atkinson et al. 2015) and in body size (reduced adult body mass in warmer conditions, Rice et al. 2014;Horne et al. 2016;Corona et al. 2021).There is also increasing evidence for major and extensive declines in copepod abundance across the North-West European shelf (Boersma et al. 2015;Capuzzo et al. 2017;Schmidt et al. 2020).How these declines relate to shifts in phenology, range and size have not yet been examined.
In this study, we focused on the extent to which shifts in phenology help to conserve the thermal niche.We benefitted from three long time-series at stations that span a latitudinal gradient across the warming North-West European shelf, comparing indices of phenology across seven major copepod species that were sampled consistently within each site.The first question we address in this paper is do species conserve their thermal niche in response to climate warming, for example by adjusting their timing so that they increase at a more thermally suitable time of year?The second question is are phenological adjustments consistently greater for some species than others, or do the degree of adjustments differ across the range of a species?The third question is how do these timing shifts relate to the observed longterm decrease in copepod abundance?
Being zooplankton abundance affected by both food and mortality, timing metrics that are determined from population abundance are therefore different from commonly used phenology indices of many longer-lived species on land (for which such timing metrics are related to, for example, bud burst, first laying day, or day of flowering).For this reason, we also examined the egg production rate (EPR) of Calanus helgolandicus, at one of the sites to obtain a more mechanistic understanding of the processes that modulate the seasonal dynamics of copepods.EPR is not affected by mortality in the same way as abundance, which may therefore make it a better phenological metric and more comparable with indices used in terrestrial species (it is however unfortunately not largely available across multiple species and years).

Sample data
The sampling station chosen for this study were Helgoland Roads (HR; Helgoland, Germany), L4 (Plymouth, UK), and Stonehaven (SH; UK) (Fig. 1).These sampling stations were chosen for their similar high temporal sampling resolution (from weekly to three times per week), for their long time-series of data (from 20 up to 43 yr in duration), and because they span a latitudinal gradient, while still sharing many of the same copepod species.The seven copepod taxa chosen for this study were Acartia clausi, C. helgolandicus, Centropages typicus, Oithona spp., Paracalanus parvus, Pseudocalanus elongatus, and Temora longicornis.These are common and important species in this region, with contrasting traits such as size, feeding, and spawning mode, while being encountered across all the three of our selected stations (see Table 1 for summary of major traits).In the case of C. helgolandicus a small number of the congener Calanus finmarchicus may have been inadvertently included, as these occur rarely at the L4 station (a median composition of 4% of C. finmarchicus, according to Maud et al. 2015), and at HR and SH (unpublished data).Of course Oithona spp.may include three different congeneric species across these sites: Oithona similis, Oithona nana, and Oithona plumifera.However, O. plumifera is an Atlantic species and it is thus found very sporadically in the North Sea, whereas O. nana is the second most abundant cyclopoid in the North Sea, but it makes up only 7% of the Oithona spp.population at L4 (Fransz et al. 1991;Cornwell et al. 2018).Therefore, the great majority of Oithona spp.individuals in our study area can be attributed to O. similis.Some key descriptors of the locations, together with methods of sampling and data collected across the three stations are summarized in Table 2.
For the L4 site (http://www.westernchannelobservatory.org.uk), sampling methods are reported in detail by Atkinson et al. (2015), but in summary, sampling has occurred every week, subject to weather conditions, since 1988.For weekly zooplankton sampling, two vertical (50-0 m depth) replicate tows of a WP-2 net (56 cm ring diameter, 200 μm mesh) are made within the $ 54 m deep water column.Each sample is then fixed in 4% buffered formalin, from which sub-samples are obtained for counting and identification of taxa.We used the most recent dataset available (McEvoy et al. 2022).C. helgolandicus female egg production experiments have been conducted since October 1992, and are based on live material collected weekly with a 710 μm mesh in the surface layers at L4 and returned to the laboratory in a cool container.Further method details are provided in Atkinson et al. (2015) and in Maud et al. (2015).For consistency of measurement across the whole time series, we used surface temperature values measured with a mercury-in-glass thermometer in a stainless-steel bucket of freshly collected surface seawater.These values corresponded well to values measured in more recent years with an electronic probe for conductivity, temperature, and depth analysis (Atkinson et al. 2015).
At the Stonehaven (SH) station (https://data.marine.gov.scot), zooplankton has been sampled weekly since 1997 by vertical hauls from 45 m to the surface with a 200 μm mesh Bongo nets (40 cm ring diameter).Fixation, counting and identification of the samples were similar to those at L4.More details on the methodology used for data acquisition can be found in Bresnan et al. (2015) (zooplankton data source: Marine Scotland Science 2018; Scottish Coastal Observatory-SH station data.doi: 10.7489/610-1).Surface temperature at SH was measured using a digital reversing thermometer fitted to the Niskin sampling bottle.
At the Helgoland Roads (HR) station (https://deims.org/1e96ef9b-0915-4661-849f-b3a72f5aa9b1) zooplankton samples have been taken three times per week (weather permitting) by two oblique tows (net mesh size: 150 μm, diameter: 17 cm) since 1975 (mostly between 6 am and 10 am).Temperature was measured on a work-daily basis.More information on characteristics and methods of sampling can be found in Wiltshire et al. (2010).Unlike the other two stations, the abundance data of P. elongatus and P. parvus were not available, since these two species are summed up into a single taxon in HR (source of the dataset: Boersma et al. 2017).Further details of the copepod sampling, analysis, and trends are provided in Greve et al. (2004), andBoersma et al. (2015), including a phenological analysis in Wiltshire and Boersma (2016).

Data analysis
We used the abundance values (of all copepodite stages) from each time-series to obtain the phenological timing indices of each species at each station using the method of Greve et al. (2005).This method uses the Julian day at which a certain cumulative percentile of abundance (CPA) is reached within a year.We used the Julian days corresponding to the 25th and the 50th CPA days in each year as a timing index.These indices approximate respectively to the early and the middle of the copepod productive season and are often used in plankton phenology studies (Greve et al. 2005;Castellani et al. 2015).We also estimated the duration of the density peak (amplitude) as the number of days between the 25th and the 75th CPA days.When looking for timing shifts related to climatic warming, we linearly regressed each individual timing index and the amplitude value against spring (April, May, and June), summer (July, August, and September), and combined spring + summer surface temperature, in all cases measured as a mean value across the months.This approach was followed because spring temperature does not strongly correlate to summer temperature across years (R 2 = 0.09; p = 0.08, and n = 31).From these regressions we obtained slope values that indicate the number of days shifted per 1 C increase in temperature (positive slopes  2 for environmental and sampling summaries for the three stations. indicate a later timing with warming, negative slopes indicate an earlier timing with warming).We also linearly regressed each individual timing index against "year" in order to assess any long-term trends in copepod phenology.We calculated an index of the "realized thermal niche" center of each species in each year and at each station in order to detect the temperature around which the annual species-specific abundance is centered (see the "center of gravity" method: Colebrook and Robinson, 1965;Edwards and Richardson, 2004;Atkinson et al. 2022a,b), using the following formula: where R is the center of the realized thermal niche, expressed in Celsius degree ( C); N T is the abundance of each species at each sampling time point with sea surface temperature T.
We performed an ANOVA on measures of thermal niche index of all species, years, and stations (aov function in R, Chambers and Hastie 1992, data were normally distributed), followed by post-hoc pairwise comparisons (TukeyHSD function in R, Yandell, 1997) to tests the effects of the species and station variables.We divided each time series into two pairs of decades: the 10 warmest and 10 coldest years and the first 10 and last 10 yr, in order to compare copepod abundance distributions between distant chronological phases and between distinct thermal regimes (these were based on the average monthly temperature between June and September inclusive).We chose to use the maximum sample size possible (this was 20 since the shortest time series of any of the stations was 20 yr in length) for all stations given the high interannual variability of both temperature increase and copepod density change.These abundance distributions and their relative timing indices were plotted for each species at each station.Finally, in order to calculate the seasonal timing difference between first and final decades and between temperature regimes, we also performed an ANOVA test on the overall (all species pooled together) z-scores ([valuemean]/ standard deviation) of the timing indices of each year and stations.

C. helgolandicus egg production-case study
At L4, C. helgolandicus egg production rate (EPR) has been measured at weekly resolution, albeit with some data gaps, since 1992 (see Maud et al. 2015).Live net samples are transported to the laboratory inside a cool box and five replicates of five healthy adult females (25 in total, where available) are picked within 2-3 h of collection.These are incubated in a 500 μm mesh-bottom Plexigas chamber inside a 2 L plastic beaker filled with 1.5 L of 0.2 μm filtered seawater, at ambient L4 surface temperature and constant darkness for 24 h.Total numbers of eggs are counted and EPR calculated (eggs female À1 d À1 ).We calculated the timing indices of EPR in the same way we did for abundance Table 1.Summary of thermal preferences and traits for each species.Realized thermal niche (average across all years and stations) and phenology shift values are obtained from the present study (average across all years for L4 and HR, we excluded SH due to non-significant shifts and for sake of semplicity).Mean body mass and temperature-size response values are from Corona et al. (2021).Temperature (T) at max EPR values are from literature (Smith and Lane, 1985;Laabir et al. 1995;Nielsen and Sabatini, 1996;Halsband-Lenk et al. 2002;Castro-Longoria, 2003;Maps et al. 2005 (the Julian days of the 25th and the 50th percentile of EPR of that given year), and linearly regressed these against the phenological indices of C. helgolandicus.We also calculated the relative change over years in each month of C. helgolandicus EPR, along with total zooplanktonic predator density (calculated as the total ind m À3 sum of fish larvae,   Habour 2021 for source data used here).To examine seasonality among these variables, we log 10 transformed C. helgolandicus abundance, C. helgolandicus EPR, total zooplanktonic predator density, and total microplankton biomass.To further provide comparability across variables with different units we transformed them into z-scores.
Fig. 4. Linear regressions between copepod abundance and years for each species at each station (HR, L4, and SH).Points and lines are color indexed as indicated by the legend at the bottom where "rest of the year" refers to the months from October to April included.R 2 and slopes (β) are reported only for the significant regressions (p < 0.05).

Seasonality
The average seasonal copepod density at the station Helgoland Roads (HR) is characterized by one single main peak (between spring and summer), whereas at L4 and Stonehaven (SH), a bimodal density distribution is more common (Fig. 2). A. clausi, C. helgolandicus, C. typicus, and Oithona spp., at SH, present on average a more pronounced peak in summer than in spring when compared to the more southerly L4 station.P. elongatus abundance generally peaks in both spring and summer at the northerly station (SH), but only in spring in the south (L4).This summer preference at SH was also detected statistically by the post-hoc pair-wise test on species-specific timing indices between stations: five out of seven species had an abundance distribution significantly later in the year in SH compared to L4 (Supporting Information Fig. S2).The 50th cumulative percentile of abundance (CPA) day of all species (pooled together) across years differed significantly between each pair of stations, but the higher difference was detected between SH and L4 (difference in means: 27.38 Julian days; 95% CI: 14.86, 39.89, p < 0.01, Fig. 3a; Supporting Information Table S1 shows results of the same analysis but on the same time period for all stations: 1998 to 2018).Species-specific timing indices at SH significantly correlated with the ones in L4 (Fig. 3b), indicating a consistent difference in timing among species.This correlation was not present when comparing the species-specific timing indices between L4 or SH with the ones in HR.

Corona et al. Phenological adaptations in marine copepods
Abundance Six out of seven species showed a significant decrease in abundance over the years (Fig. 4; Supporting Information Fig. S1 also shows the long-term trend of abundance over years for each month).The only taxon that did not decrease at any station was Oithona spp.Looking in more detail, copepod abundances decreased at slightly different times throughout the year at different stations.As shown in Supporting Information Fig. S1, the decrease of abundance in HR occurred mainly in summer and autumn, whereas at Fig. 6.Change over years of different variables (indicated above each plot on the left-hand side) by each month at location L4.On the left-hand side, points indicate the slope value on the y-axis, calculated from a linear regression of z-scores against years, each point presents a vertical band that indicates the 95% CI (thus, bands that do not intersect the 0 [horizontal dashed lines] indicate a significant slope).Background colors refer to seasons (blue: winter, green: spring; yellow: summer; orange: autumn).Predators are calculated as the sum of fish larvae, chaetognaths, and cnidarians.Food concentration (mg C m À3 ) is calculated as the sum of diatoms, dinoflagellates, flagellates, ciliates, coccolithophorids, and Phaeocystis (the nanoflagellate functional group, which was not included, is also decreasing overall at L4 as described by Atkinson et al. 2022a,b).The plots on the right-hand side correspond to the variables on the left-hand side and show the seasonal distribution change over time in absolute terms across the first and last decade.

Corona et al.
Phenological adaptations in marine copepods L4 we found the strongest decreases often in the period between spring and summer, and between summer and autumn.At SH, significant decreases of abundance were observed in only a few species in a few periods (mainly P. elongatus during the cold part of the year and in September).

Timing of abundance with warming
At HR, three species (A.clausi, C. helgolandicus and Oithona.spp.) had at least one timing index (either the 25th or the 50th CPA days) negatively correlated with increasing seasonal temperature (T.longicornis showed a negative correlation, but this was slightly non-significant; p = 0.068).At L4, P. parvus showed a significantly delayed 25th CPA day with increasing surface temperature.At SH, no significant correlations were detected except for C. helgolandicus (for which the 50th CPA day got earlier with increasing spring surface temperature).The number of days between the 25th CPA and the 75th CPA days was never found to be significantly increasing for any species at any station, instead, it significantly decreased with warming in C. typicus and P. parvus at L4, in A. clausi at HR, and in T. longicornis at SH.There was no clear evidence that springabundant species tended to increase earlier in warmer years, or indeed that autumn species appeared later in warmer years (regressing timing shift in days C À1 ) vs. average day of appearance (average timing index day).C. helgolandicus showed a consistent and significant earlier occurrence of its seasonal timing of abundance against years at all stations.A. clausi, C. typicus and T. longicornis also showed timing advance but only at HR (Supporting Information Fig. S3).The amplitude of abundance (days between the 25th and the 75th CPA day) decreased over years in two species at HR (A. clausi and Oithona spp.) and one in SH (P.elongatus), and it increased in one species at L4 (C.helgolandicus).Overall, timing indices (25th and 50th CPA days) of all species at all stations changed more (became earlier) from the first to the last decade of the respective time series (F 1,377 = 23.35:p < 0.01, n = 378 for the 25th CPA day and F 1,349 = 24.59;p < 0.01; n = 378 for the 50th CPA day), than between the 10 coolest and the 10 warmest summers (25th CPA day: F 1,377 = 0.20; p = 0.65; n = 363 and 50th CPA: F 1,349 = 0.18; p = 0.66; n = 363; Fig. 5b-d).

Thermal niches
Realized yearly thermal niches (across all species together) differed significantly both among stations (F 2,454 = 67.40,p < 0.01) and among species (F 4,454 = 19.55,p < 0.01), but the difference among stations depended on the species (station Â species interaction term: F 8,454 10.62, p < 0.01).Not all seven species could be used, as HR lacks data on P. elongatus and P. parvus.The posthoc pairwise test detected a significant difference in the realized thermal niche means between L4 and SH for all species except Oithona spp., Supporting Information Fig. S2).Tukey post-hoc pairwise comparison showed that all realized thermal niches differed more between SH and L4 (difference in means = À1.9;95% CI: À2.77, À1.04; p < 0.01) than between L4 and HR (difference in means = À0.4;95% CI: À1.12, 0.29; p = 0.36) or SH and HR (difference in means = À2.3;95% CI: À3.13, À1.49; p < 0.01, Fig. 3c).Species-specific thermal niches in SH correlated significantly with the ones at L4 (Fig. 3d), indicating a consistent difference in thermal niches among species (absolute distribution of yearly realized thermal niches at each station for each species can be seen in Supporting Information Fig. S4).This correlation was not present when comparing the species-specific thermal niches between L4 or SH with the ones in HR.Over time, realized thermal niches did not change over time except in Oithona spp. in which it significantly increased (Supporting Information Fig. S5).

Timing of egg production in C. helgolandicus
To examine potential mechanisms behind phenology shift, we examined available time series data in more detail in C. helgolandicus at L4. Figure 6 shows long-term seasonal changes for all the four summary variables relating to its population dynamics: C. helgolandicus abundance, egg production rate (EPR), biomass of potential food, and abundance of its potential predators.C. helgolandicus abundance decreased over the 30 yr in May, June, and July: exactly a month delayed to the long-term decrease of C. helgolandicus EPR, which occurred in April, May, and June.Food (as total microplankton [mg C m À3 ]) decreased significantly in July and slightly significantly in April, whereas predators abundance (ind m À3 ) increased significantly between winter and spring and in late autumn.

Thermal niche conservation
One important, but often untested assumption of species distribution models that utilize the thermal niche to predict shifts in space and time, is that the relationship of a species to its environment (and particularly to temperature) is broadly fixed (Cheung et al. 2008;Brun et al. 2017).We tested this assumption using the species realized-thermal-distributions at the three different stations and found that this assumption was not upheld among stations.For this reason, it would have not been optimal to use the lab-measured thermal niches as a reference (however their values are reported in Table 1).Comparison between L4 and Stonehaven (SH) shows, generally for all species, a more summer-centered abundance distribution in SH, whereas in L4, warmer temperature allows them to take advantage of spring too.This difference can be seen both in the seasonal abundance distributions and the seasonal relative percentage abundances and in their average timing indices (the 25th and the 50th cumulative percentiles of abundance (CPA) days), which are later in the year in SH than in L4 (Fig. 3 and Supporting Information Fig. S2).From this, we may assume that if the timing of seasonal abundance can change from one latitude to another, as a local plastic response, it can also potentially change over time within the same latitude as sea temperature rises, although this depends on how close the species live to their thermal tolerance limits (Pinsky et al. 2019).Furthermore, variables other than temperature (such as food and salinity, for instance) may explain some of the interlatitudinal variance in seasonal abundance distribution (Uriarte et al. 2021).However, this seasonal timing adjustment between different latitudes is not enough to maintain the same realized thermal niche, because SH copepods still experience colder temperature than at L4 (Fig. 3).One reason for this is that the maximum surface temperature in SH is not as high as that at L4, and the time gap to take advantage of warmer conditions is shorter at higher latitudes.The distribution of the thermal niches within the thermal limits of each station (Supporting Information Fig. S4) clearly shows how copepods at L4 manage to exploit a much larger thermal range (they can thrive at both above and below the local average surface temperature), whereas in SH they can only maximize their abundance in the warmer part of the local thermal range.This narrower thermal range in SH may explain (besides other possible biotic and abiotic factors) why the mean yearly abundance was lower than in L4: all species (except for A. clausi) reached significantly higher numbers at L4 than SH (Supporting Information Fig. S4).Interestingly, the way copepods change their seasonal timing between L4 and SH remains consistent interspecifically.In other words: "cold species" and "warm species" at L4 remain relatively the same "cold species" and "warm species" at SH, and the same applies to "early species" vs. "late species" too (Fig. 3).This pattern could mean that, despite having a narrower season length at SH, of the copepod species examined, they cannot greatly overlap their thermal niches with each other, following the principle of competitive exclusion (Gause 1934).Helgoland Roads (HR) has higher and lower temperatures than at both L4 and SH, so that copepods tend to have a larger thermal range at their disposal.This could explain why copepods at HR have thermal niches intermediate between those at L4 and SH.Additionally, copepods at HR present a different distribution of seasonal abundances (uni-modal) than at L4 and SH (bi-modal), which could obscure the comparisons (interspecific differences in thermal niches were inconsistent between HR and the other two stations).
Overall, all species at all stations have conserved the same realized thermal niche over time within each station, with the only exception being Oithona spp., whose thermal niche has significantly increased over time in HR (Supporting Information Fig. S5).However, it has to be acknowledged that possible phenological and thermal niche adjustments of Oithona species may have been obscured at the three sites: first because identification was to genus level, rather than species level; and second, given the relatively small size of this taxon, juvenile abundance, and timing may have been underestimated at L4 and SH, which was sampled with a coarser mesh size net (200 μm) than at HR (150 μm).An increase in abundance in summer (relative to other periods) indicates an increase in the temperature of the realized thermal niche of that year, and yet, species which instead decreased in abundance in summer did not show an opposite trend (i.e., a thermal niche, which gets colder over time due to fewer individuals in summer, instead, it remained the same throughout the whole time series).One explanation for this is that while the individuals during the warm season have decreased, their overall realized thermal niche has not changed because temperatures have increased over time.This could mean that the decrease of abundance in summer has been a response to the warming, which has led to these species maintaining the same thermal niche (measured as the yearly center value of the overall temperatures experienced by all individuals).From this simple concept come two important considerations: first, the seasonal timing changes at HR may simply be a consequence of the decrease of abundance in late summer/early autumn.This could also explain why we observed fewer or no significant shifts in timing index at the other two stations where the decrease of abundance occurred differently (in two distinct parts of the year instead of one; Fig. 5).The second consideration is that this thermal niche conservation was associated with possible costs, because although similar temperatures are still experienced by copepods, their abundance has significantly decreased over the time series.

Varying influence of temperature on phenology across the range of a species
Our second hypothesis, based around the predictions of Beaugrand and Kirby (2018), is that a species should have a different response in their seasonal timing to temperature at different latitudes or in different parts of their thermal range.Specifically, the "width" of the abundance peak of a species should broaden with warming (with a resulting increase in the annual abundance) at the cold part of its thermal range, and get narrower with warming at the warm part (with a resulting decrease of annual abundance).By contrast, the same species at the center of its thermal range is predicted to reduce in abundance during summer under warming, with a tendency to increase its abundance earlier in the year (see Fig. 4 in Beaugrand and Kirby 2018).The species examined here did not show longer duration of high abundance with warming at any station (there were actually cases of decrease), and did not show an increase of annual abundance over time (rather they mostly show a decrease instead).This could indicate that the stations considered do not represent the northern range (the minimum temperature limit) of the species, not even the coldest station (SH).In fact, all the species in our study can be found at higher latitudes than SH (Halvorsen and Tande 1999;Bucklin et al. 2000;Potters et al. 2004;Beaugrand et al. 2007;Evjemo et al. 2008).This may indicate that our stations are somewhat within the central part of the thermal niche of our species, where the calculation of phenological index should reveal an earlier phenology, and the annual abundance should not change substantially, according to the Beaugrand and Kirby model.In our study, instead, we found significant changes in timing indices with temperature almost exclusively in HR (Supporting Information Fig. S3), and we found longterm decrease of abundance at all stations.
It is also difficult to determine whether the absence of multiple and strong significant timing shifts in species at SH is due to limited number of years of sampling, or due to intrinsic environmental characteristics of the region.Uriarte et al. (2021), for instance, suggests that phenology at SH is driven more by chlorophyll than by temperature.Stronger timing changes in HR could also be due to the unimodal seasonal distribution of abundance there compared with the bimodal distributions in L4 and SH.In fact, having a single and wide peak of copepod density may makes abundance a better proxy for phenology, whereas having a decline in density between the two peaks could imply a greater seasonal mortality, which can potentially obscure the real phenological changes (and the thermal niche) of copepods.
Overall, the timing indices across all species and stations changed more (species abundance was distributed earlier in the year) from the first to the last decade of the respective time series, than between the top 10 coolest to the top 10 warmest summers.Therefore, it is likely that variables other than temperature may be regulating copepod seasonal abundance distribution (as found in Fanjul et al. 2018), or that the effect of temperature increase is indirect and can be detected only on a very long-term scale.Surface temperature has a very high interannual variability, as illustrated by regressions of annual "surface temperature vs. years" regressions in all the stations, which have shallow slopes and high scatter in HR (β = 0.03 C yr À1 , R 2 = 0.28) and L4 (b = 0.02 C yr À1 , R 2 = 0.21), albeit significant.Moreover, at SH the surface temperature shows no obvious temporal warming in the 1988-2019 time-series (β = 0.004 C yr À1 , R 2 = 0.01, p = 0.07), despite warming occurring globally (IPCC, 2021).Another confounding factor is the constant decrease of copepod abundance over years, which does not occur homogeneously throughout the year, but more in some months than others.It should also be acknowledged that fixed sampling stations are also vulnerable to great variability and error in the data due to the advection of "foreign" individuals.Therefore, studies that rely on this type of data source must assume that each station is representative of the mean population dynamics over a broad regional area.

Drivers of population dynamics and their influence on phenology
Abundance of copepods across the North-West European shelf and North-East Atlantic have declined greatly over the last 60 yr (Boersma et al. 2015;Capuzzo et al. 2017, Schmidt et al. 2020).Are these changes a cause or a consequence of the lack of clear phenological adjustment with temperature that we see?The abundance distributions of each species at L4 have changed greatly: from a uni-modal distribution to a bi-modal one (Fig. 5), as a result of abundance decrease between spring and autumn (Supporting Information Fig. S1).This change can be seen both between thermal regimes and between the first the final decades, but not in the much shorter and more recent time series at SH, which has a bi-modal distribution throughout its whole time series.At HR however, copepods decreased in abundance mainly in summer, but maintained a uni-modal distribution for each species both over time and over different thermal regimes.Having two peaks of abundance during the year can significantly increase the complexity and difficulty in analyzing seasonal timing, especially when the number of peaks also changes over time, as seen in L4.Therefore, this probably could also explain why we see stronger seasonal timing changes in HR, where the peak was uni-modal for the whole time-series, making abundance distribution a more reliable proxy for phenology.
Given the intrinsic complexity of abundance alone as a metric of phenology, we tested the relationship of timing of abundance with the timing of a life cycle event that is less influenced by mortality: egg production rate (EPR).For C. helgolandicus, the seasonal timing indices of abundance did not correlate with those of EPR, as also found by Maud et al. (2015) and Cornwell et al. (2018).This lack of correlation may indicate that the drivers of EPR and abundance are different, probably because of the strong effect of mortality on abundance (and not on per capita EPR).However, Fig. 6 shows that the C. helgolandicus abundance decrease over time occurs in a period of the year, which is a month delayed from the decrease of the EPR.This delay (of $ 1 month) roughly coincides with the average development time of this copepod species, indicating a possible (although speculative) relationship between the decline of abundance and egg production on a long-term scale.There does not seem to be a relationship between seasonal increase in EPR and increase in adult abundance, however, as found by Maud et al. (2015) and Cornwell et al. (2018) for C. helgolandicus, a surplus of eggs is not guaranteed to provoke an increase of adult abundance, due to mortality.Conversely, a deficit of eggs is an obvious prerogative of adult deficit, therefore, what is causing this seasonal decrease in EPR? Again, we speculate that this may be due to a decrease in the environment quality for C. helgolandicus in almost all months of the year, with particularly abundant predators and scarce food in the warmer summer months.This is manifested in reduced spring/summer EPRs and lower abundance, but the shift to later egg production may have sustained the population density in autumn and winter.Similarly, Schmidt et al. (2020) showed from the Continuous Plankton Recorder survey dataset, that microplankton significantly decreased over a period of 60 yr during the months from May to September.However, at L4 we noticed that the decrease of EPR over the study period started in April, which coincides with a decrease of microplankton and a significant increase of predators (Fig. 6).Indeed, the month of April appears to be a crucial time at L4, as the beginning of the temperature-dependent phase, when copepods tend to show strong negative correlations between adult body size and temperature, most likely because of high ratios of predators to available food, and potentially high and varying degrees of competition for food (Corona et al. 2021).Overall, the compound effects of warmer summer temperatures, and changes in food and predators, are likely to have negative effects on copepods across the North-West European shelf, but their relative roles will likely vary regionally.

Conclusion
Both temperature-related phenological shifts and to some extent a degree of thermal niche conservation were evident within our study sites.However, the effects were far from clear, and obscured by major within season changes in abundance, these changes are likely related to station-specific balances of food and predation controls.These non-temperatures related drivers of population dynamics severely challenge the ability to understand and predict shifts in plankton phenology in relation to temperature.Moreover, the earlier seasonal timing was demonstrated more clearly across years than with temperature.Therefore, it is difficult to determine if the actual effect of global warming on plankton phenology is obscured (and possibly underestimated), or simply not directly correlated with seasonal distribution of species abundance.However, time series are lengthening worldwide and are becoming ever-better networked across the North-West European shelf.Linking these into coherent analyses will unravel the drivers of the profound, climate-related changes that we are seeing in the pelagic food web.

Fig. 1 .
Fig. 1.Location of the three stations (L4, HR, and SH) in the North Sea and English Channel.See Table2for environmental and sampling summaries for the three stations.
cnidarians, and chaetognaths) and total microplankton biomass (calculated as the total mgC m À3 sum of diatoms, dinoflagellates, flagellates, ciliates, coccolithophorids, and Phaeocystis).These microplankton biomass values are based on Lugol's preserved water samples from 10 m depth, with inverted microscopy used to determine taxonomic groups and compute carbon from biovolume estimates specific to the station(Widdicombe et al. 2010; seeWiddicombe and

Fig. 3 .
Fig. 3. (a) Phenological indices values (50th cumulative percentile of abundance (CPA) day) of all years and species pooled (except P. elongatus and P. parvus, since these two species are not available for HR station) for each station (HR, L4, and SH), in forms of boxplots (horizontal thick lines indicate the median value, the lower and upper end of the box correspond to the 25th and 75th percentiles, whiskers indicate approximately the 95% confidence interval of the median).The p values are reported for the pairwise significant differences between stations according to the ANOVA test.(b) Reduced major axis regressions between the L4 and SH stations of the specific mean phenological indices: the 50th CPA day.Symbols indicate different species (see legend), horizontal error bars indicate the standard deviation of the specific mean values, red line is the fitted line of regression (95% CI is indicated by the red dashed lines), correlation coefficient (R) and p values are reported within both graphs.The black dashed line is a 1 : 1 line added for better interpretation.(c, d) are respectively analogous to (a, b), but they report realized thermal niche values ( C).

Fig. 5 .
Fig. 5. (a,b) Box-plots represent the distribution of the 50th cumulative percentile fo abundance (CPA) day timing index (z-scores scaled by species and stations), showing the direction between different chronological phases (top) and different thermal regimes (bottom): negative values on the y axis mean earlier phenology, positive values mean a later phenology (only the difference between chronological phases is statistically significant for both the 25th and the 50th CPA day, unlike the difference between thermal regimes in both timing indices).(b-d) Seasonal abundance of each species denoted with curved lines which are smoothed averages from the "loess" function with a span level of 0.75, for each station (HR, L4, and SH).Vertical lines denote 50th CPA day timing index.Color indexed denotes either the chronological phase (top graph) or the thermal regime (bottom graph).Gray vertical dashed reference lines indicate winter, spring, summer, and autumn.

Table 2 .
Environmental characteristics of each station and corresponding sampling methods used.
Fig.2.Average yearly copepod density distributions are indicated by black lines (smoothed average) for each species at the three stations SH, L4, and HR.Density values are reported on the y axis in a log scale.The background colors of each vertical band refer to the four seasons (blue: winter, green: spring; yellow: summer; orange: autumn).White boxplots indicate median and distribution of monthly density values (white boxes indicate the 50% quartile and whiskers indicate the 95%).The numbers reported in the boxes in each season above the curve indicate the seasonal percentage of mean density relative to the whole year mean (the color of the boxes correlate to the value of these percentages, as indicated by the color scale at the bottom), thus, the values in the boxes show how successful (in terms of density) copepods are in a given season compared to other seasons.Thermal niche conservation is demonstrated in the always higher summer relative abundance in SH than in L4, and the opposite relationship for the spring relative abundance.