A century of fish growth in relation to climate change, population dynamics and exploitation

Marine ecosystems, particularly in high‐latitude regions such as the Arctic, have been significantly affected by human activities and contributions to climate change. Evaluating how fish populations responded to past changes in their environment is helpful for evaluating their future patterns, but is often hindered by the lack of long‐term biological data available. Using otolith increments of Northeast Arctic cod (Gadus morhua) as a proxy for individual growth, we developed a century‐scale biochronology (1924–2014) based on the measurements of 3,894 fish, which revealed significant variations in cod growth over the last 91 years. We combined mixed‐effect modeling and path analysis to relate these growth variations to selected climate, population and fishing‐related factors. Cod growth was negatively related to cod population size and positively related to capelin population size, one of the most important prey items. This suggests that density‐dependent effects are the main source of growth variability due to competition for resources and cannibalism. Growth was also positively correlated with warming sea temperatures but negatively correlated with the Atlantic Multidecadal Oscillation, suggesting contrasting effects of climate warming at different spatial scales. Fishing pressure had a significant but weak negative direct impact on growth. Additionally, path analysis revealed that the selected growth factors were interrelated. Capelin biomass was positively related to sea temperature and negatively influenced by herring biomass, while cod biomass was mainly driven by fishing mortality. Together, these results give a better understanding of how multiple interacting factors have shaped cod growth throughout a century, both directly and indirectly.


| INTRODUC TI ON
Climate change is affecting marine ecosystems worldwide and there have been measurable and increasing consequences to populations over the past century, affecting individual physiology and survival as well as population dynamics, species distribution, productivity or ecological diversity (Cheung et al., 2013;García-Reyes et al., 2015;Hoegh-Guldberg & Bruno, 2010). Temperatures are predicted to increase in all projected scenarios (Stocker et al., 2014) and ocean warming has consequently been at the center of research.
However, its future consequences for fish populations remain difficult to predict given the scarcity of long-term biological chronologies for certain marine ecosystems or communities (Poloczanska et al., 2013). Additionally, many fish populations are subject to other significant stressors, especially fisheries exploitation. Selective harvesting can gradually reduce population complexity (Hilborn, Quinn, Schindler, & Rogers, 2003), resulting in changes in life history (Jørgensen, 1990;Law, 2000) and demography (Ottersen, Hjermann, & Stenseth, 2006). These alterations, in turn, erode population resilience to environmental changes and magnify their impacts (Hidalgo et al., 2011;Morrongiello, Sweetman, & Thresher, 2019;Planque et al., 2010). In this context of significant stressors, there is a need to evaluate the past variability of fish populations to identify, quantify and contextualize the causes of observed changes and forecast their future impacts accordingly (Reid & Ogden, 2006).
Fish somatic growth is an ideal proxy which responds to both environmental changes and exploitation pressure. Spatial or temporal changes in growth have been associated with different environmental factors such as temperature (Campana, Mohn, Smith, & Chouinard, 1995) or prey availability (Graeb, Dettmers, Wahl, & Cáceres, 2004) but also to population-specific effects such as density dependence (Lorenzen & Enberg, 2002) or harvesting (Enberg, Jørgensen, Dunlop, Heino, & Dieckmann, 2009). Small changes in growth rates within a population can not only influence individual fitness but also population mortality, productivity and reproductive success (Audzijonyte, Kuparinen, Gorton, & Fulton, 2013;Hixon, Johnson, & Sogard, 2014;Lorenzen, 2016). As a result, long-term growth chronologies can be valuable indicators of population response to the combined effects of environment and exploitation (Morrongiello et al., 2019) and associated changes in fitness over time (Morrongiello, Walsh, Gray, Stocks, & Crook, 2014).
Fish growth data are often derived directly from size-at-age data as a population index, using growth functions to estimate theoretical maximum size and the rate at which it is reached (e.g., Pilling, Kirkwood, & Walker, 2002). However, growth processes are controlled by both intrinsic (individual-specific) and extrinsic (environment-specific) factors affecting the accessibility and allocation of necessary resources. In particular, resources acquired by an individual are distributed between basal metabolism, structural growth and reproduction (Enberg et al., 2012), which means growth rates are directly dependent on multiple processes fluctuating throughout an individual's life. Methods relying solely on population data such as age-length keys may therefore only represent the "final" state of fish individuals and fail to capture the importance of inter-and intra-individual variability in growth trajectories. Linking individual growth histories from different cohorts can instead provide a better understanding of the different sources of growth variation and aid in interpreting historical growth patterns. Individual growth patterns can then be applied to population-level questions through biochronological reconstructions. Such growth chronologies based on calcified tissues have been successfully developed for both bivalves and fish, and show promising potential to study long-term, ecologically relevant data at fine-scale resolution (Morrongiello, Thresher, & Smith, 2012). Incremental growth of these hard structures is typically closely related to somatic growth (Black et al., 2019;Doubleday et al., 2015) and often identifiable at fixed temporal scales (i.e., daily, seasonal or yearly), thus providing continuous time-resolved growth histories of individuals and populations (Black, von Biela, Zimmerman, & Brown, 2013;Morrongiello, Crook, King, Ramsey, & Brown, 2011).
In particular, large numbers of fish otoliths are collected worldwide every year (Campana, 2001), making it a widespread and easily accessible resource for developing long-term individual-based growth chronologies across a wide range of environments and fish taxa (Morrongiello et al., 2012). Otolith-derived biochronologies have been developed in diverse contexts, with a strong focus on the influence of climatic factors on growth. For example, classical detrending approaches directly derived from tree rings analysis (dendrochronology) have been successfully applied to fish studies to reconstruct population-level variability in growth based on otolith increment widths (Black et al., 2013;Matta, Black, & Wilderbuer, 2010;Ong et al., 2018). However, more modern statistical approaches such as mixed-effect modeling allow for in-depth investigation of how fish growth varies in response to different factors (Weisberg, Spangler, & Richmond, 2010). A key advantage of mixed-effects models is that they do not detrend or standardize growth time series to maximize relationships at the population level. Rather, they partition and investigate both biological and environmental effects simultaneously by utilizing all available biological information and accounting for the hierarchical structure of the data within individuals, years and cohort (year-classes, i.e., year of hatching; Morrongiello & Thresher, 2015). This provides a broader and more comprehensive understanding of how fish respond to particular factors of interest, as well as being more appropriate for the study of species where growth variations might be subject to high inter-individual variability (Black et al., 2019).
The Barents Sea is a highly productive shallow area situated to the north of Norway, bordering the Norwegian Sea to the west and the Arctic Ocean to the north (Jakobsen & Ozhigin, 2011). Described as the "Arctic warming hotspot" (Lind, Ingvaldsen, & Furevik, 2018), this region at the interface between the Atlantic and the Arctic is an ideal system to study the effects of global warming on ecosystems due to its fast changing hydrography. It is also home to one of the largest populations of the commercially important Atlantic cod (Gadus morhua; ICES, 2018), the Northeast Arctic (NEA) cod, a large migratory population with a long history of human exploitation. NEA cod lives and feeds in the Barents Sea, where it is considered an apex predator with an essential role in the food-web dynamics, notably through its high predation on small pelagic fish such as capelin (Holt, Bogstad, Durant, Dolgov, & Ottersen, 2019). Every winter it migrates southward to spawn around the Lofoten archipelago (Höffle et al., 2014), where the main fishery has also historically taken place (Sundby & Nakken, 2008). Spawning has recently been taking place north of the traditional area along the north-western Norwegian coast, suggesting a possible reaction to changes in the environment and demography (Opdal & Jørgensen, 2015). Multiple studies have highlighted the significant influence of environmental changes and exploitation on the biology and demography of cod across the North Atlantic (e.g., Brander, 2000;Eero, MacKenzie, Köster, & Gislason, 2011;Mieszkowska, Genner, Hawkins, & Sims, 2009). NEA cod is therefore a particularly attractive population to study the long-term impacts of climate change and exploitation, especially in light of the multiple collapses that occurred in the western Atlantic cod stocks since the 1990s (Bavington, 2011;Pershing et al., 2015).
In this study, we used a large archive of otoliths collected from commercial catches and scientific surveys to develop a near century-long biochronology for the NEA cod. Using a hierarchical mixed-effects modeling approach, we related the observed growth variability to selected climatic, demographic and fishing factors which have been considered to have an important influence on cod growth. We hypothesized that warming conditions and increasing prey availability likely had a positive influence on cod growth. We also hypothesized that fishing can selectively remove either fast or slow growing fish from a population (Enberg et al., 2012;Sinclair, Swain, & Hanson, 2002) but also alter density-dependent processes that, in turn, influence growth (Lorenzen & Enberg, 2002). We then applied path analysis to test hypothesized causality between the different factors tested within the growth modeling and identify their respective direct and indirect effects on cod growth. The results provide a long-term and integrated perspective on the different factors that have influenced NEA cod growth during the 20th and beginning of the 21st centuries.

| Sample selection
An extensive archive of cod otoliths exists at the Institute of Marine Research in Bergen (Norway), where about one million Atlantic cod otoliths from both research surveys and commercial fishing have been archived since 1920. Age readers have routinely aged cod from otoliths since 1932, and location, date, fishing gear and biological information (length, weight and sex) has been collected for most of the fish. For this study, sagittal otoliths of adult fish were randomly selected for each year between 1933 and 2015. Following the multiple cohort approach advocated by Morrongiello et al. (2012), these subsamples consisted, whenever possible, of 50 fish per year.
Sampling of overlapping cohorts can better capture the range of growth phenotypes as well as within-population variations. Sample selection was limited to fish caught by bottom trawl, longline and seine, to avoid selectivity bias toward bigger, faster-growing fish in gillnet catches. To sample otoliths from sexually mature individuals only, our selection was limited to fish collected from the Lofoten spawning grounds and consisted of a mixture of survey and fisherycaught fish (Figure 1). Mean age-at-maturity has decreased from 10.8 to 6.1 since 1920 due to fishing selectivity (Jørgensen, 1990;Nash, Pilling, Kell, Schön, & Kjesbu, 2010), and as a result the archive for the recent decades comprises relatively few fish older than 9.
Conversely, the archive for the earliest decades has fewer samples available from fish caught during the spawning season, which limited the age classes available for collection of a reasonable sample depth.
We consequently collected mainly fish of age 8 or older to provide enough samples of mature fish in the more recent decades while trying to minimize the effects of faster growth and maturation in the earlier period. The resulting growth chronology was nonetheless potentially biased toward earlier maturing fish and faster juvenile growth for the earlier years. Due to the limited availability of older fish in 1988, 1990, 1991 and 1996, some samples of age 7 were included (respectively 22, 2, 1 and 15 otoliths). A total of 4,096 fish of age 7-21 were retrieved, of which 3,023 were age 8.

| Otolith processing and measurement
Otoliths were embedded in black epoxy (NM Laminering 275,Nils Malmgren AB with NM Svart Pasta pigment paste) to enhance contrast around the edges. They were sectioned transversely (~800 µm thickness) through the core, which was identified during embedding by a distinctive narrow and hourglass-shaped depression near the sulcus, on the proximal surface of the otolith. Images of the sectioned F I G U R E 1 Map of Northeast Atlantic area and Barents Sea ecoregion. Arrows represent the main circulation of water masses. Dark grey area represents the main Northeast Arctic cod spawning grounds along the Norwegian coast. Blue line represents the Kola transect of temperature measurements. Top left insert shows the location of the study area [Colour figure can be viewed at wileyonlinelibrary.com] otoliths were digitalized using a high-resolution (5,616 × 3,744 pixels) mounted DSLR system comprising a Canon EOS 5D Mark II body and a Canon MP-E 65 mm f/2.8 1-5× lens. The imaging process was standardized at 3× magnification, 1/50 s exposure, f/6.3 and ISO 100 to ensure consistency in the captured pictures. The otolith images were then enhanced in Adobe Photoshop CC 2019 using a standardized macro that converted them to greyscale to remove color aberrations, adjusted levels/brightness/contrast to enhance the transition between opaque and translucent zones and added a sharpening mask to increase clarity and readability.
An open-access, dedicated set of ObjectJ macros was developed for the software ImageJ (Schneider, Rasband, & Eliceiri, 2012) to annotate the otolith sections following a standardized workflow (Supporting Information B). Otolith radius is often linearly correlated with fish total length (Francis, 1990;Harvey, Loughlin, Perez, & Oxman, 2000), and a strong relation between otolith and fish growth has been experimentally confirmed in Atlantic cod (Hüssy & Mosegaard, 2004;Li et al., 2008). Changes in increment widths can therefore be used directly as a proxy for changes in somatic growth (Campana, 1990). Because our study used mostly fish of age 8 and older, we verified this relationship by including a wider range of age and length classes from Høie et al. (2009) in the correlation analysis ( Figure S2).
Northeast Arctic cod form an annual otolith growth increment comprising a translucent zone formed from December to April and an opaque zone formed from May to November (Høie et al., 2009), and increment growth can thus be used as a proxy for annual somatic growth. Increments widths (µm) were measured on the distal side of the otolith along a single linear axis. Because the quality of the core was variable depending on the precision of the sectioning, we marked the largest diameter of the first increment and defined our measurement axis as the line in the distal direction that intersected perpendicularly the maximum number of increments from the otolith edge to the intersection with the drawn diameter ( Figure 2). Each increment width was assigned to a year of formation by counting back from date of capture and accounting for marginal increment interpretation. These measurements then provided subsequent estimates of annual growth, yearly age, age-at-capture and year-class for each individual. One otolith from each pair had previously been broken for age estimation by experienced otolith readers, and the age provided was consequently used as verification. The high agreement (95.76%) and low coefficient of variation (3.70%) between both readings indicated no important differences between age estimates of sectioned and previously broken otoliths. Age-bias cross-tabulation (Table S1) and plot ( Figure S4) indicated that most of the discrepancies were mainly associated with a difference of 1 year (positive or negative). Measurements from aberrant, crystalline or damaged otolith sections were removed from the analysis. Additionally, otoliths from fish older than 12 were removed after initial model exploration, as the scarcity of data points for these age classes created convergence issues during modeling. A total of 3,894 samples were consequently retained in the final analysis ( Figure S1). Finally, both the core (growth between hatching and first winter ring) and the marginal (residual growth between last winter ring and outer edge) increments were excluded from the analysis as they did not comprise a whole year of growth. In total, the measurements of 28,504 increments ranging from age 2 to 12 were included in the analysis ( Figure S3).

| Growth factors
Possible sources of inter-annual growth variation were identified and selected for analysis (Table 1). Intrinsic variables representing individual-specific factors included Age (to account for age-related trends in growth) and Sex (to account for possible sex-specific variability in growth rates), as well as an interaction (Age × Sex).

Furthermore, Year (calendar year of increment formation) and
Cohort (year-class) were included to investigate inter-annual and between year-class growth differences.
Extrinsic variables were then included to investigate the influence of the environment on cod growth. Since temperature is known as one of the most significant factors influencing fish growth due to its effects on metabolism and prey availability (Brander, 1995;Michalsen, Ottersen, & Nakken, 1998) Bochkov, 1982;Tereshchenko, 1996). Additionally, two large-scale hydro-climatological indices were included: the winter North Atlantic Oscillation (NAO) calculated as the standardized NAO from December to March (Hurrell, 1995), and the Atlantic Multidecadal Oscillation (AMO; Enfield, Mestas-Nuñez, & Trimble, 2001;Kerr, 2000). Both indices are derived from temperature anomalies but usually reflect larger changes in the hydrology, the circulation and the distribution of the water masses of the North Atlantic, which have been shown to cause large-scale fluctuations in many fish populations including cod, affecting for example spatial distribution, prey availability or recruitment (Alheit et al., 2014;Drinkwater et al., 2003;Stige, Ottersen, Brander, Chan, & Stenseth, 2006). Given the temporal range covered by our chronology, the AMO index could additionally provide insights on the response of cod growth to lower frequency changes (20-40 years) in the environmental conditions.

| Statistical analysis
Hierarchical linear mixed-effects modeling (Zuur, Ieno, Walker, Saveliev, & Smith, 2009) was used to model increment width (Increment) in relation to the selected factors. This approach helps partition different sources of variability while taking into consideration repeated measurements within individuals (Morrongiello & Thresher, 2015;Weisberg et al., 2010). Prior to analysis, increment width (referred as growth in the manuscript) and Age were log transformed to linearize the relationship and meet model assumptions.
All intrinsic and extrinsic factors were mean-centered to aid model convergence and interpretation (Morrongiello & Thresher, 2015;Morrongiello et al., 2014). Models with random intercept and age slope for each individual (FishID) were tested to allow individualspecific growth trajectories. These models were then extended to incorporate the effects of sex (intrinsic fixed effect), calendar year of formation and cohort (random effects), and the effects of both climatic and population variables (extrinsic fixed effects). The most complex model can be represented by: where y ijkl is the increment width relative to age at formation (x ijkl ) for the ith fish at age j from year k and cohort l; α 0 and β 0 are the fixed mean intercept and slope describing population-wide decline in increment width in relation to age; F i and F are the fish-specific random intercept and slope describing individual growth trajectories in relation to age; Y k and Y are the year-specific random intercept and slope in relation to age; C l and C are the cohort-specific random intercept and slope in relation to age; and f(·) represents additional fixed factors (e.g., intrinsic factor sex, and extrinsic factors described in Section 2.3).

| Base model and intrinsic sources of growth variation
Initially, different random effects were tested in combination to optimize the random structure of the base model (Table S2).
Random intercepts for individual fish (FishID), year of increment formation (Year) and fish year-class (Cohort) were tested to control for correlations among growth increments within individual fish, year or cohort. Random by-FishID, by-Year and by-Cohort slopes for the Age effect were also examined to test for differences in age-dependent growth trajectories among individuals, years and  (2002).
year-classes. Models were fitted to the entire growth time series from 1924 to 2014 (n = 28,504 increments, 3,894 FishID) with the full intrinsic fixed-effect structure (Age × Sex) and restricted maximum likelihood estimation (REML; Morrongiello & Thresher, 2015;Zuur et al., 2009). Models were ranked using Akaike's information criterion corrected for small sample sizes (AICc) and the optimal model was selected (Burnham & Anderson, 2004). If the difference in AICc (ΔAICc) between the highest and a second highest ranked model was <2, the two were considered to be equally supported (Burnham & Anderson, 2004).
To explore the optimal intrinsic fixed-effect structure, combinations of Age, Sex and Age × Sex interaction were fitted with the maximum likelihood (ML) method with the previously selected random structure and compared using AICc (Table S3). The optimal model was then refitted with REML to obtain unbiased estimates (Zuur et al., 2009). For each model, the marginal (variance explained by intrinsic effects alone) and conditional (variance explained by both intrinsic and extrinsic effects) R 2 metrics were calculated (Nakagawa & Schielzeth, 2013). Growth synchrony between individuals from a given year or a given year-class was estimated by calculating the intra-class correlation coefficients (ICC) for optimal Year-only and Cohort-only models. These coefficients estimated the correlation of fish growth within years and year-classes by measuring the relative similarity as proportions (Koch, 2006).
The master growth chronology was developed by extracting the best linear unbiased predictor (BLUP) for the Year random effect from an optimal intrinsic-only model to visualize mean inter-annual variations in growth at the population level (Morrongiello & Thresher, 2015). This Year random effect essentially represents the yearly variation in growth which can be associated with the extrinsic environmental factors, after accounting for age and sex effects. BLUP for the Year random effect was compared with mean size-at-age and growth rates derived from survey data from 1981 to (ICES, 2018. A sequential t test analysis (STARS) was applied to the extracted chronology to identify significant shifts in mean population growth. This method provides an estimation of statistically significant changes in the mean level and magnitude of fluctuation in a time series and their probability (Lindegren, Diekmann, & Möllmann, 2010;Rodionov, 2004Rodionov, , 2006. The cut-off length (i.e., the minimum length of a shift for which its magnitude remains intact) was set to 5 years with a significance level to .05 so that shifts in cod growth were identified only when the shift lasted at least 5 years.
"Prewhitening," a process that eliminates or reduces short-term autocorrelation to enable the detection of trend changes, was implemented prior to the sequential t test analysis to remove the red noise component of the growth time series (Rodionov, 2006;Smoliński & Mirny, 2017).

| Influence of extrinsic factors on cod growth
External factors influencing fish growth were assessed through a series of model comparisons by adding potential extrinsic factors to the best intrinsic structure previously selected. First, the influence of environmental conditions on annual growth was investigated by fitting combinations of sea temperature (KolaT), NAO and AMO index to the entire growth time series from 1924 to 2014 (n = 28,504 increments, 3,894 FishID; Table S4a). Second, the effects of cod biomass (Stock), capelin biomass (Cap_ts) and fishing mortality (Fbar5.10) were evaluated by comparisons of optimal growth models refitted to a more limited time period, from 1973 to 2014, due to the temporal limitations of the capelin and fishing mortality time series (n = 12,171 increments, 1,837 FishID; Table S4b). Finally, global extrinsic models incorporating combinations of both environmental, population and fishing factors were refitted to data from 1973 to 2014 (n = 12,171 increments, 1,837 FishID; Table S4c) and compared to investigate the influence of all extrinsic factors on cod growth. Potential collinearity between explanatory variables in the final models was investigated by calculating their variance inflation factors (VIF), with a more conservative threshold defined at VIF < 2 (Borcard, Gillet, & Legendre, 2018). Assumptions of normality and homoscedasticity were evaluated visually by inspecting model residuals. Linear trends in mean growth over the years and agedependent effects of extrinsic variables on growth were also explored but were found to be negligible and are therefore not presented.

| Mechanistic modeling of cod growth and extrinsic factors
To Confirmatory path analysis based on the piecewise fitting of component hierarchical models was conducted to model hypothesized causal relationships between climate, population dynamics and fish growth (Lefcheck, 2016;Shipley, 2009 where indirect effects were described as one variable affecting another through a simultaneous response of a third. All analyses were conducted using the R scientific computing language (R Core Team, 2019) with the packages tidyverse (Wickham, 2017), AICcmodavg (Mazerolle, 2019) and MuMIn (Bartoń, 2019). Linear mixed-effect modeling was performed using the lme4 (Bates, Mächler, Bolker, & Walker, 2015) and effects (Fox, 2003;Fox & Weisberg, 2018 packages. SEM was developed using the package piecewiseSEM (Lefcheck, 2016) and the result flowcharts were constructed using DiagrammeR (Iannone, 2019).

| Temporal variations in growth
The 28 (Figure 3b). As commonly found in fish otoliths, increment width declined significantly with age, indicating that growth was age-dependent ( Figure 4a).
The growth time series extracted from the BLUP for the Year random effect showed significant interannual variation of NEA cod population growth (Figure 5a). Periods of peak growth were seen in 1937, 1969, 1985, 1990 and 2000, while periods of lowest growth were observed in 1942, 1948, 1987, 1993

| Intrinsic sources of growth variation
Age had the greatest influence on growth among the tested variables, with predicted increment widths strongly declining as fish  (Table 2). A significant Age × Sex interaction showed that males older than 6 years old grew slightly slower than females, which corresponds to the onset of maturity (Figure 4b). Calculated AICc values for the different model random structures indicated the strongest support for the incorporation of a random intercept and Age slope for FishID and Year, and a random intercept for Cohort (Table S2).
The optimal intrinsic fixed-effect structure consisted of Age in interaction with Sex (Table S3). The associated conditional R 2 (.66) indicated that intrinsic and random factors together explained around 66% of the observed variance in cod growth.

| Influence of extrinsic factors on cod growth
The inclusion of climatic factors (temperature, NAO and AMO indices) to the optimal base model fitted to data from 1924 to 2014 revealed significant effects on growth. Model predictions revealed that growth was positively correlated with temperature and weakly correlated with NAO, while AMO had a significant negative effect on cod growth (Table S5). Comparisons of AICc showed the strongest support for models incorporating temperature, and the optimal identified model incorporated both temperature and AMO index (Table S4a). Because there was similar strong support for the model including all three factors (ΔAICc < 2), it was used to estimate the predicted effects on growth. Second, the addition of population factors (cod and capelin biomass) and fishing mortality to the optimal base model refitted to data from 1973 to 2014 significantly improved model fit and revealed significant relationships to growth (Table S4b). In models restricted to population factors for the period 1973-2014, both cod biomass and fishing mortality had a significant negative influence on growth, while capelin biomass had a positive effect (Table S6).  (Table 3). Variability of the AMO was also associated with a significant decrease in growth (−7.20%). Increased growth was associated with variations of temperature (8.48%) and capelin biomass (4.67%). Comparison of AICc showed the best fit for the model comprising all factors except fishing mortality and NAO (Table S4c). However, the inclusion of fishing mortality had a similar statistical support and the variable was significant in the models previously refitted without climatic factors (Table S4b). Fishing mortality was consequently retained in the final model, where we found a significant but weak negative influence on growth (−2.78%; Table 2; Figure 6a-c). The associated marginal R 2 (.48), the proportion of variance explained by fixed effects alone, indicated that the inclusion of environmental variables explained more of the growth variance than the intrinsic-only model (.45) refitted to the same period (Table S3).

| D ISCUSS I ON
We used an otolith-based approach to reconstruct a 91-year-long continuous and annually resolved growth history of the NEA cod population. After accounting for intrinsic effects such as the agedependent decline in growth, we identified density dependence associated with cod population size as the strongest extrinsic factor influencing cod growth throughout the time period. Our analysis also highlighted significant variability in cod growth which was associated with fluctuations in prey availability, as well as with environmental changes at both the local (i.e., temperature) and large (i.e., AMO) scales. Potential negative effects of fishing mortality on cod growth were also identified, although the relationship was weak and hard to disentangle from density dependence. full growth history for every individual and therefore account for individual-specific variability in growth rates, as opposed to size-atage data which only provides information on the sizes at capture.
Fish somatic and otolith growth is generally age-dependent (Campana, 1990;Lee, 1920) and this has been shown for different species in similar otolith-derived biochronologies (Doubleday et al., 2015;Izzo et al., 2016;Smoliński & Mirny, 2017). In this study, age was also, as expected, the strongest factor influencing growth of NEA cod, and growth decreased significantly as fish got older. We also highlighted small age-specific sexual differences in growth rates where males grew slightly slower than females as they got older.
Because sexual differences in cod growth were small, this probably reflected the decrease in growth following sexual maturation since males tend to mature earlier and at a smaller size (Brander, 1994). Growth synchrony was higher across years, irrespective of age, than across cohorts (year-classes) but remained relatively low in both cases, which is a consistent pattern seen in other otolith chronologies based on mixed-effect modeling (Martino, Fowler, Doubleday, Grammer, & Gillanders, 2019;Morrongiello & Thresher, 2015;Smoliński, 2019). The lower level of inter-individual growth synchronicity in fish and more specifically NEA cod, as opposed to relatively static organisms such as bivalves or corals, may be caused by largescale migrations undertaken once they reach sexual maturity (Höffle et al., 2014).
Our growth chronology highlighted significant alterations of NEA cod growth that can be correlated to important fluctuations in fish populations, and in particular that abrupt shifts in cod growth since 1973 seemed to closely reflect variations in cod and capelin biomasses. Density-dependent decrease in growth as a result of enhanced competition for food and resources has been identified in many species and is known to be one of the key factors regulating fish populations (Beverton & Holt, 1957;Lorenzen & Enberg, 2002).
A recent study showed that the importance of density-dependent processes in NEA cod had changed throughout the last century, and in particular that models based on contemporary data  showed better support for strong density-dependent growth (Eikeset et al., 2016). Despite a period of heavy exploitation and decreasing stock size between the 1950s and the 1980s, NEA cod biomass in the Barents Sea has since been rising steadily and in 2014 reached the highest biomass ever recorded . Our model predictions also suggest strong density-dependent growth, as cod biomass was identified as the strongest extrinsic factor and was significantly correlated with a large decrease in cod growth (more than 20% between the lowest and highest recorded biomass). Since our final growth model was refitted to data from 1973 to 2014, it therefore likely captured a period of strong density dependence correlated to the significant increase in stock biomass.
Temporal variations in the strength of density-dependent growth have been observed when ecosystem conditions change significantly, and especially when resources become scarcer (Casini, Rouyer, Bartolino, Larson, & Grygiel, 2014;Rueda, Massutí, Alvarez-Berastegui, & Hidalgo, 2015). The significant positive relationship found between capelin biomass and cod growth within our models is in line with recent studies on the importance of capelin biomass and cod growth within our models is in line with recent studies on the importance of capelin in the diet of NEA cod (Bogstad, Gjøsaeter, Haug, & Lindstrøm, 2015). After the collapse of the NSS herring population at the end of the 1960s and the subsequent increase in capelin biomass (Lees, Pitois, Scott, Frid, & Mackinson, 2006), capelin has become an essential item in the diet of cod (Holt et al., 2019).
The three capelin collapses that occurred in 1985, 1993(Gjøsaeter, Bogstad, & Tjelmeland, 2009) may have negatively affected cod feeding, which is likely reflected in our chronology by sudden drops in cod growth during these periods. In addition, it has been hypothesized that particularly strong herring year-classes might be a key factor in capelin recruitment failure due to their predation on capelin larvae (Hallfredsson & Pedersen, 2009;Hjermann et al., 2010). Our study did not include herring biomass as a direct factor influencing cod growth but it was tested as a factor of capelin biomass within the SEM, where we found a strong and significant negative effect. It appears likely that large herring year-classes were partly responsible for capelin collapses, which, in turn, led to a lower prey availability and a significant decrease in cod growth.
Furthermore, a collapse of the capelin population and a low herring biomass occurred simultaneously in the mid-1980s and resulted in a temporary shift in cod diet toward less-preferred prey such as Euphausiids and Hyperiids (Reid, Battle, Batten, & Brander, 2000).
The abrupt decrease in cod growth we observed between 1985 and 1988 likely reflects the reduction in growth associated with this temporary sub-optimal diet, although the duration was too short to be identified as a significant shift in growth within the STARS analysis. In addition, cannibalism is known to occur in many cod populations and is especially prevalent in NEA cod (Bogstad, Lilly, Mehl, Palsson, & Stefánsson, 1994;Yaragina, Bogstad, & Kovalev, 2009).
The observed variability in cod growth after 1980 could thus reflect changes in prey availability and cannibalism consistent with large fluctuations in capelin abundance (Eriksen, Skjoldal, Gjøsaeter, & Primicerio, 2017;Gjøsaeter et al., 2009;Holt et al., 2019). Density dependence would therefore be important even during the period of relatively lower cod biomass due to increased cannibalism and lack of alternative prey.
Temperature was the second strongest extrinsic factor identified in our study. Cod growth increased significantly with sea temperature in all of our models, which is in agreement with our initial hypothesis that it is positively influenced by warming conditions. Sea temperature is known to significantly influence fish growth and may correlate with growth rates through changes in direct (i.e., metabolism, individual fitness, duration of the growth season) or indirect (i.e., food availability) processes (Brander, 1995;Neuheimer & Grønkjaer, 2012). Temperature signals in growth have already been identified in multiple otolith-based biochronologies, although the strength and direction vary between species and locations (Gillanders, Black, Meekan, & Morrison, 2012;Izzo et al., 2016;Martino et al., 2019).
We are confident in the relationship found here since our analysis is based on the Kola section temperature time series, which has been shown to accurately represent temperature variability in the southwestern regions of the Barents Sea where Atlantic water masses are prevalent (Dippner & Ottersen, 2001). We also used in situ measure- The positive relationship with temperature is also evident in the significantly higher growth observed in the earliest parts of the chronology. During the 1920s and 1930s, the North Atlantic underwent a significant warming that had important environmental and ecosystem effects (Drinkwater, 2006;Johannessen et al., 2004).
The Barents Sea in particular showed a drastic increase in primary and secondary production, which likely created favorable conditions for fish growth (Drinkwater, 2006). This also coincided with high recruitment, high population size, and a northward expansion of the feeding and spawning grounds observed for NEA cod during the same period (Hylen, 2002;Sundby & Nakken, 2008). However, while our analysis indicates a positive correlation between growth and warming conditions, it must be noted that cod growth has decreased despite record warming of the Barents Sea in the most recent decade (Levitus, Matishov, Seidov, & Smolyar, 2009;Lind et al., 2018).
The ongoing warming trend has been associated with significant changes in the ecosystems and the hydrography of the region (see e.g., Fossheim et al., 2015;Lind et al., 2018;Smedsrud et al., 2013), which could eventually affect the productivity of cod under the continued warming projected for the next decades (Stocker et al., 2014).
The fast alteration of the Barents Sea conditions combined with the highest cod biomass recorded in more than 50 years could consequently explain the sharp decline in growth observed since 2008, although edge effects were likely in effect as these years only comprised fish from the same cohort and thus were represented by increasingly older individuals.
In addition to local changes in the environmental conditions, we identified significant shifts in cod growth which could be linked to large-scale changes in the Barents Sea, especially between 1993 and 1999. Significant shifts in fish population biomass and life history in the Northeast Atlantic have been historically correlated to largescale climatic variations (Alheit et al., 2014;Hjermann, Stenseth, & Ottersen, 2004;Lehodey et al., 2006). The concept of regime shift has, in turn, grown out of a number of studies showing significant ecosystem-level changes during the late 1980s and the mid-1990s in the Northern Hemisphere (Auber, Travers-Trolet, Villanueva, & Ernande, 2015;Beaugrand et al., 2015;Reid, Borges, & Svendsen, 2001). The mid-1990s ecosystem shift in particular has been associated with a complex series of atmosphere-ocean changes characterized by a significant weakening of the NAO and a switch of the AMO from a cold to a warm phase (Alheit et al., 2019;Hughes, Holliday, & Gaillard, 2012;Robson, Sutton, Lohmann, Smith, & Palmer, 2012).
Yet, despite periods of substantial alterations of both the climatic conditions and the abundances of several species during the last decades, no persistent ecological regime shift has been identified for the Barents Sea (Johannesen et al., 2012). Our model reveals that AMO had a strong negative influence on cod growth and that the abrupt decrease in growth rates between 1993 and 1999 correlates with the onset of its warmest phase. However, while AMO is defined as a detrended indicator of sea temperature anomalies (Kerr, 2000), its influence on the marine ecosystems seems to be primarily linked to large-scale changes in the strength and direction of the water masses circulating in the North Atlantic (Alheit et al., 2019). In particular, the dynamics of many small pelagic fishes vary significantly with the AMO phases, not directly due to the temperature anomalies but through complex changes in the coupled atmosphere-ocean system (Alheit et al., 2014). The significant negative relationship we found between AMO and cod growth, which contrasts with the positive influence of warming conditions discussed earlier, may thus be more representative of its indirect effect through large fluctuations in small pelagic fish assemblages and prey availability for NEA cod.
Finally, we found a significant but weak negative relationship between growth and fishing mortality. Fisheries-induced evolution is well studied and has been associated with changes in growth and maturity in many populations (Kuparinen & Merilä, 2007;Swain, Sinclair, & Mark Hanson, 2007), which makes growth signals an effective proxy to identify its presence in a given population. While the present relationship may indicate that multi-decadal fishing have caused a decline in cod growth through selective harvesting of faster-growing individuals (Enberg et al., 2012;Heino & Godø, 2002;Sinclair et al., 2002), the relationship was very limited and models that did not include fishing mortality had a similar statistical support.
In contrast, it is possible that fishing mortality had an indirect positive influence on growth through its significant effect on cod biomass and consequently on density-dependent effects, as suggested by the optimized SEM. In the present study, it is therefore hard to disentangle the effects of fishing mortality and density dependence, and fisheries-induced effects on cod growth are ultimately difficult to interpret.
In addition to the interannual variability, our results showed cod growth variations at the year-class level, which indicate a persistent growth signal in fish born during the same spawning season.
Cohort-specific growth variability is an important aspect of temporal variations in growth, as it likely reflects intrinsic differences in the systematic response to environmental conditions (Morrongiello & Thresher, 2015). This is often associated with carry-over effects (Murphy, Jenkins, Hamer, & Swearer, 2013) wherein the growth trajectory might be influenced by the early conditions, for example, individuals born in a particular year-class where environmental conditions were detrimental (or beneficial) for future growth.
Our results showed that year-specific and cohort-specific growth extracted from the BLUP for each random effect were often negatively correlated, and cohorts born in years with poorer than average growth often displayed higher than average growth trajectories throughout their entire life. The two signals are not directly comparable: one represents the overall growth deviation from the population average pooled across all age-classes for a single year, the other the deviation from the population average in the entire growth trajectories of individuals born in a given year. However, this negative relationship could be the marker of some form of compensatory growth response, indicating that fish born in poor conditions grew faster following intervals of slow growth. Such compensatory growth is widespread in juvenile fish especially in an aquaculture context (e.g., Urbinati, Sarmiento, & Takahashi, 2014), and this response was shown to sometimes overshoot the "average" growth trajectory (Hayward, Noltie, & Wang, 1997). This hypothesis seems plausible for a near-arctic population of Atlantic cod, as Schultz, Lankford, and Conover (2002) showed that high-latitude populations in highly seasonal environments were prone to strong compensatory growth responses. It is however difficult to quantify this effect with precision without experimental data, especially as the growth synchrony within cohorts was ultimately minimal. Year-class strength can also have persistent effects on growth through competition for resources among same-age individuals (Amundsen, Knudsen, & Klemetsen, 2007). These results suggest that years with unfavorable conditions may have been correlated with weaker year-classes and therefore relaxed early-life competition, which is in line with the importance of density-dependent effects highlighted in our study.
The combination of mixed-effects models and SEM has only been used recently in fish studies (e.g., Taylor In addition, many model-derived variables such as fishing pressure or stock biomass carry a non-negligible level of uncertainty as we progress back in the time series beyond the period covered by survey data, due to increased model assumptions. The biological basis for the assumed proportionality between somatic and otolith growth is well accepted and has been experimentally verified (Hüssy & Mosegaard, 2004;Li et al., 2008;Yamamoto, Ueda, & Higashi, 1998). However, as shown in this study, this relationship can be more subtle due to inter-individual differences in otolith growth and mineralization rates that stem from different individual environmental and feeding conditions (Grønkjaer, 2016). One advantage of the mixed modeling approach over traditional biochronology studies is that we partly controlled for this by including random FishID effects to account for this type of inter-individual variation (Morrongiello & Thresher, 2015;Weisberg et al., 2010). In the future, it would nonetheless be interesting to conduct experimental studies to quantify how these effects may influence sclerochronological reconstructions. Finally, the sampling strategy we adopted due to the uneven availability of older fish throughout the whole period means that our study could be biased toward earlier maturing, potentially faster growing individuals during the earliest parts of the growth chronology. This could have affected some of our results and more specifically the growth variation found during the earliest years. We nonetheless remain confident in our interpretations, as our hypotheses were based on relevant ecological knowledge of the NEA cod dynamics in the Barents Sea and our core results were focused on the period 1973-2014, where the available data and therefore our modeling approach were much more robust.

| CON CLUS ION
Altogether, we provide evidence that NEA cod growth has responded to the joint effects of climate change, population dynamics and harvesting throughout the 20th and beginning of the 21st centuries. We found that density dependence had the strongest influence on cod growth, but that climatic factors also had impor-

ACK N OWLED G EM ENTS
The authors thank Erlend Langhelle (IMR) for his help with retrieving and processing the otoliths, as well as Anders Thorsen (IMR) and Norbert Vischer (University of Amsterdam) for their collaboration on the plugin used for imaging and annotating the samples. We also thank Hildegunn Mjanger (IMR) for her assistance with age reading and validation of the otolith sections.

CO N FLI C T O F I NTE R E S T
The authors declare they have no conflict of interest.

CO M PLI A N CE WITH E TH I C A L S TA N DA R DS
The authors declare they have complied with ethical standards.