Potential sources of bias in the climate sensitivities of fish otolith biochronologies

Analysis of growth increments in the hard parts of animals (e.g., fish otoliths) can be used to assess how organisms respond to variability in environmental conditions. In this study, mixed-effects models were applied to otolith data simulated for two hypothetical fish populations with assumed biological parameters and known growth response to environmental variability. Our objective was to assess the sensitivity of environment–growth relationships derived from otolith biochronologies when challenged with a range of realistic ageing errors and sampling regimes. We found that the development of a robust biochronology and the precision of environmental effect estimates can be seriously hampered by insufficient sample size. Moreover, the introduction of even moderate ageing error into the data can cause substantial underestimation of environmental sources of growth variation. This underestimation diminished our capacity to correctly quantify the known environment–growth relationship and more generally will lead to overly conservative conclusions concerning the growth response to environmental change. Careful study design, reduction of ageing errors, and large sample sizes are critical prerequisites if robust inferences are to be made from biochronological data. Résumé : L’analyse des incréments de croissance préservés dans les parties dures d’animaux (comme les otolites de poissons) peut être utilisée pour évaluer les réactions des organismes aux variations des conditions ambiantes. Des modèles à effets mixtes ont été appliqués à des données d’otolites simulées pour deux populations de poissons hypothétiques en intégrant des paramètres biologiques présumés et une réaction connue de la croissance à la variabilité des conditions ambiantes. L’objectif consistait à évaluer la sensibilité des relations entre le milieu et la croissance obtenues de biochronologies reposant sur les otolites étant donnée une fourchette d’erreurs réalistes de détermination de l’âge et de régimes d’échantillonnage. Nous constatons qu’un échantillon trop petit peut sérieusement entraver l’élaboration d’une biochronologie robuste et limiter la précision des estimations des effets du milieu ambiant. L’introduction dans les données d’erreurs ne serait-ce que modérées de détermination de l’âge peut en outre entraîner une sous-estimation considérable des sources ambiantes de variation de la croissance. Cette sous-estimation réduit notre capacité de quantifier correctement la relation milieu-croissance connue et, plus généralement, se traduira par des conclusions trop modérées concernant la réaction de la croissance à des changements des conditions ambiantes. Une méthodologie judicieuse, la réduction des erreurs de détermination de l’âge et des échantillons de grande taille sont des conditions préalables essentielles à l’établissement d’inférences robustes à partir de données biochronologiques. [Traduit par la Rédaction]


Introduction
Otoliths are calcified structures located in the inner ear of fish that commonly show growth bands corresponding to daily, seasonal, or annual growth patterns (Campana and Thorrold 2001). For many fish species, otolith size is highly correlated with fish size (e.g., Harvey et al. 2000); therefore, the width of an increment can be used as a proxy of somatic growth and reflect individual responses to environmental variability (Morrongiello et al. 2011). The somatic growth changes estimated with such a biochronological approach integrate the direct and indirect impacts of environmental conditions on individual fish over time (Rountrey et al. 2014). Somatic growth rate is an important driver of population dynamics (Smoliń ski 2019), as it affects an individual's size and age at maturity (Heino et al. 2002;Uusi-Heikkilä et al. 2015), its reproductive output (Barneche et al. 2018), and size-dependent processes such as recruitment, predation, resource acquisition, and competition (Bergenius et al. 2002;Post 2003). Understanding the drivers of growth variability can therefore inform the sustainable management of harvested stocks (Whitten et al. 2013;Stawitz and Essington 2019).
More than 800 000 otoliths are collected worldwide annually (Campana and Thorrold 2001), many of which are subsequently stored in the archives of fishery institutes. These archives thus present an enormous source of material for sclerochronological studies on fish growth across a wide range of spatial and temporal scales (Morrongiello et al. 2012). Traditionally, most studies in fishery science focused on the back-calculation of body size made possible by measurements of increment widths in fish otoliths (Campana 1990;Francis 1990). More recently, advances in analytical techniques (Black et al. 2005;Morrongiello and Thresher 2015) have yielded an increasing number of studies focusing on the environmental drivers of year-to-year variability in fish growth. These studies generally aim to extract population-wide estimates of temporal growth variability from the available otolith collection and subsequently relate this to variability in the environment (e.g., Black et al. 2013;Morrongiello et al. 2019;Smoliń ski and Mirny 2017;van der Sleen et al. 2018).
Generating accurate predictions of how fish respond to different scenarios of environmental change (e.g., climate change and fishing pressures) can enhance our ability to manage marine resources (Hoegh-Guldberg and Bruno 2010). Despite the wide availability of biological materials, the predictive potential of otolith biochronologies remains largely untapped (Morrongiello et al. 2012; Barrow et al. 2018). Since most of the studies in the field of otolith biochronology are based on archived materials, decisions regarding sample size and individual otolith selection affect the general robustness of the analysis. Moreover, it has been suggested that some commonly used sampling strategies (e.g., selection of oldest individuals) may bias ecological inferences that are based on the otolith increment analyses (Morrongiello et al. 2012). Ageing accuracy and precision are also a challenge in otolithbased studies and can affect ecological inferences (Campana 2001). Errors caused by falsely added or missed increments can propagate back through time and eliminate high-frequency variability in a time series. This reduced variability can prevent the identification of environmental effects on organisms' growth , which ultimately can lead to incorrect conclusions concerning a species' vulnerability to environmental change (Rowland et al. 2011). Simulation studies offer great potential to evaluate the effects of sample size, sampling strategy, and underlying measurement errors, as well as their interactions, on the statistical power of biochronological inferences.
Here, we applied a mixed-effects modeling framework (Weisberg et al. 2010;Morrongiello and Thresher 2015) to analyze simulated otolith-based growth data from two hypothetical species with assumed underlying biological parameters (i.e., overall mean annual growth, coefficients of age and environmental effects, and intra-individual covariance structures). We assessed the sensitivity of these models to a series of scenarios where we imposed ageing error and alternate sampling regimes. We hypothesized that our ability to robustly detect the true growth response to environmental variability would be most sensitive to ageing errors, but that large sample size could partly ameliorate the impacts of these errors, at least at lower-frequency time scales. We conclude with recommendations for the analysis of otolith growth chronologies in future studies.

Simulation of otolith growth data
A simplified mixed-effects model with a random slope and intercept for age was assumed for each fish and used to model growth from otolith data: (1) where y ij is the annual otolith growth y for fish i at age j; ␣ 0 is the overall mean annual growth intercept; ␣ i F is a random intrinsic effect for fish i; ␤ A is a coefficient for the fixed effect age; b Ai F is a random slope for age for fish i, correlated with ␣ i F ; ␤ E is a coefficient for the fixed environmental effect, ⌺ is the covariance matrix between random intercept and random slope, and is assumed to be independent errors with mean zero and common variance 2 (Morrongiello and Thresher 2015).
We used the mean annual sea surface temperature for southwestern Icelandic waters (60°N-65°N, 27°W-17°W) from the Had-ISST model (Rayner et al. 2003) as the environmental variable influencing fish growth. This allowed us to explore growth responses to a time series with typical climatic variance and autocorrelation structures (Fig. 1a). The age distribution for the species considered relatively short-lived (Fig. 1b) reflected catch numbersat-age of Icelandic cod (Gadus morhua) for the period 1955(ICES 2018. Age of the second hypothetical species was assumed to be uniformly distributed between 91 and 100 years (Fig. 1e). The latter was intended to simulate sampling focused on the collection of old individuals from a population of a relatively long-lived species (e.g., from the family Trachichthyidae; Thresher et al. 2007Thresher et al. , 2014. We assumed stationary (remaining the same across years) age distributions for both species.
Parameter values for the fixed and random effects in the simulation were obtained from the mixed-effects model (eq. 1) fitted to real otolith growth data from the Icelandic cod stock (S. Smoliński, unpublished data; Table 1). The same parameters were used for both species to make simulation results more comparable. A negative effect of age on annual growth was simulated, which is typical for most fish species and often dominates the otolith growth patterns (Figs. 1c,1f). The magnitude of the environmental effect was set to be ϳ13% of the age effect (Figs. 1d,1g), consistent with the value observed in an informal review of scientific studies in which mixed-effects modeling of otolith growth data was applied (see online Supplementary materials for the list of publications 1 ). The simulated data set included 500 individuals caught per year for the period 1901-2000 for the short-lived species (Fig. 1h) and 200 individuals caught per year for the period 1981-2000 for the long-lived species (Fig. 1i). The data for the long-lived species included more observations per individual (annual growth increments), and therefore the number of individuals was limited due to computational costs.

Simulation of otolith ageing errors and different sampling designs
Five scenarios of ageing error were introduced to each species' growth data to reflect different degrees of accuracy and precision in the ageing methods. The age of all fish used in the "perfectly aged" scenario was correct, and therefore all increments were correctly dated. Age error matrices for the remaining four scenarios were simulated following McBride (2015). In the "no bias" scenario, random error (without systematic bias) was introduced by adding to the true age a normally distributed error with mean = 0 and standard deviation = 0.1 of the true age; values were rounded to the integer. This resulted in an average coefficient of variation (CV) between the interpreted and the true age of 4.26% for the short-lived and 5.71% for the long-lived species. The "bias (-1)" scenario reflected consistent underestimation of age by an average of 1 year. In this case, a normally distributed error with mean = -1 and standard deviation = 0.1 of the true age was added (CV = 24.58% for short-lived and CV = 5.86% for long-lived species). The estimated age in the "bias (-10%)" scenario was 10% less than the true age, and a normally distributed error was added with mean = 0 and standard deviation = 0.1 of the true age (CV = 7.72% for short-lived and CV = 9.47% for long-lived species). Two types of biases were combined for the "bias (+1, -10%)" scenario by adding normally distributed error (mean = 1 and standard deviation = 0.1 of the true age) and subtracting 10% of the true age (CV = 11.19% for short-lived and CV = 8.86% for long-lived species; Fig. 2). For the incorrectly aged individuals under different ageing error scenarios, the overestimation of age in the individual time series of Fig. 1. Visualization of the data used in the study: (a) time series of mean annual sea surface temperature for the area 60°N-65°N, 27°W-17°W from the HadISST model, (b, e) assumed static fish population age structures, where AAC is age at capture, (c, f) simulated annual increment measurements, (d, g) interaction plot of simulated mean annual growth in relation to age and the environmental variable, (h, i) temporal range of simulated annual increment time series for short-lived (b, c, d, h) and long-lived (e, f, g, i) species. otolith growth was simulated by the addition of false rings that split a randomly chosen annual increment into two. Underestimation of age was simulated by combining two randomly chosen neighboring annual increments into one. Therefore, all increments prior to the introduced error were incorrectly dated.
A series of sampling regimes were applied to the simulated data sets. First, we subsampled 1, 5, 10, 20, or 50 individuals per collection year. The ageing error introduction and subsampling were repeated 1000 times for each combination of our five ageing error schemes and five sampling regimes. In total, this simulation resulted in 25 000 data sets per species.
Many otolith collections are composed of sporadically collected samples. Hence, it is appealing to select older individuals for further analysis, as their longevity helps to maximize the length of the biochronology and enlarge the number of observations (increment measurements). In the absence of cross-dating, such sampling may, however, be associated with a greater age reading error (see Fig. 2a). We explored the impact of this possibility only for the short-lived species by sampling 10 individuals from different age groups, ranging from 3 to 10 years (Fig. 3). This simulation was repeated 1000 times for five ageing error schemes and all possible ranges of age groups (36 combinations), resulting in 180 000 data sets in total. We omitted similar tests for the long-lived species due to the expected low differences in the ageing errors within the simulated age ranges (91-100 years) and the high computational and data storage cost.

Evaluation of model performance
Two separate linear mixed-effects models were fitted to each data set. The first model (M1) followed eq. 1. Distributions of 1000 parameter estimates from independent simulation runs of M1 were visualized, and the variance of these estimates and root mean squared error (RMSE, with "true" parameter as reference) were calculated to assess accuracy and precision.
The fixed effect of the environmental factor was excluded from the second model (M2) and replaced by an additional year random effect. The following formula was used for M2: (2) where y ijk is annual growth y for fish i at age j from year k, and ␣ k Y is a random extrinsic environmental effect for year k (intercept).
Other terms are consistent with those presented in eq. 1. When the fixed environmental effect coefficient is excluded from the model (as in M2), random effects of year reflect interannual changes in the mean fish growth associated with pooled extrinsic sources of variation, after accounting for intrinsic fac-tors (Weisberg et al. 2010;Morrongiello and Thresher 2015). These estimates of random effects of year are treated as a biochronology, which is often of primary interest to ecologists. Best linear unbiased predictors (BLUPs) for effects of year were extracted from M2 to obtain temporally resolved model estimates of pooled environmental effects on fish growth. This allowed us to investigate covariance between biochronological and environmental signals with cross-wavelet analysis. Wavelet transforms are often used for analyzing a time series characterized by different power at different frequencies. Cross-wavelet coherence analyses are used for detecting covarying patterns between two time series both across time and frequencies. The cross-wavelet plots depict coherence, which is a measure of the intensity of covariance, and the phase between two time series. The statistical significance of coherence shown on plots is estimated against a red noise background. We used cross-wavelet coherence analysis with the Morlet wavelet function and 1000 Monte Carlo simulations to obtain the significance level (see Grinsted et al. 2004 for details on wavelet coherence analysis).
All combinations of the error matrices and sampling designs were tested with 1000 random permutations of the independent ageing error introduction, sampling, and model fitting (Fig. 2). To make results comparable, we used the same sets of randomly chosen individuals in parallel runs under different ageing errors. All analyses were conducted with the R scientific computing language (R Core Team 2018). Data were simulated and linear mixedeffects models were fitted with lme4 using restricted maximum likelihood (Bates et al. 2019). Wavelet coherence analysis was performed with biwavelet package (Gouhier et al. 2018).

Precision and accuracy of environmental parameter estimates
Our results show a gradual increase in the precision (decrease in the variance) of environmental parameter estimates, with increasing sample size (n per year) for both short-lived and longlived species (Figs. 4 and 5). Increasing precision was also observed with reduction of ageing error in the simulated data sets (Figs. 4 and 5). The short-lived species had lower parameter estimate precision for a given age scenario and sample size than the long-lived species (Figs. 4 and 5). Environmental effects were underestimated by between 56.8% and 58.3% in inaccurately aged long-lived data sets ("no bias", "bias (-1)", "bias (-10%)", "bias (+1, -10%)"; Fig. S1 1 ; Table S1 1 ). In the "bias (+1, -10%)" ageing error scenario, environmental effects for the short-lived fish were underestimated on average by 27.3% (lowest row of Fig. 4, Fig. S1 1 ), while in the remaining ageing error scenarios estimates were biased on average from -8.4% to +8.2% (Fig. S1 1 ; Table S1 1 ).

Ability of the model to reconstruct environmental time series
The correlation between time series of the year random effect (BLUP) extracted from the biochronological models and time series of the environmental variable used for the simulation of data increased with increasing sample size, with the greatest increase occurring from n = 1 to n = 10 fish per year (Fig. 6). Correlation decreased with the addition of an ageing error. The lowest mean correlation for the short-lived species was observed under the scenario "bias (-1)" and under the scenario "bias (-10%)" for the long-lived species (Fig. 6).
Note: Terms are defined in eq. 1; var = variance, corr = correlation. Fig. 2. Visualization of ageing errors introduced to the simulated growth data of the short-lived (a) and long-lived (b) species. Points indicate mean difference ±2 SD between simulated interpretation and the true age. Notice different x-axis ranges for the short-lived (a) and long-lived (b) species. CVcoefficient of variation between true and simulated age interpretation. Headings of the panels correspond to the different age error scenarios (see Materials and methods section for the detailed description).  Wavelet analysis showed poor coherence between growth biochronologies (BLUP of year random effect) and environmental time series when one fish per year was sampled (Fig. 7). The model fitted for the long-lived species under this sampling regime (n = 1 per year), and the scenario "perfectly aged" was able to recreate only low-frequency variations of the environmental time series. A gradual increase in the coherence was observed for both species with increasing sample size (Figs. S3 1 and S4 1 ). The addition of the ageing error significantly reduced the coherency of the highfrequency signals. The increase of sample size to ≥20 individuals per year resulted in partial coherence of the low-frequency signals (wavelet period > 16 years) between developed biochronologies and environmental time series (Figs. 7, S3 1 , and S4 1 ).
Sampling 10 correctly aged individuals of short-lived species per year resulted in a correlation between developed biochronologies and environmental time series of between R = 0.15 and R = 0.89, depending on the age group ranges tested (Fig. S5 1 ). Under all the ageing error scenarios, sampling of the oldest individuals contributed to the highest correlation between biochronological and environmental time series used in the test. For example, with the error scenario "bias (-1)", similar sampling of 10 individuals from the age group 3 only or age groups within the range 8-10 resulted in the mean correlation R = 0.15 and R = 0.36, respectively.

Discussion
Accurate prediction of the impacts of future environmental change on fish populations is largely dependent on our ability to quantify their climate sensitivity (Edwards et al. 2010). A growing number of studies have utilized the biological information naturally archived in otoliths to identify past environmental drivers of fish somatic growth (Mazloumi et al. 2017;Martino et al. 2019;Tanner et al. 2019). While many sclerochronological studies have successfully related interannual changes in fish growth to environmental factors based on otolith increment measurements (Morrongiello et al. 2012), there is a limited number of critical assessments of the robustness of inferences drawn from biochronological data (but see Black et al. 2016). Simulation-based studies can give useful insight into the possible shortcomings and strengths of analytical approaches (Bolker 2008;Kain et al. 2015;Allegue et al. 2017). Even though our simulations are specifically oriented to the analysis of fish otoliths, the results obtained can be   . 7. The wavelet coherency (colour gradient) and phase between biochronology obtained for short-lived (a) and long-lived (b) species and simulated environmental time series under selected ageing error scenarios and sampling schemes. Complete results are presented in the online Supplementary materials (Figs. S3 1 and S4 1 ). Time and wavelet periods (in years) are indicated on the x and y axes, respectively. The 5% significance level of coherence is shown with the thick black contour. When R 2 ≥ 0.7, the relative phase relationships are indicated with arrows (in-phase pointing right; the anti-phase pointing left; and lead of biochronology by 90°, pointing straight down). [Colour online.] of general importance for sclerochronological studies based on growth data collected from the hard parts of organisms.
Our simulations clearly identified that the development of a robust biochronology can be seriously hampered by insufficient sample size. Even in the case of perfectly aged individuals, the collection of one to five otoliths per year was not enough to obtain reliable estimates of environmental effects. The consequences of using fewer than ten individuals per year in the analyses were especially noticeable for the short-lived fish because of the small number of increments per individual relative to long-lived individuals. Moreover, a low number of otoliths limits the ability to properly recreate the variability of the environmental time series. The cross-wavelets analysis showed that the coherence of highfrequency signals in biochronology and environmental time series was significantly diminished when the number of otoliths was fewer than six per year. Thus, some minimum number of increment measurements per calendar year is necessary to discriminate population-wide growth pattern from individual-level variability. In dendrochronology, the minimum number of increments per calendar year is estimated through the expressed population signal, which is a measure of how well the sample set represents the theoretical, infinite population from which it was drawn (Wigley et al. 1984). Also, the related subsample signal strength indicates the extent to which a decreasing sample size back through time represents the full data set, assuming that the number of contributing samples diminishes back through time (Wigley et al. 1984). Both metrics are a function of the degree of correlation among samples (their synchrony or covariance) as well as sample size. Data sets with stronger synchrony among samples, as may be the case if environment is strongly limiting and growth-increment boundaries are clearly defined, require fewer samples to retain population-wide signals. Such criteria for defining minimum sample sizes must also be considered for biochronologies to ensure that environmental signals are fully retained along the length of the time series.
The introduction of even moderate ageing error (CV 5%-10%) into the biochronological data caused substantial underestimation of the environmental parameter in the mixed-effects models and thus decreases the probability of detecting an environmental signal when one is present. Falsely added or missed increments caused the frameshift of individual time series. Misaligned growth patterns between individuals can result in growth "peaks" cancelling out "troughs" and reduced synchrony among individual growth histories (Boehlert et al. 1989). Frameshifts also eliminate high-frequency growth variability in the biochronology, leaving only low-frequency variability, which can become offset in time if the degree of ageing error is sufficiently large . The loss of high-frequency signals is likely to increase with the longevity of the species (Boehlert et al. 1989;Black et al. 2016) because a single mistake can induce a frameshift that extends much farther back in time (Black et al. 2019). However, the ageing error of the "bias (+1, -10%)" scenario also caused noticeable systematic underestimation of the environmental parameters in the short-lived species.
Ageing errors can have different origins and appear in biochronological data by counting additional false increments (e.g., formed during major biological events and stresses) or by omission of true growth increments (e.g., due to their limited visibility) during visual investigation of the otolith structure (Panfili et al. 2002). While we did not investigate the specific nature of the ageing error, our results highlighted that even moderate ageing errors (CV 5%-10%) introduced in the random part of the individual growth time series can hamper the development of a robust biochronology, regardless of species longevity. The ageing errors introduced in the simulation are considered to be realistic, since a previous review of the literature across all ageing structures (otoliths, scales, rays, spines, and vertebrae, at both annual and daily scales) showed a median CV of age estimation of 7.6%, while the modal CV was 5% (Campana 2001).
Mixed-effects models overcome many of the shortcomings of other modeling techniques, such as fixed-effects linear models or detrending each measurement time series separately, by allowing for the proper partition of growth variation (Weisberg et al. 2010;Morrongiello and Thresher 2015). However, an important limitation of mixed models and other frequentist approaches is that they are not able to accommodate ageing errors in parameter estimation. Indeed, the majority of sclerochronological studies do not even quantify and report ageing error. Although Bayesian hierarchical models can incorporate ageing error (Shelton et al. 2013), statistical corrections cannot substitute for accurate and precise ageing procedures (Campana 2001). Careful preparation of correctly aged and dated growth measurements is the primary requirement prior to the sclerochronological analysis of organismal responses to environmental variability (Black et al. 2005).
One solution to eliminating or greatly reducing ageing error is use of cross-dating, which is a standard practice in tree-ring research and for which there is no substitute. If environmental variability influences growth, it will induce a synchronous, timespecific pattern among all individuals of a given species and location. Cross-dating is the process of matching these patterns among individuals, beginning with the known year of capture and working back through time. If an increment has been skipped or falsely added, the synchronous growth pattern will be offset in that individual, indicating that an error has occurred . The ability to apply cross-dating declines with decreasing species lifespan, given that the synchronous growth "bar code" has too few degrees of freedom to confidently date in short-lived (<10 years) individuals. However, cross-dating can still be useful, especially in helping to accurately date the increment at the otolith edge (Matta et al. 2010). A notable caveat to this approach is that the power of cross-dating to detect errors diminishes as the innermost otolith increments are approached. Thus, some ageing error may persist in the early years of life. Moreover, cross-dating is not possible across fish that experience different environmental regimes, as may be the case for migratory species ). Yet wherever it can be applied, cross-dating is a crucial step in the development of highly accurate biochronologies that fully capture environmental signals .
Precision and accuracy of estimated environmental effects was maximized by sampling younger fish under "no bias" and "bias (+1, -10%)" ageing error scenarios. However, our simulations also revealed that the mean correlation between developed biochronologies and the environmental time series was consistently maximized by sampling the oldest individuals. These findings emphasize that the precise and accurate quantification of average population-level responses to environmental variability, and the reconstruction of past environmental conditions, constitute two separate goals that can be accomplished with differently optimized sampling (Fritts 1976;Morrongiello et al. 2012;Klesse et al. 2018). To maximize ecological and evolutionary information within the biochronological data set and establish unbiased climate-growth relationships for the population, we should apply random selection of individuals (Fritts 1976;Klesse et al. 2018). In contrast, for the development of environmental reconstructions to hindcast climate prior to the start of the instrumental record, sampling of the individuals most likely to have sensitivity to the environmental variable in question would be more appropriate (Fritts 1976;Black et al. 2009;Butler et al. 2013;Reynolds et al. 2016). By intent, such sampling methods minimize ecological, physiological, and genetic sources of variation to maximize a common interannual environmental signal (Morrongiello and Thresher 2015).
"Regional curve standardization" is a technique used in dendrochronology (often for very long-lived individuals) in which all measured increments across all sampled individuals are aligned with respect to the age of formation and detrended with a single function (Briffa et al. 1992;Butler et al. 2010;van der Sleen et al. 2016), which is the equivalent of linear model with one fixed age slope for all fish (Weisberg 1993). In dendrochronological research, time series of growth can be also individually detrended using a mathematical function (e.g., negative exponential curves) to remove age-related and other long-term trends and then averaged to form a master chronology (Black et al. 2005). Such individual detrending corresponds to the mixed-effects model with random age intercept and slope for each fish (Weisberg et al. 2010). Year random effects BLUPs of these models are often used in sclerochronological studies to represent combined extrinsic environmental effects on fish growth. It has been shown that these year random effect BLUPs are equivalent to the biochronology developed in dendrochronological research (Helser et al. 2012;Matta et al. 2016). We used BLUPs to assess, in a simple and intuitive way, the association between the time series of our biochronology and environmental variables used to drive fish growth variation in the data simulations. However, it should be noted that BLUPs in the frequentist approach are point estimates, which can be drawn towards an overall average when a given year's increment sample size is low (see Hadfield et al. 2010 for further discussion). BLUPs are predicted with error, and their use in secondary analyses can be nonconservative because the error inherent in their prediction is excluded from these further tests (Houslay and Wilson 2017).
We assumed age of the long-lived fish to be uniformly distributed. We applied stationary (remaining the same across years) age distributions and the same parameters of growth model for both species. These simplifications made results more comparable and helped to isolate effects of ageing errors and sampling regimes on the sensitivity of environment-growth relationships derived from simulated otolith biochronologies. In reality, the age structure of the population can change through time under environmental or human-induced pressures (Ottersen 2008). Irregular exploitation due to differences in management regimes or market preferences, alterations in the fishing gear selectivity, or scientific program objectives could result in time-biased estimates of growth rates, which has implications for our ability to detect growth-environment relationships (Ricker 1969;Morrongiello et al. 2012). Since growth model coefficients (e.g., age parameter) characterize growth trajectories, they may also affect the reliability of the biochronological information. For example, strongly asymptotic growth (rapid growth during early life and very slow growth in the older ages) may limit the biochronological value of the older increments due to reduced width and elevated measurement errors. Since these effects were not investigated in our work, we encourage further simulation studies to test the potential consequences of varying age structures and growth model parameters on the accuracy and precision of the otolith biochronologies.
Here, we explored a simplified ecological system with one extrinsic variable affecting fish growth. Real-world ecological data sets often include many confounding and correlated variables (Legendre et al. 2002) and nonlinear relationships and complex interactions (Guisan and Zimmermann 2000). Moreover, contrasting environmental drivers of adult and juvenile growth (Ong et al. 2015), altered growth sensitivity to environmental conditions through time , regime shifts in the growth of organisms Smoliń ski and Mirny 2017), and regime-specific relationships (Lin and Petersen 2013) can substantially complicate the proper inference. It is therefore likely that our capacity to detect the environmentgrowth relationships by the application of the sclerochronological approaches in the real ecological systems can be even more diminished than in this simulation.
In conclusion, we show that realistic ageing error and commonly employed sample designs can cause a strong underestimation of environmental effects on fish growth in biochronological models and substantially reduce the coherence between a biochronology and associated environmental time series. These problems call for strict quality control of age data and usage of ageing validation techniques (Campana 2001), including cross-dating procedures whenever possible (Black et al. 2005. The use of short-lived species' otoliths, which are less biased by ageing errors, are beneficial when these errors are difficult or impossible to reduce. We also show that increasing sample sizes, up to at least ten individuals per year, can improve the precision of environmental parameter estimates (consistency of results from repeated estimations). Collection of more than ten individuals in a given year did not, however, result in substantially improved accuracy (proximity to the true value) of these estimates, which means that under ageing error, even a large sample size does not prevent biases. Age reading errors are naturally present in many biochronological data sets, and they may cause consistently biased estimates and diminish our ability to recognize true environmental effects. Such biases can affect more general conclusions on biological responses to environmental changes (e.g., in the metaanalytical studies ;Parmesan 2007;Poloczanska et al. 2013). For these reasons, reduction of ageing errors and high sample size in sclerochronological research are considered critical to reach robust conclusions on biological responses to environmental variability and generate accurate prediction of effects of global change.