Transient demographic dynamics of recovering fish populations shaped by past climate variability, harvest, and management

Large‐scale commercial harvesting and climate‐induced fluctuations in ocean properties shape the dynamics of marine populations as interdependent drivers at varied timescales. Persistent selective removals of larger, older members of a population can distort its demographic structure, eroding resilience to fluctuations in habitat conditions and thus amplifying volatility in transient dynamics. Many historically depleted marine fish stocks have begun showing signs of recovery in recent decades following the implementation of stricter management measures. But these interventions coincide with accelerated changes in the oceans triggered by increasingly warmer, more variable climates. Applying multilevel models to annual estimates of demographic metrics of 38 stocks comprising 11 species across seven northeast Atlantic ecoregions, this study explores how time‐varying local and regional climates contributed to the transient dynamics of recovering populations exposed to variable fishing pressures moderated by management actions. Analyses reveal that progressive reductions in fishing pressure and shifting climate conditions discontinuously shaped rebuilding patterns of the stocks through restorations of maternal demographic structure (reversing age truncation) and reproductive capacity. As the survival rate and demographic structure of reproductive fish improved, transient growth became less sensitive to variability in recruitment and juvenile survival and more to that in adult survival. As the biomass of reproductive fish rose, recruitment success also became increasingly regulated by density‐dependent processes involving higher numbers of older fish. When reductions in fishing pressure were insufficient or delayed, however, stocks became further depleted, with more eroded demographic structures. Although warmer local climates in spawning seasons promoted recruitment success in some ecoregions, changing climates in recent decades began adversely affecting reproductive performances overall, amplifying sensitivities to recruitment variability. These shared patterns underscore the value of demographic transients in developing robust strategies for managing marine resources. Such strategies could form the foundation for effective applications of adaptive measures resilient to future environmental change.


| INTRODUC TI ON
Excessive human exploitation that persisted through the 20th century depleted marine living resources worldwide, threatening biodiversity, food and nutrition security, and other ecosystem services that the oceans provide (Crowder et al., 2008;Pauly et al., 2005;Worm et al., 2009).Long-lived, late-maturing vertebrate predators in particular were disproportionately depleted by industrial fishing as target or incidental catch (Christensen et al., 2003;Juan-Jordá et al., 2022).
As part of the efforts to halt and reverse the losses, national government and intergovernmental organizations began applying the precautionary principle to the management of marine fisheries (Hilborn et al., 2001;Punt, 2006), following international agreements such as the Code of Conduct for Responsible Fisheries adopted by the Food and Agriculture Organization of the United Nations (FAO, 1995;Melnychuk et al., 2021).Since the implementation of stricter management measures to regulate fishing operations in the early 1990s, an increasing number of historically overexploited fish stocks began showing signs of recovery (Hilborn et al., 2020;Juan-Jordá et al., 2022).
These management interventions, however, coincide with accelerated changes in the oceans (Barange et al., 2018).Global climate is becoming increasingly warmer and more variable, modifying regional and local habitat conditions that support the biomass production of fish stocks and possibly shifting ecologically sustainable levels of harvesting (Chassot et al., 2010;Free et al., 2019).
The dynamics of marine resource populations are shaped by drivers operating interactively at varied timescales, including low-frequency variability in ocean properties forced by climate and high-frequency perturbations like fishing (Goto, Filin, et al., 2022;Rouyer et al., 2012).
When exposed to such persistent, nonstationary perturbations, nearterm (transient) dynamics of an exploited population, which depend on both variable vital rates (survival, growth, and fecundity) and demographic (age and stage) structures, may deviate from asymptotic (equilibrium) dynamics (Hastings, 2004;Koons et al., 2016).Past research instructs us that large-scale shifts in ecosystem productivity under alternative climate regimes can propagate through the system to influence fluctuations in fisheries catches (Chassot et al., 2010;Cheung et al., 2013).Climate effects on marine species are often integrated over a series of oceanic, ecological, and biological processes over space and through time (Ottersen et al., 2001;Stenseth et al., 2003).
Fluctuations in population structure and size can, for example, emerge with some time lag from age (or stage)-dependent demographic processes like birth and death (Jenouvrier et al., 2022) modulated by seasonal and interannual variations in weather and habitat conditions (Hutchings & Reynolds, 2004;Lindegren et al., 2013).These transient fluctuations in the population could sustain over generations through delayed feedbacks (Hastings, 2004).
In contrast, fishing operations and management actions that regulate the operations can directly modify short-term trends and variability in population structure and demographic rates, thereby modulating the transient dynamics of exploited species (Gamelon et al., 2019).
Large-scale commercial fishing often selectively removes larger, older members of the population, truncating age (and size) structure, modifying life-history traits, reducing the capacity to absorb the effects of perturbations, and thus amplifying year-to-year variation in population size (Hsieh et al., 2006;Rouyer et al., 2012).Fishing-induced age truncations (juvenescence) are widely reported for harvested marine populations (Hsieh et al., 2010;Rouyer et al., 2012), many of which suffer from reduced reproductive capacity (Hixon et al., 2014).
Climate-induced changes in the variability of ocean habitat conditions thus may indirectly promote or impede the recovery of depleted populations with weakened demographic resilience to perturbations (Harris et al., 2018;Lindegren et al., 2013).
In commercial fisheries, a harvested population is commonly managed through the regulation of fishing pressure based on its population status (Beddington et al., 2007).Because of high socioeconomic values, the exploitation of marine resource populations is regularly monitored through scientific surveys and assessed based on population size estimates (Hilborn & Walters, 1992).These population assessments form the scientific basis for the implementation of management measures like catch limits.When population size declines below a predefined threshold (a reference point set under the precautionary principle, ICES, 2021;Goto, Devine, et al., 2022), for example, fishing pressure may be adjusted to mitigate the risk of overharvesting.This adaptive management cycle (human-natural system feedback through monitoring, assessment, and management action) thus can be considered a natural experiment (Hilborn & Walters, 1992;Jensen et al., 2012), providing opportunities to gain insights into the transient dynamics of managed populations in variable environments.Understanding what drives these transients would help develop robust strategies to sustainably manage human exploitation of marine resource populations as we face increasingly uncertain conditions in the oceans.
This article reports a study that empirically explores how past fishing, management action, and climate variability influenced the transient dynamics of recovering fish populations in the northeast Atlantic Ocean (FAO Major Fishing Area 27, Figure 1a).These populations experienced wide ranges of fishing pressure, climate, and management action in recent decades (Rouyer et al., 2011;Zimmermann & Werner, 2019).Annual fishing pressures applied to many northeast Atlantic stocks began to be eased in the late 1990s to rebuild biomass that supports sustainable harvesting (Lassen et al., 2014;Sparholt et al., 2007).Following the reductions in fishing pressure, these intensively managed stocks showed considerable changes in population size, demographic structure, and productivity (Britten et al., 2016;Hidalgo et al., 2012;Rouyer et al., 2012).During this period, the regional climate also became increasingly warmer and more variable as greenhouse gas emissions rose, which was reflected especially in upper ocean heat content and transport (Delworth et al., 2016;Zanna et al., 2019).Warming patterns, however, highly varied among geographic subregions (ecoregions) owing to local airsea feedbacks (Delworth et al., 2016;Roemmich et al., 2015), triggering a series of oceanic and ecological responses at variable rates (Beaugrand, 2009;Delworth et al., 2016).To glean insights into the recovery patterns of historically depleted fish stocks, this study specifically asks how local and regional climate variability contributes to variations in the reproductive performance and transient growth of marine species exposed to variable fishing pressures moderated by management actions.

| MATERIAL S AND ME THODS
To investigate the transient dynamics of northeast Atlantic fish stocks at local and regional scales, a series of multilevel statistical models were applied to the time series of population and demographic metrics along with intrinsic and extrinsic covariates.

| Data
Time series  of annual estimates of age-specific information (mean body mass, abundance, maturity rate, and natural and fishing mortality rates) and management reference points (which define management targets, as described below in Management action for fisheries) of marine fishes in the northeast Atlantic Ocean were extracted from the 2017 stock assessment reports by the International Council of Exploration of the Sea (ICES; Table 1).
Climate indices and sea surface temperature (SST) anomalies were used as proxies that represent variations in regional (North Atlantic) and local (ecoregion) climate conditions (during 1946-2016), which are filtered through oceanic processes fish populations directly and indirectly experience over interannual and decadal timescales (Ottersen & Stenseth, 2001;Stenseth et al., 2003).SST was chosen as a proxy for local habitat conditions because its reconstructed data can provide adequate spatial and temporal coverage for this study.Furthermore, the climate data were used as covariates in multilevel models with early life-history processes as response variables (as described below in Covariates of recruitment success and Covariates of elasticity to recruitment success).Climate effects on these processes can, for example, be mediated by lower trophic level (like zooplankton) production responding to regional climate-driven SST anomalies (Beaugrand & Kirby, 2018) and other ocean properties (like heat transport), which may in turn manifest as delayed effects in young fish survival (Hjermann et al., 2007).
Monthly time series of three regional climate indices were downloaded and averaged over winter months (when the atmosphere is most active in the region, Stenseth et al., 2003) for analysis (Figure 1b https://psl.noaa.gov]and averaged over the sources), and upper (0-700 m) ocean heat content (OHC) anomaly in the North Atlantic (NOAA National Centers for Environmental Information; https:// www.ncei.noaa.gov).The NAO index reflects the dominant mode of atmospheric variability between the Azores and Iceland, linked with a range of large-scale oceanic features (such as surface temperature and heat content) in the North Atlantic that regulate biological and ecosystem processes (Stenseth et al., 2003).The AMO index, computed by averaging SSTs over the North Atlantic (25°-60° N at 7°-70° W) using the method by van Oldenborgh et al. (2009), also reflects large-scale climate variability that may influence the movement and migration of marine species (Faillettaz et al., 2019).OHC measures the amount of excess heat (from increasing greenhouse gases) trapped in the ocean, integrating basin-wide heat storage and transport (Zanna et al., 2019).Because of high signal-to-noise ratios, OHC reflects more consistent decadal trends (less influenced by natural variability) than other climate indices (Roemmich et al., 2015).
Reconstructed gridded monthly SST data were downloaded from and averaged over three sources to account for systematic uncertainties in each source: HadISST1, COBE sea surface temperature (SST; NOAA physical sciences laboratory), and ERSSTv5.I first averaged the SST data over ecoregion and for each ecoregion over season; winter (December in the previous year to February, wSST anom ), spring (March-May, spSST anom ), summer (June-August, smSST anom ), and autumn (September-November, aSST anom ).I then computed SST anomalies by subtracting the mean across years  for each season and ecoregion (Figure 1c).

| Transient population growth and elasticity
I used transient population growth as an integrated measure that captures time-varying demographic rates and age structures responding to non-stationary extrinsic factors (fishing pressure, management action, and climate) and estimated the growth rates for each stock using an age-structured matrix projection model (Caswell, 2001).Leslie (state transition) matrices (A) were first constructed with annual estimates of age-specific survival and reproductive rates.Survival rates (S a,y ) were estimated as S a,y = exp(−Z a,y ), where Z a,y = F a,y + M a,y , and F a,y and M a,y are the fishing and natural TA B L E 1 Stock-specific information and management measures used in the assessments of 38 fish stocks in the northeast Atlantic Ocean included in the study.reference points set for each stock as target fishing mortality rates that produce maximum sustainable yields and biomasses that satisfy the precautionary approach; and status indicates recovery status assigned based on spawner biomass relative to B pa for analysis (R: recovering and NR: non-recovering) in this study.

Ecoregion
a F msy was not defined for Barents Sea Saithe in 2017 and an alternative reference point applied (F mp ) were used for analysis in the study. b The reference points for North Sea turbot were not defined until 2018 and the stock was excluded from multilevel modeling analysis.where n y and A y are a vector of initial age-specific numbers and an a × a state transition matrix in year y, and T gen,y is the generation time (years) in year y.To account for variations in life-history strategy among stocks, the transient growth rates of each stock were estimated over its generation time, which is defined as the mean age of parents of a cohort (Caswell, 2001), The transient growth rates here were thus defined in timescales of 2.6-6.7 years (Figure S1), relevant for fisheries management.
To evaluate the contributions of demographic variations to transient growth I did transient life table response experiments (LTRE) by numerically estimating transient elasticities (e y ) to survival and fecundity rates through matrix projections (Koons et al., 2016); demographic parameters (fecundity or survival) of all ages were first perturbed by a small amount (1%), and stock numbers were projected to re-compute λ trans,y (Durant et al., 2013).Elasticity (relative sensitivity), defined as the proportional change in λ trans relative to a proportional change in a demographic rate (Caswell, 2001), was computed as e y = Δlog trans,y

Δlog y
, where α y is the survival or fecundity rate in year y that is perturbed (Coulson et al., 2004).Age-specific elasticities were then summed for analysis; changes in elasticities to juvenile and adult survival rates were evaluated separately by summing age-specific elasticities weighted by maturity rates.To quantify how the transient dynamics deviate from expected long-term equilibriums (asymptotic dynamics assuming stable environments and age structures) across years I also computed year-specific asymptotic (geometric) growth rates and stable age distributions from the logarithms of the dominant (largest) eigenvalues and the right eigenvectors (respectively) of the annual state transition matrices (Caswell, 2001).
To evaluate how variations in demographic structure had contributed to the reproductive performance and recovery patterns of fish where p a,y is the proportion of a-year-old spawner number in a given age group in year y (Marteinsdottir & Thorarinsson, 1998).

| Management action for fisheries
Management action (here adjustment of fishing pressure) is based on perceived current stock status and harvest control rule set by the ICES (which is designed to achieve sustainable harvesting without compromising the reproductive capacity of the stock).The current status is assessed using scientific monitoring and catch data before the provision of annual catch advice (total allowable catch).Under the control rule, when the estimated SB remains above a fixed threshold (B pa ) at the start of the year following the assessment, the catch limit is computed based on target exploitation rate (F msy ).
These two parameters (B pa and F msy ) are designed to prevent overharvesting by accounting for uncertainty in population and harvest dynamics (ICES, 2019a).When the SB falls below B pa , exploitation rate is adjusted to F msy scaled to the proportion of SB relative to B pa , thereby allowing the population to rebuild.For analysis, fishing pressure and SB relative to stock-specific reference points (F /F msy and SB/B pa ) were used as proxies to evaluate the influence of the management action on the demographic metrics of fish stocks.
Fishing pressure is defined here as the instantaneous rate of fishing mortality (F), which is estimated by averaging age-specific fishing mortality rates over dominant ages in fishery catch (Table 1).

| Data analysis
To explore how management-guided changes in fishing pressure and climate change influenced the vital rates, transient dynamics, and recovery patterns of the northeast Atlantic fish stocks, I applied multilevel (hierarchical) modeling to derived demographic metrics to account for nonindependence in data structure (Gelman & Hill, 2006) with stock, ecoregion, and year as grouping factors (random intercepts).Data exploration including identification of collinearity (through visual inspection correlation plots and computation of variance inflation factors) (1) | 6025

GOTO
and outliers (>three standard deviations) was done before analysis (Zuur et al., 2009).The data exploration revealed two groups of stocks with diverging (increasing and declining) patterns in stock status based on relative spawner biomass (SB/B pa ).To explore correlations underlying the differences I applied generalized additive mixed effects models (GAMMs, Wood, 2017) to these groups separately.In the following analyses the stocks were grouped based on stock status in the last 15 years of the time series (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016); if spawner biomass is larger than 80% of a biological threshold (B pa ) (Costello et al., 2016) for at least 50% of the period (8 years) a stock is classified as "recovering," otherwise "non-recovering" (sensitivity analyses showed that the diverging patterns in SB/B pa are robust to variations in the criteria used, Figures S2 and S3).This classification resulted in 27 recovering and 10 non-recovering stocks (Table 1).for year; stock i , ecoregion j , and year jk are normally distributed random intercepts with mean zero and standard deviation σ for stock, ecoregion, and year nested within ecoregion.The structure of random effects was selected based on likelihood ratio tests applied to 10 structures fitted with a restricted maximum likelihood method (REML, Table S1).
Abundance and demographic parameters in stock assessments are estimated with several known sources of uncertainty (including observation, process, and estimation errors, ICES, 2019b), which may bias post hoc analyses when used as data (Brooks & Deroba, 2015).2a), the following model was fitted to the data:  where ф is a parameter (>0), and β 0 , f 1 -f 5 , stock , ecoregion , and year , are the same as above (here tensor product construction for time-varying covariate smoothers was modeled using the te function in the R package mgcv).
All numeric covariates were standardized by subtracting the mean and dividing by the standard deviation (z-scores).In model selection intrinsic covariates were first added, followed by indirect fishing, and then climate.With the random effect structure selected (as described above, Table S1) the models were tested sequentially for covariates by fitting with a maximum likelihood method, followed by parameter (re) estimation with REML (Zuur et al., 2009) for the model with the lowest Akaike information criterion (AIC) score and highest model weight (w).
AIC is computed as AIC = −2ln(L) + 2k, where L is the likelihood of the model and k is the number of parameters, and the model weight of model i (w i ) in a set of J models is computed as As diagnostics, model residuals were checked for homogeneity and regressed against covariates that were both included and not included in the selected models (Zuur et al., 2009).All analyses were done in R (ver- Because fish populations in different ecoregions experience a wide range of climate conditions and often have adapted to local environments, fish-covariate patterns may deviate from shared patterns captured by the regional models.To test this hypothesis, ecoregion-scale analyses were also done using the same model structures for RS total and e rec and the model selection steps as above except for random effects.Because random effect estimates become less stable when the number of groups is fewer than eight (Bolker, 2015), the random effect of the ecoregion-scale models consisted of year only, with stock being a fixed effect, except for the Celtic Seas; the random effects of the Celtic Seas model consisted of year and stock.

| RE SULTS
Analyses were designed to primarily address (a) whether long-term changes in spawners (biomass and age structure) shared by the Although climate effects varied among ecoregions, most showed shifting responses in recruitment (from positive to negative) and in elasticity to recruitment (from negative to positive) to increasingly warmer than average local and regional climates.
In contrast, relative fishing pressures on the non-recovering stocks kept rinsing until the late 1990s (to F/F msy = ~3.7 on average) before declining (~3.0%per year, ΔAIC = 63.3), but their relative spawner biomass continued to decline (ΔAIC = 92.7, Figure 2b).Although the mean spawner age of the non-recovering stocks began rising with declining relative fishing pressures, no change in age diversity was detected (Fig-

ure 2b
). Direct effects of selective fishing consistently emerged as negative correlations between relative fishing pressure and spawner biomass and age structure in both stock groups (ΔAIC = 15.1-165.0, Figure 2c).
Contrasting temporal patterns between relative fishing pressure and spawner features were more pronounced in some ecoregions as captured in year-specific (nested within ecoregion) random effect (variance) estimates, especially the Norwegian-Barents Seas, which showed abrupt shifts in the early 1990s (Figure S5).
The time series of transient population growth rates reflected the contrasting patterns between the recovering and the non-recovering stocks, and variations among ecoregions (Figure 3; Figure S6).Sensitivity analyses showed that multilevel model parameter estimates in the trend analyses for both stock groups were robust to uncertainties in select abundance and demographic rates used as input (simulated parameter values < two standard deviations of the estimates, Fig-

ure S7
).The regional average growth rates of the recovering stocks linearly increased over the study period (ΔAIC = 15.0, Figure 3).Although the relative spawner biomass of the non-recovering stocks continued to decline, their average growth rates began rising in the late 1990s (ΔAIC = 2.8, Figure 3) when the relative fishing pressures began declining (Figure 2b).These transient dynamics differed from asymptotic dynamics; the transient growth rates were lower than asymptotic growth rates in most stocks; and the deviations varied over years, with ~50% (19) of the stocks showing increasingly less differences in the last decade (Figure S8).The observed age distributions of many of these ( 8) stocks also differed from theoretical stable age distributions and were skewed toward younger ages (Figure S9) owing to disproportionately higher annual mortality rates of older fish by fishing (Figure S10).
Regionwide analyses detected declining trends in per capita recruitment success rate in both recovering and non-recovering stocks (ΔAIC = 2.1 and 2.7, Figure 4a), with high among-stock variances ( 2 stock = 0.31 and 0.48, Figure 4f-h).Ecoregion-specific random intercepts also revealed generally higher year-to-year fluctuations in the non-recovering than in the recovering stocks (Figure S6).After accounting for variations among stocks and ecoregions (and years within each ecoregion) the recruitment rates of both stock groups declined with increasing relative spawner biomass (Figure 4c; Table 2), which was more supported by data than absolute biomass (Table S2).But the temporal patterns of density dependency differed between the groups; the density-dependent effect in the recovering stocks increased as the biomass rose, whereas that in the non-recovering stocks emerged in the last two decades but declined over time (Figure 4c).Likewise, the contribution of spawner age structure to recruitment variability differed between the groups; the relationship between recruitment success and mean age in the recovering stocks shifted from positive to negative in the early 1990s, whereas that in the non-recovering stocks remained largely positive (Figure 4b).An indirect positive fishing effect on recruitment rate was detected in the recovering stocks until the late 1980s but gradually diminished since (Figure 4d).The models for both stock groups showed weak support for climate effects (Table S2).Of the climate indices tested for the non-recovering stocks (whose model had slightly more support than the recovering stocks) the model with ocean heat content anomaly in the year of spawning showed shifting (positive to negative) climate effects in the recent decades (Figure 4e).
Ecoregion-scale analyses on recruitment success revealed patterns largely consistent with the regional analyses with some systemspecific deviations (Figure S11).In ecoregions with non-recovering stocks like the Baltic, Celtic, and North Seas, the density-dependent effect on recruitment on average became weakened over time; in contrast, the effect in the Norwegian-Barents Seas (with only recovering stocks) became strengthened (Figure S11c).Correlations with spawner age structure metrics were detected in all but the Bay of Biscay and Faroes (Figure S11b); in the Baltic Sea, a positive correlation with age diversity became detectable in the early 1990s when fishing pressures began declining (Figure 2).The time-varying effect of spawner age was especially pronounced in the Norwegian-Barents Seas stocks, showing a shift from positive to negative in the 1990s (Figure S11b).An indirect fishing effect was (weakly) supported in the North Sea stocks only (Table S3; Figure S11d).Climate anomalies (in the year of spawning) contributed to variability in recruitment success in all ecoregions except the North Sea but the indices with the most statistical support varied among ecoregions (Figure 5; Table S3) despite high correlations among climate indices (Figure S12).The ways climate anomalies contributed to recruitment variability also differed (Figure 5).In the Baltic Sea and Faroes correlations with climate anomalies remained largely positive during the past few decades (Figure 5a,d).In contrast, the positive effect in the Norwegian-Barents Seas began diminishing (Figure 5f).In the Bay of Biscay and Icelandic waters the climate effect shifted from positive to negative (Figure 5b,e).
Trend analyses revealed consistent temporal patterns in elasticities to vital rates for the recovering and non-recovering stocks.
Both showed declining elasticities to (post-recruitment) juvenile survival rate (ΔAIC = 15.6 and 18.9) and increasing elasticities to adult survival rate over time (ΔAIC = 17.0 and 6.7; Figure S13).fluctuated over years, but on average elasticities to adult survival showed higher deviations from the regional trends than the other vital rates (Figure S13).Covariates selected for each type in the models (Table 2; Table S2) are indicated in the y-axis labels.Relative spawner biomass is defined as estimated spawner biomass relative to management reference points (B pa , Table 1) defined by the International Council of Exploration of the Sea (ICES).Color gradients and contour lines indicate the partial effects of covariates estimated (as tensor interaction terms) in the models (only the covariate types that are supported by data are shown).
Warmer colors indicate higher values, and the ranges of the partial effects differ among covariates and are indicated as values on the contour lines.Gray circles indicate backtransformed input data.(f-h) Random effects estimates (±95% confidence intervals) for stock and ecoregion for the recovering and non-recovering stocks (ecoregion-specific estimates for the non-recovering stocks are negligible; random effect estimates for year nested within ecoregion are given in Figure S6b).[Colour figure can be viewed at wileyonlinelibrary.com] indicated that more diverse ages increasingly contributed to elasticities in more recent years, whereas the contributions of mean ages gradually increased in the non-recovering stocks (Figure 6b).
Elasticities to recruitment declined with higher spawner biomass in both stock groups, but the density-dependent effect in the nonrecovering stocks diminished over time (Figure 6c).The model with an indirect fishing effect in both stock groups showed a declining contribution to variations in the elasticities after the reductions in TA B L E 2 Parameter estimates and test statistics from regionwide analyses using the generalized additive mixed effects models selected for the northeast Atlantic fish stocks with per capita recruitment success (RS total ) as a response variable (Table S2).fishing pressure (Figure 6d).In contrast, the models for both stock groups indicated an increasingly more positive response to climate anomalies (the AMO index and winter sea temperature anomaly) over time (Figure 6e).
Ecoregion-scale analyses on elasticities to recruitment success detected declining trends in all but the Bay of Biscay and the Norwegian-Barents Seas (Figure S14a), which showed little change.The analyses also indicated spatial heterogeneity in the ways time-varying intrinsic and extrinsic drivers contributed to variations in the elasticities.For example, shifting contributions of spawner age structure was detected in most ecoregions, showing both diminishing (e.g., Faroes and the North Sea) and increasing (e.g., the Baltic Sea and Icelandic waters) trends (Figure S14b).
Negative correlations with spawner biomass were consistently detected except in the Baltic Sea (Figure S14c).This maternal effect increased over time except in the Norwegian-Barents Seas, which showed the opposite trend (Figure S14c).The models with an indirect fishing effect had statistical support for Faroes, and the Celtic, North, and Norwegian-Barents Seas; the fishing effect declined and/or shifted over time except for the Norwegian-Barents Sea, which showed an increasing contribution of fishing since the 1990s (Figure S14d).3; Table S2) are indicated in the y-axis labels.Relative fishing pressure and relative spawner biomass are defined as estimated fishing mortality rate and spawner biomass relative to management reference points (F msy and B pa , Table 1) defined by the International Council of Exploration of the Sea (ICES).Color gradients and contour lines indicate the partial effects of covariates estimated (as tensor interaction terms) in the models.
Warmer colors indicate higher values, and the ranges of the partial effects differ among covariates and are indicated as values on the contour lines.Gray circles indicate backtransformed input data.(f-h) Random effects estimates (±95% confidence intervals) for stock and ecoregion for the recovering and non-recovering stocks (ecoregion-specific parameters for the non-recovering stocks are negligible; random effect estimates of year nested within ecoregion are given in Figure S6c).

| Transient dynamics of populations recovering from intense fishing pressures
Fishing-induced changes in the demography and reproductive capacity of an exploited population can be reversible when management measures are effectively applied to regulate fishing pressure.This pattern is displayed especially in the Norwegian-Barents Seas, where fishing pressures were reduced by more than 60% (compared to the regional average of ~50%) since the early 1990s, followed by sharp rises in the average age and biomass of reproductive fish.But when a population experiences long-term intensive harvesting, recovery time can become longer (Hutchings, 2015) and more uncertain (Neubauer et al., 2013), with the population remaining in long transient states (Hastings et al., 2018).Although fishing pressures on many northeast TA B L E 3 Parameter estimates and test statistics from regionwide analyses using the generalized additive mixed effects models selected for the northeast Atlantic fish stocks with elasticity to recruitment success (e RS ) as a response variable (Table S2).which could further diminish fishery yields (Hilborn et al., 2020).
Removal of some biomass from a population by harvesting can promote its productivity (Neubauer et al., 2013) by attenuating densitydependent effects on early life stage processes (competition for food and space, e.g., Anderson et al., 2008).Intense fishing pressures before the implementation of stricter management measures-mediated through reproductive fish-may have led to higher recruitment rates in depleted stocks in the northeast Atlantic Ocean.Fishing not only truncates age structure but also induces changes in life-history traits that may influence recovery (Hutchings & Kuparinen, 2020;Neubauer et al., 2013).The generation times of these depleted northeast Atlantic stocks were shortened (by 14-50%), which is often associated with high year-to-year variability in recruitment (Minto et al., 2008) and high average recruitment rates (Bjørkvoll et al., 2012).Higher growth rates of depleted (and age-truncated) fish populations also likely benefited during the early stages of recovery as more reproductive fish contributed to spawning events, thereby increasing the probability of recruitment success.These initial benefits, however, diminished over time following further reductions in fishing pressure.Despite the signs of recovery in reproductive fish, per capita recruitment success rates on average declined in the northeast Atlantic stocks, as also reported in other Atlantic ecoregions (Britten et al., 2016).But these apparent patterns were shaped by different processes among stocks.In the recovering stocks density dependency remained present in recruitment success throughout the past decades and its strength increased as reproductive fish biomass rose, likely slowing recoveries as fishing pressures declined.
In contrast, in the non-recovering stocks density dependency in recruitment success emerged only after the reductions in fishing pressure but its strength continued to decline with reproductive fish biomass, suggesting contributions of other (possibly extrinsic) drivers.
Density dependency in early life also perhaps explains discontin- Greater sensitivities to adult survival also may make harvested populations more manageable.Nonlinear responses to age structure in reproductive fish are, however, also evident in sensitivity to recruitment variability.Although the buffering effect of increasingly age-diverse reproductive fish emerged in many ecoregions especially in the early stages of recovery, this relationship also shifted as the recoveries further progressed (as seen in the Norwegian-Barents Seas fishes), likely responding to other changes in intrinsic or extrinsic drivers (such as catch limits being raised as stock biomass rebuilds).

| Climate regulation of demographic variations
Fish populations experiencing persistent selective depletions of older adults often become sensitive to variations in early life-history stage processes that reflect signals in environmental conditions like climate variability (Ottersen et al., 2006;Rouyer et al., 2011).Warmer than average climate conditions in the northeast Atlantic Ocean increasingly contributed to observed patterns in fish recruitment success, reflecting shifting climate effects on recruitment processes.For example, stratification becomes intensified as climate warms, limiting nutrient movement to deeper waters and food accessibility for the early life of demersal and benthic species (Moore et al., 2018).At the region scale, however, the models captured only weak responses in recruitment to climate signals likely because of heterogeneous warming patterns among ecoregions, as also seen in other parts of the world (Barange et al., 2018).Region-scale analyses may mask differential effects of local climate and oceanographic features (e.g., climate-induced shifts in ocean heat content and primary productivity) that marine species are adapted to and experience (Conover, 1992).
Ecoregion-scale variability in recruitment success may reflect how fish reproductive behavior responds to regional climate mediated by local habitat conditions, as hypothesized in past research (Myers, 1998;Ottersen & Stenseth, 2001).For stocks in high latitude systems like the Norwegian-Barents Seas and Icelandic waters, for example, climate conditions in winter-spring (December-May) that regulate phonologies like timings of spawning migration (Sundby & Nakken, 2008) contributed to recruitment variability.Warming-induced reductions in sea ice extent in these ecoregions promoted ecosystem productivity in the late 1990s and early 2000s and may have benefited upper trophic animals like fish (Stenevik & Sundby, 2007), likely through changes in growth rate and life stage duration of young fish (Houde, 2016).When spawning times coincide with the period of high prey abundance, pre-recruit fish can grow more and experience less size-dependent mortalities (like predation), contributing to greater recruitment success (Peck et al., 2012).In contrast, in semi-enclosed systems like the Baltic Sea warming-induced changes in nearshore habitat conditions in summer (June-August; lower salinity caused by higher precipitation and freshwater inflows, e.g., MacKenzie et al., 2007) may have played a greater role in regulating recruitment success through biophysical processes during the early years of life.
Demographic responses by fish populations to changing climate conditions also may vary discontinuously (Brander, 2007).Local and regional climate conditions (both favorable and unfavorable) that early life-history stage fish experience can persist through time (Ottersen et al., 2001) and may emerge as delayed effects on transient population dynamics.In some northeast Atlantic ecoregions the ways that climate affects fish recruitment success shifted from positive to negative and may have slowed their recoveries.These shifts suggest that even in the ecoregions where warming climates promoted productivity in the early 20th century, fish began to experience adverse (direct and indirect) effects.Mobile species can adjust their movement (deeper and poleward) to mitigate physiological stress from warming upper water columns (Engelhard et al., 2014;Pinsky et al., 2013).But asynchronous shifts in spatial distribution by spawning fish and seasonal events in lower trophic groups like spring plankton blooms (phenological mismatches) also could trigger a series of events unfavorable for early life fish (Asch & Erisman, 2018), shrinking the probability of recruitment success.
Climate variability that regulates local ocean conditions may propagate through vital rates to shape transient population dynamics of marine species (Jenouvrier et al., 2018).Depleted populations in particular may have heightened sensitivities to adverse environmental conditions (Hidalgo et al., 2011;Rouyer et al., 2011), which are becoming increasingly more variable in recent decades (Rathore et al., 2020).Although the sensitivities of transient growth to recruitment variability in the northeast Atlantic fish stocks declined as they recovered, the analyses reveal that increasingly warmer than average climates began altering how recruitment variability propagates through transient dynamics.For example, in high-latitude systems like Faroes, where fish stocks became less sensitive to recruitment variability under a warmer climate, this stabilizing effect gradually diminished (and shifted to the opposite) as the climate further warmed.
In contrast, in the Baltic Sea, where fish stocks became more sensitive to recruitment variability under a warmer climate, this amplifying effect became strengthened as the climate further warmed.Climate effects on transient dynamics may, however, also depend on population size and demographic structure.In the Norwegian-Barents Seas, where fish stocks experienced higher reductions in fishing pressure early and thus higher biomass rebuilding rates than in the other ecoregions, the stocks responded to recruitment variability independently of climate conditions.Emerging non-negligible changes in demographic responses to rising climate variability nevertheless may influence the efficacy of fisheries management.Although historical climate conditions promoted the productivity of fish stocks in the northeast Atlantic Ocean, some stocks in this region and elsewhere in the Northern Hemisphere already are experiencing adverse effects of changing climates (Free et al., 2019;Pershing et al., 2015).Because many biological responses to climate-driven changes in surface temperature and other ocean properties are often nonlinear, we would expect continued shifts in how climate influences marine resource populations and management (Britten et al., 2016).Accounting for such climate-induced changes thus may improve the quality of assessments and the performance of management measures (Bahri et al. 2021;Rouyer et al., 2011).
Exploited marine fish species are an integral part of an ecosystem and changes in their abundances and demographic structures may propagate through community and ecosystem dynamics (Jennings & Kaiser, 1998), triggering a series of changes in ecologically connected species (Goto, Devine, et al., 2022;Pérez-Rodríguez et al., 2022).As

CO N FLI C T O F I NTE R E S T S TATE M E NT
The author declares no conflict of interest.
In the following, I describe (a) the development of matrix population models to quantify transient population growth rates and elasticities through transient life table response experiments, (b) trend analyses to characterize the dynamics of reproductive fish biomass and age structure, transient growth, and fishing pressure, followed by analyses to test for direct fishing effects, and (c) evaluation of the relative F I G U R E 1 Study systems and climate conditions in the northeast Atlantic Ocean.(a) Geographic locations of the ecoregions, which are defined by the International Council of Exploration of the Sea (ICES): the Norwegian and Barents Seas (62°0'N-82°0'N at 11°0'W-68°30′ E), the Greater North Sea (56°0'N at 2°55′ E), the Baltic Sea (53°352032′N-65°54′N at 11°57′ E-30°14′ E), Faroes (60°0′N-63°0′N at 15°0′W-4°0′W), Icelandic waters (64°38′N at 18°55′W), the Bay of Biscay and the Iberian Coast (42°39′N at 8°53′W), and the Celtic Seas (48°0′N-60°30′N at 18°0′W-1°21′W).Map lines delineate study areas and do not necessarily depict accepted national boundaries.(b) Regional climate indices (the North Atlantic Oscillation (NAO) index, the Atlantic Meridional Oscillation index, and upper (0-400 m) ocean heat content (OHC) anomaly) during 1946-2016.(c) Seasonal sea surface temperature anomalies in each ecoregion.The time series of temperature anomalies in each ecoregion are truncated to show only the years that the fish stock data and assessments included in the study are available.[Colour figure can be viewed at wileyonlinelibrary.com] contributions of time-varying intrinsic (reproductive fish biomass and age structure) and extrinsic (fishing pressure and climate) drivers to variations in recruitment success and in elasticity to recruitment success.
The final dataset consists of 38 stocks comprising 11 species (cod Gadus morhua, haddock Melanogrammus aeglefinus, saithe Pollachius virens, herring Clupea harengus, plaice Pleuronectes platessa, turbot Scophthalmus maximus, sole Solea solea, whiting Merlangius merlangus, sprat Sprattus sprattus, megrim Lepidorhombus spp., and four-spot megrim Lepidorhombus boscii) across seven ecoregions in the northeast Atlantic Ocean defined by the ICES: the Norwegian and Barents Seas (combined as a single ecoregion for analysis), ): the North Atlantic Oscillation (NAO) and Atlantic Multidecadal Oscillation (AMO) indices (both computed from the Hadley Centre Sea Ice and Sea Surface Temperature dataset [HadISST1, https://www.metoffice.gov.uk] and the extended reconstructed SST [ERSSTv5, NOAA Physical Sciences Laboratory; Stock ID indicates a unique ID for each stock assigned by the International Council of Exploration of the Sea (ICES); period and years indicate the years and the total number of years that stock data are available; age indicates the age classes included in the assessments; F range indicates the age range used to compute mean fishing mortality rates for assessment; F msy and B pa indicate the management

c
Stock information is taken from ICES Scientific Reports: (a) ICES (2017a), (b) ICES (2017e), (c) ICES (2017b), (d) ICES (2017c) (e) ICES (2017f), (f) ICES (2017d).TA B L E 1 (Continued)(non-fishing like predation and starvation) mortality rates for a-yearolds in year y.Egg production rates are often not estimated in the assessments of marine fish stocks, and here I estimated recruitment success (integrating over reproduction and early life-history stage survival) rate per spawner (reproductive fish) for a-year-olds in year y (RS a,y;Rouyer et al., 2011) as a proxy for per capita fecundity rate as where R y is the recruit number in year y, m a,y−a R and N a,y−a R are the maturity rate (fraction of juveniles becoming adults) and abundance of a-year-olds in year y-a R (spawning year), a R is the age of recruits, which accounts for variations in age at recruitment (0-3 years old) among stocks reported in the assessments.a min and a max are the youngest and oldest ages in the stock.I then numerically estimated deterministic transient growth rates (λ trans ) by projecting age-specific abundances using the annual state transition matrices and the annual estimates of age-specific numbers taken from the assessments as 2.4.3 | Covariates of elasticity to recruitment successThe transient dynamics of depleted fish stocks with age truncation often become recruitment-dependent and thus sensitive to variations in climate and habitat conditions.As the stocks recover (increases in abundance and age diversity), the transient dynamics are expected to become less sensitive to recruitment variability.To test for the time-varying contributions of the extrinsic and intrinsic factors to variations in elasticity to RS (e rec , Figure2a), the following model was fitted to the data, assuming e rec follows a Beta distribution with expected values μ (0 < μ < 1) and a logit link function ase rec,ijk ∼ Beta ijk , ф E e rec,ijk = ijk var e rec,ijk = ijk 1− ijk ∕ (ф + 1) F I G U R E 2Factors affecting the transient dynamics of northeast Atlantic fish stocks.(a) Schematic diagram of hypothesized mechanisms involving intrinsic (spawner biomass and age structure) and extrinsic (fishing and climate) factors tested for demographic response variables (recruitment success and elasticity to recruitment variability) in the study.(b) Temporal trends in relative fishing pressure, relative spawner biomass, mean spawner age, and spawner age diversity of the recovering and non-recovering stocks during 1946-2016.Vertical red dashed lines indicate the reference year (2002) used to define the recovery status of fish stocks based on relative spawner biomass in analysis.(c) Relationships between relative fishing pressure and relative spawner biomass or age structure metrics.Solid lines and ribbons indicate regional average trends in time (b) or relative fishing pressure (c) and 95% confidence intervals estimated by generalized additive mixed effects modeling with year (b) or relative fishing pressure (c) as a fixed effect, and ecoregion, year within ecoregion, and stock as random effects.Circles indicate partial residuals estimated from the models.Relative fishing pressure and relative spawner biomass are defined as estimated fishing mortality rate and spawner biomass relative to management reference points (F msy and B pa , Table 1) defined by the International Council of Exploration of the Sea (ICES).[Colour figure can be viewed at wileyonlinelibrary.com] sion 4.2.1, R Development Core Team, 2022) with the R packages gamm4 (v.0.2-6,Wood & Scheipl, 2020) and mgcv (v.1.8-40,Wood, 2022) for model fitting, mgcv and gratia (v.0.7.3, Simpson & Singmann, 2022) for model diagnostics, and popbio (v.2.7, Stubben et al., 2020) and popReconstruct (1.0-7, Wheldon, 2019) for computations of demographic parameters and projections of population dynamics.
recovering and the non-recovering fish stocks reflect changes in fishing pressure, (b) how intrinsic (spawner biomass and age structure) and extrinsic (fishing and climate) drivers contribute to variations in recruitment success, and (c) how the intrinsic and extrinsic drivers influence the sensitivity of transient dynamics to recruitment variability.Overall, the diverging patterns in spawner biomass between the recovering and the non-recovering stocks also emerged in age structure and transient population dynamics, reflecting differential reductions in fishing pressures the stocks experienced.The differences in spawner biomass and age structure between the groups were further reflected in variations in recruitment success and elasticity to recruitment success, indicating divergence in the strength of density dependency and the contribution of older spawners.
Figure 6a) with less among-stock variances than recruitment success (Figure 6f-h).Ecoregion-specific variances in all the models Regionwide analyses on elasticities to recruitment success largely reflected the time-varying responses in recruitment to intrinsic and extrinsic drivers.For the recovering stocks the model F I G U R E 4 Regional trends in and covariates of per capita recruitment success (fecundity) of the recovering and non-recovering northeast Atlantic fish stocks during 1946-2016.(a) Temporal trends.Solid lines and ribbons indicate temporal trends and 95% confidence intervals estimated by generalized additive mixed effects modeling with year as a fixed effect, and ecoregion, year within ecoregion, and stock as random effects.Red circles indicate partial residuals estimated by the models.Vertical red dashed lines indicate the reference year (2002) used to define the recovery status of fish stocks based on relative spawner biomass in analysis.(b-e) Time-varying effects of spawner age structure, spawner biomass, fishing pressure, and climate condition.

FIGURE 6
FIGURE 6 Regional trends in and covariates of elasticity to per capita recruitment success of the recovering and non-recovering northeast Atlantic fish stocks during 1946-2016.(a) Temporal trends.Solid lines and ribbons indicate temporal trends and 95% confidence intervals estimated by generalized additive mixed effects modeling with year as a fixed effect, and ecoregion, year within ecoregion, and stock as random effects.Circles indicate partial residuals estimated by the models.Vertical red dashed lines indicate the reference year (2002) used to define the recovery status of fish stocks based on relative spawner biomass in analysis.(b-e) Time-varying effects of spawner age structure, spawner biomass, fishing pressure, and climate condition.Covariates selected for each type in the models (Table3; TableS2) are indicated in the y-axis labels.Relative fishing pressure and relative spawner biomass are defined as estimated fishing mortality rate and spawner biomass relative to management reference points (F msy and B pa , Table1) defined by the International Council of Exploration of the Sea (ICES).Color gradients and contour lines indicate the partial effects of covariates estimated (as tensor interaction terms) in the models.Warmer colors indicate higher values, and the ranges of the partial effects differ among covariates and are indicated as values on the contour lines.Gray circles indicate backtransformed input data.(f-h) Random effects estimates (±95% confidence intervals) for stock and ecoregion for the recovering and non-recovering stocks (ecoregion-specific parameters for the non-recovering stocks are negligible; random effect estimates of year nested within ecoregion are given in FigureS6c).[Colour figure can be viewed at wileyonlinelibrary.com] FIGURE 6 Regional trends in and covariates of elasticity to per capita recruitment success of the recovering and non-recovering northeast Atlantic fish stocks during 1946-2016.(a) Temporal trends.Solid lines and ribbons indicate temporal trends and 95% confidence intervals estimated by generalized additive mixed effects modeling with year as a fixed effect, and ecoregion, year within ecoregion, and stock as random effects.Circles indicate partial residuals estimated by the models.Vertical red dashed lines indicate the reference year (2002) used to define the recovery status of fish stocks based on relative spawner biomass in analysis.(b-e) Time-varying effects of spawner age structure, spawner biomass, fishing pressure, and climate condition.Covariates selected for each type in the models (Table3; TableS2) are indicated in the y-axis labels.Relative fishing pressure and relative spawner biomass are defined as estimated fishing mortality rate and spawner biomass relative to management reference points (F msy and B pa , Table1) defined by the International Council of Exploration of the Sea (ICES).Color gradients and contour lines indicate the partial effects of covariates estimated (as tensor interaction terms) in the models.Warmer colors indicate higher values, and the ranges of the partial effects differ among covariates and are indicated as values on the contour lines.Gray circles indicate backtransformed input data.(f-h) Random effects estimates (±95% confidence intervals) for stock and ecoregion for the recovering and non-recovering stocks (ecoregion-specific parameters for the non-recovering stocks are negligible; random effect estimates of year nested within ecoregion are given in FigureS6c).[Colour figure can be viewed at wileyonlinelibrary.com] Atlantic fish stocks precipitously declined through stricter management measures applied since the late 1990s, demographic restorations and stock recoveries on average progressed gradually.Reducing fishing pressures can allow more adults to survive and contribute to the production of new recruits but the restoration of the demographic structure takes place over generations and may take longer (a decade or more) for severely depleted stocks.Furthermore, although the declining trends in reproductive fish biomass halted in the early 2000s in most ecoregions, reduction levels in fishing pressure varied among ecoregions and stocks and were reflected in recovery patterns.When reductions in fishing pressure were insufficient or delayed, the declining trends in reproductive capacity and deterioration of demographic structure in some stocks continued.Slow recoveries of severely depleted fish stocks are reportedly widespread(Hutchings, 2000).Regionwide analyses on the non-recovering stocks in the northeast Atlantic reveal that sustaining high fishing pressures for even a few more years can result in vastly different recovery patterns, with further age truncations and reductions in reproductive capacity,

4. 3 |
Transients as an integrated measure for assessing managed populationsAnalysis of transient population dynamics can reveal how timevarying demographic structure and vital rates contribute to recovery patterns of managed resource populations when also exposed to low-frequency, nonstationary drivers like climate.Declined sensitivities to highly variable demographic parameters like recruitment, as displayed by intensely managed northeast Atlantic fish stocks in this study, can provide a buffer against environmental variability.Increasing sensitivities to adult survival further suggest that the stocks also are expected to become more responsive to management actions through regulation of fishing pressure (enforcing positive feedback), and future management plans would benefit from maintaining healthy (age-diverse) reproductive populations of exploited species(Pinsky & Byler, 2015).Combined with other management measures (like use of marine protected area and fish size limit,White et al., 2013;Hixon et al., 2014), restoring and maintaining the demographic structure of an exploited species thus may enhance resilience against destabilizing effects of environmental variability amplified by climate change.
adult survival improves in historically overharvested populations, ecological processes including resource competition and predation also may indirectly modify vital rates and demographic structures of non-target species (van Gemert & Andersen, 2018).Improving our understanding of how warming and more variable climates continue to reshape the transient dynamics (through changes in life-history traits and demographics) of marine populations would help develop management plans to mitigate the risk of overexploitation and safeguard ecosystem services against rising environmental variability.ACK N OWLED G M ENTS This work was partially supported by the Institute of Marine Research's Reduced Uncertainty in Stock Assessment (REDUS) project.Some figures use images from the IAN Symbols, courtesy of the Integration and Application Network, University of Maryland Center for Environmental Science (ian.umces.edu/symbols/).
stocks I used mean age and age diversity of spawners as indices for demographic structure.Mean age (SA mean ) is the weighted (by abundance) average of spawner age as The Shannon diversity index was used as an indicator of age diversity (SA diversity ) and computed as SA diversity,y = − the assessments) and propagated the estimation uncertainties through the computation of state transition and transient population growth rates (see FigureS4for more detail).This step was then followed by multilevel model-fitting with transient population growth as a response variable and year as a covariate (as described above) applied to each simulated dataset.
To address this issue I did sensitivity analyses by Monte Carlo simulations to illustrate how uncertainties in the estimates can propagate through matrix population and multilevel modeling analyses.I first generated 1000 simulated datasets with Gaussian noise in recruit numbers and fishing mortality rates using the stock-and year-specific coefficients of variation derived from the assessments (uncertainty in age-specific parameters was not consistently reported in jik ,τ, β 0 , stock , ecoregion , and year are the same as above.f 1 is a smoother function for relative fishing pressure.The contributions of fishing, management action, and climate to the transient dynamics of fish stocks with variable age structures were evaluated in two steps.I first tested for time-varying contributions of these extrinsic factors to variations in vital rates and then evaluated how these factors influenced the elasticities of transient population growth:2.4.2 | Covariates of recruitment successBecause natural mortality rates were assumed time invariant in most of the assessments, variations in survival rates are in effect driven by fishing.Thus, RS was selected as a demographic response variable (which integrates probabilities of spawning and early life-history stage survival) for analysis.Because RS can vary with variations in spawner age structure and biomass, these intrinsic factors were also tested.To test for the time-varying contributions of the extrinsic and intrinsic factors (Figure Note: EDF, F, 95% CI, t, and SD indicate the effective degrees of freedom, F-statistics, 95% confidence intervals, t-values, and standard deviations.All time (year)-dependent smoother terms for covariates are modeled as tensor product construction (using the te function in the R package mgcv).