Is telomere length a molecular marker of past thermal stress in wild fish?

Telomeres protect eukaryotic chromosomes; variation in telomere length has been linked (primarily in homoeothermic animals) to variation in stress, cellular ageing and disease risk. Moreover, telomeres have been suggested to function as biomarker for quantifying past environmental stress, but studies in wild animals remain rare. Environmental stress, such as extreme environmental temperatures in poikilothermic animals, may result in oxidative stress that accelerates telomere attrition. However, growth, which may depend on temperature, can also contribute to telomere attrition. To test for associations between multitissue telomere length and past water temperature while accounting for the previous individual growth, we used quantitative PCR to analyse samples from 112 young‐of‐the‐year brown trout from 10 natural rivers with average water temperature differences of up to 6°C (and an absolute maximum of 23°C). We found negative associations between relative telomere length (RTL) and both average river temperature and individual body size. We found no indication of RTL–temperature association differences among six tissues, but we did find indications for differences among the tissues for associations between RTL and body size; size trends, albeit nonsignificant in their differences, were strongest in muscle and weakest in fin. Although causal relationships among temperature, growth, oxidative stress, and cross‐sectional telomere length remain largely unknown, our results indicate that telomere‐length variation in a poikilothermic wild animal is associated with both past temperature and growth.


Introduction
The presence of environmental stress (i.e. the response to mostly unfavourable physical environmental variables) is central to both ecological and evolutionary research (Bijlsma & Loeschcke 2005;Wingfield & Boonstra 2013). Ecologically, physical environmental stress emerges at the boundaries of species-specific niches, and evolutionarily, it contributes to total phenotypic fitness. Thus, both ecologists and evolutionary biologists are interested in determining and quantifying environmental stress (Menge et al. 2002;Wingfield & Boonstra 2013;Schulte 2014). Whereas experimental laboratory settings allow quantification of contemporary stress by sampling individuals during direct stress exposure, stress often occurs in field-sampled individuals prior to sampling. Furthermore, stress is often quantified most meaningfully after stressful periods (Bijlsma & Loeschcke 2005). This poses a methodological challenge in determining and quantifying cumulative effects from previous stressful periods in both field and laboratory studies.
One major ecological parameter that can also become an environmental stressor, and which has increased in research focus during the recent decades, is environmental temperature. Temperature governs physical, chemical and biological reactions. In poikilothermic animals, increases in environmental temperature lead to increases in metabolic rates and, if food is unlimited, to accelerations in somatic growth (Allan & Castillo 2007). However, this potentially beneficial temperature effect on growth is expected only within certain species-specific ranges (Brett 1979). If the environmental temperature gets too high or too low, organisms may reach physiological limits (Fry 1947), and oxidative stress may occur as a consequence (Finkel & Holbrook 2000;P€ ortner 2002;Lushchak 2011). However, oxidative stress, which is defined as an imbalance between pro-and antioxidants that may disrupt intracellular signalling and damage proteins, lipids and DNA (Jones 2006), may also emerge under high somatic growth conditions (Finkel & Holbrook 2000;Carney Almroth et al. 2010;Pauliny et al. 2015). Accordingly, high environmental temperature may increase stress in poikilothermic animals either directly and/or indirectly via accelerated growth.
Useful biomarkers for past experienced stress may be the telomeres (von Zglinicki 2002;Houben et al. 2008;Monaghan 2014;Nussey et al. 2014;Bateson 2016). Telomeres are nuclear structures made up of noncoding DNA repeats and associated proteins that protect eukaryotic chromosome ends from accidental chromosome fusion as well as from loss of coding information during cell replication, in which chromosome ends are shortened due to the so-called 'end-replication problem' (Olovnikov 1973;Blackburn 2000). Every cell-replication cycle thus shortens the noncoding repetitive sequences of the telomeres; this process has been shown to limit the replication capacity of cells (Olovnikov 1973). However, repetitive telomeric sequences are also extremely sensitive to damage by reactive oxygen species, which are generated under oxidative stress (Oikawa & Kawanishi 1999). The extent of telomere shortening through oxidative stress damage is assumed to exceed that of telomere shortening through cell replication by somatic growth, even under mild but chronic oxidative stress (von Zglinicki 2000;Houben et al. 2008). The extent of past environmental and growth stresses may be an important determinant of telomere length, and telomere length may thus serve as a potential biomarker for quantifying those cumulative stress effects.
Compared to numerous studies that have focused on telomere-length dynamics imposed by growth or oxidative stress-elevating conditions in homoeothermic mammals or birds (reviewed by Houben et al. 2008;Monaghan 2014), similar studies in aquatic poikilothermic animals, such as fishes, are limited (reviewed by Bateson 2016). In contrast to mammals and birds, most fishes not only are poikilothermic (and therefore are much more affected by environmental temperature) but also are indeterminate growers (i.e. they exhibit lifelong somatic growth). Moreover, fishes, similar to many other poikilothermic animals, express telomerase (a ribonucleoprotein) across tissues and throughout all life stages (Klapper et al. 1998;Elmore et al. 2008;Hartmann et al. 2009;Lund et al. 2009;Anchelin et al. 2011). Telomerase is expressed only in some tissues and under specific conditions in homoeothermic animals, as compared with poikilothermic animals (reviewed by Gomes et al. 2010;Armanios & Blackburn 2012), and it is able to extend the telomeres and cell-replication capacity (Kim et al. 1994).
Even though telomerase is expressed throughout the entire life cycle of poikilothermic fishes, telomere attrition of many species has been reported after phases of rapid growth in early and late stages, though not at intermediate stages (Hatakeyama et al. 2008(Hatakeyama et al. , 2016Anchelin et al. 2011). Importantly, this life stage-specific telomere attrition in fishes is paralleled by life stagespecific telomerase expression variation, thus leading to life stage-specific telomere lengths (Hatakeyama et al. 2016). This phenomenon may explain, to some extent, why previous studies have yielded very inconsistent results with respect to growth-or age-related telomere attrition among strains of the same species (telomerelength constancy vs. shortening: Hartmann et al. 2009) or among organismal life stages of individuals from the same species (telomere-length extension vs. reduction : N€ aslund et al. 2015). An extreme case may be the much stronger telomere attrition observed in extremely rapidly growing fish transgenically expressing growth hormone, compared with their slow-growing wild-type maternal half-siblings . Finally, growth in fishes appears to affect tissues differently, similarly to observations in humans, and less severe effects are observed in the brain and heart relative to the muscle and blood (Hatakeyama et al. 2008).
Temperature effects on telomere attrition in poikilothermic vertebrates have, to our knowledge, been assessed in only two cross-sectional laboratory studies. A subtropical fish (Gambusia holbrooki, the Eastern mosquitofish) study has revealed shorter telomere length in individuals kept at a colder temperatures than the species-specific preferred range, but high-temperature stress effects were not tested (Rollings et al. 2014). A second study has found 15% shorter telomeres in 18 Siberian sturgeon (Acipenser baerii) individuals kept for 1 month at 30°C, 3°C above their normal temperature range, relative to 18 individuals kept at 20°C (Simide et al. 2016). Unfortunately, it remains unclear whether the latter study had confounding tank and temperature effects, because the presence of environmental replication was not reported. Together, the conditions under which growth affects telomere attrition in fishes have largely been inconclusive, and the suitability of telomeres as a biomarker for thermal stress in poikilothermic animals remains to be further investigated, especially whether telomere length might be used as a biomarker under natural conditions in which past growth and temperature effects may be confounded.
In this study, we tested for associations between cross-sectional relative telomere length (RTL) and summer water temperatures in 112 wild young-of-the-year (YOY) Salmo trutta L. (brown trout), a cold-water salmonid fish species, that were collected from 10 natural rivers. We tested the suitability of RTL as a biomarker for past environmental high-temperature stress while controlling for the association of RTL with the size of the YOY in autumn, which reflects the cumulative effects of rapid growth during their first summer. For both associations, we also tested for differences among six tissues (blood, fin, kidney, liver, muscle and skin). For poikilothermic animals (i.e. animals whose temperature covaries with the ambient temperature), we expected relatively uniform effects in response to environmental thermal stress because all tissues are affected similarly by the ambient temperature. In contrast, we hypothesized that the growth effects on the telomere length might differ among tissues because some tissues might either physiologically contribute to growth more than others (such as blood, kidney and liver relative to fin or skin) and/or experience stronger growth (such as muscle relative to skin or fin). Altogether, our study provides empirical evidence for the association between the telomere length, experienced past water temperature, past growth and for variation of these associations across six tissues in a wild poikilothermic aquatic species.

Materials and methods
All 10 sampled rivers drain into the Baltic Sea ( Fig. 1). To obtain river-specific temperature records, we planted loggers (HOBO 8K Pendant Temperature/Alarm Data Logger; Onset Computer Corporation) in nine rivers in May 2014, and for one river (Vodja), we obtained records from a hydro-meteorological station (Estonian Weather Service; http://www.ilmateenistus.ee/?lan-g=en). The loggers were set to record data directly at, or downstream of, the fish sampling sites at midnight and noon (Fig. 1).

Sampling
Between 11 and 14 September 2013, we caught 112 YOY individuals (9-13 per river) by electrofishing in 50-to 100-m river stretches (special fishing permit #82/2013; Estonian Ministry of Environment). We euthanized fish by an anaesthetic overdose with tricaine methanesulfonate (MS-222; Sigma-Aldrich) and measured individual body size (fork length AE 1 mm). Immediately thereafter, we collected blood samples by tail ablation with heparinized glass capillaries (0.5-0.6 mm diameter, Marienfeld), which we froze in liquid nitrogen. To obtain other tissue samples, we cut each individual dorsoventrally, anterior and posterior of the dorsal fin, with a scalpel. Whole cross-sections (containing tissues of the dorsal fin, skin, muscle, and kidney) were stored individually in 96% ethanol. We also dissected each liver, which we froze in liquid nitrogen.

DNA isolation
We extracted DNA for the quantitative PCR (qPCR) analyses. We used a commercial kit for the blood samples, following the manufacturer's protocol (DNeasy Blood and Tissue Kit; Qiagen), but added 190 lL of PBS to 10 lL of blood. For other tissues, we used the method described by Aljanabi & Martinez (1997), but we added 400 lL of NaCl instead of 300 lL after incubation, transferred 700 lL of the supernatant after centrifuging, added cold isopropanol instead of ethanol for washing and resuspended the mix in 50-200 lL of ddH 2 O. We assessed DNA purity and concentration spectrophotometrically (ND-1000; Thermo Scientific) and retained only samples with absorption ratios (260:280 nm) of approximately two. We standardized DNA concentrations with ddH 2 O to 2.5 ng/lL.

Quantitative PCR assay
We used qPCR to determine relative telomere length (RTL) as proposed by Cawthon (2002), which allows estimating telomere length as a total quantity across chromosomes relative to the quantity of a nuclear control gene. We employed primers for the telomere target (T) gene (Epel et al. 2004), which hybridize with vertebrate telomeric repeat sequences and have previously successfully been used in our study species (N€ aslund et al. 2015), and vertebrate 18S rRNA primers for the nuclear control (C) gene (Maeda et al. 2002). We performed amplifications using a 384-well 7900HT Fast Real-Time PCR system (Life Technologies). All reactions contained 10 lL 29 SensiFAST SYBR Hi-ROX mix (Bioline Reagents Ltd.), 400 nM of each primer, 4 lL template DNA and ultrapure PCR-grade H 2 O totalling 20 lL. Cycle conditions encompassed 5 min at 95°C, 40 cycles of 94°C for 10 s, 58°C for 15 s and 72°C for 20 s, followed by dissociation curve analysis. Both primer pairs amplified a single product (Appendix S2, Supporting information). We assigned samples randomly to one of 13 plates with both genes amplified on the same plate to reduce systematic biases. Six samples were included on all plates to reliably determine among-plate bias. Samples and no-template controls (NTCs; on all plates) were run in triplicate for each gene. Some NTCs resulted in unspecific amplifications, detected a minimum of 10 (T) or 15 (S) cycles later than samples. Due to the large differences in amplification cycles, we expect the effects of unspecific amplifications to be small in the targeted quantification (occurring in earlier cycles).

Statistical analysis
River temperature model. Temperature data were recorded between 31 May and 22 September 2014. Because some of the Preedi River data were compromised, we needed to construct a prediction for the 2014 Preedi temperatures: we modelled a time-series using the 2013 water temperature data from the closest studied river (Vodja), the 2013 Preedi water temperatures and the corresponding air temperatures available from a nearby weather station (25 km, V€ aike-Maarja, Fig. 1). The missing Preedi 2014 data were then predicted on the basis of the 2014 weather-station air and Vodja water temperature data (Appendix S1, Supporting information).
We estimated the differences in average summer water temperatures among the studied rivers on the basis of average daily water temperatures (for 2014). To this response of temperatures, we fit a linear mixed model (LMM) with fixed terms for rivers ('River'; our factor of interest), second-order orthogonal time polynomials (i.e. mean-centred linear and quadratic time effects: 'Time' and 'Time 2 ', modelling the major seasonal trend), interactions of time polynomials and rivers (modelling seasonal trend differences among rivers) and random terms for an autoregressive function of order one (ar1) across time ('Time.fac'; modelling the temporally correlated daily deviation from the smooth trends across rivers) and ar1-correlated residuals (modelling the temporally correlated, river-specific daily deviations from the common daily deviations), whose variance was conditionally fit for each river to provide a more reliable statistical testing (homo-vs. heteroscedastic residual variance for rivers: X 2 9 = 400.34, P < 0.001, range of s 2 = 0.20-2.38). The model was specified as follows (random terms, here and elsewhere, in italics): Quantitative PCR data processing and standardization. Based on raw qPCR-machine-recorded fluorescence data, we first estimated both plate-by-well-specific quantification threshold values (C q ) and amplification efficiencies (E; subsequently averaged, see below) using the Real-time PCR Miner (Zhao & Fernald 2005) online tool. The implemented algorithm estimates C q and E for each replicate from the kinetics of each amplification curve while accounting for noise during the exponential amplification phase (Zhao & Fernald 2005).
Before modelling the C q data, we removed technical outliers identified by melting temperature (T m ) values obtained by our dissociation curve analysis (Appendix S3, Supporting information). We found that pruning the T m -residual outliers prior to C q data modelling resulted in more consistent estimates of E and, therefore, facilitated subsequent C q data modelling (because C q is standardized by averaged E values, which are affected by outliers; see also Schefe et al. 2006). After the removal of the T m outliers, we calculated gene-byplate-specific E values (range: E (T) = 1.81-1.83, E (C) = 1.83-1.84). All C q values were first efficiency-standardized to a relative DNA copy quantity as E Cq (E being specific for each plate-by-gene combination) and subsequently standardized to a common log2-C q scale as log 2ðE Cq Þ With this method, E values do not need to be equal between genes or among plates, and all gene contrasts (C-T) represent the logE (standardized to E = 2) of the traditionally used T/C ratio, (E T /E C ) À1 (Cawthon 2002), that is gene contrasts in our study represented differences in E-standardized C q and thus RTL on the log2 scale. This approach is similar to the flexible and powerful method of Steibel et al. (2009) and allows jointly accounting for the technical and biological variations of both genes, following strongly recommended but outside of the agricultural sciences rarely implemented practices (Nussey et al. 2014).
Quantitative PCRaccounting for technical among-plate bias. Using the standardized C q data, we first estimated gene-specific plate effects (enabling a removal of this technical bias) on the basis of the subset of the six samples that were replicated on each plate. We used a LMM with fixed gene effects and random terms for plate-by-gene (our term of interest), individual, geneby-individual and well. This model also enabled the estimation of among-plate measurement precision (before removing the technical bias), which was calculated as the among-individual variance divided by the total variance per gene: intraclass correlation (ICC AE SE), ICC among = 0.96 AE 0.02 (T) and 0.98 AE 0.01 (C). Finally, we predicted plate-by-gene effects from this model, which we subtracted from all respective sample C q values to remove this technical bias (analogous to the method suggested by Pfaffl 2001 that uses ratios). Predicted plate effects were in the range of AE 0.1 C q for both genes.
Telomere-length model. We were interested in determining whether (i) RTL was associated with river-specific summer water temperatures or with individual body size, (ii) RTL differed among tissues and (iii) temperature or size associations differed among tissues. Accordingly, we modelled the response of standardized C q values jointly with the predictors of the results from both our summer-temperature model and field measurements on individual size. In this model, we included a fixed gene term, which we interacted with all other fixed model terms. The inclusion of these gene contrasts was necessary because we modelled RTL directly on the basis of C q values for both genes. Through this approach, only gene-contrast terms (C-T) represent RTL, whereas simple effects, in the absence of significant respective interactions, represent effects that are common to both genes (i.e. variation in nuclear DNA amount).
We thus fit an LMM with fixed terms for the two continuous covariates of river-specific summer temperature ('Temperature') and individual body size ('Size'; ln-transformed), and terms for the factors for the six tissues ('Tissue') and the two genes ('Gene'), and all possible interactions. Both continuous covariates were meancentred and variance-scaled to set intercepts at average values and to make estimates more meaningful. We also included random effects at several study-design levels to represent the two-phase 'multiple randomization experiment' that made it necessary to account for treatment design, blocking and randomizations at both the field and laboratory phases, a statistical necessity that is overlooked in many studies (Brien & Bailey 2006). We accounted for correlations at the river level and the associated among-tissue correlation by including random terms for river ('River') and river-by-tissue. We accounted for individual data correlations ('Individual') by adding an individual-by-tissue term that we fit with a 6 9 6 unstructured covariance structure among tissues to allow for (i) different among-individual variance for each tissue and (ii) different correlations between tissue pairs. Because we wanted to make inferences about RTL (i.e. the gene contrast), we also specified all random terms for the gene contrast. We accomplished this by recoding the gene term as binary ('Gene.bin'; C = 0, T = 1) that we interacted with all random terms. We also fitted a covariance between each gene-contrast term and either the river or the river-by-tissue term. In addition to plate effects, qPCR data can suffer from technical bias due to well effects (e.g. Eisenberg et al. 2015), which may be spatially correlated on the plate. Thus, we fitted this technical bias among wells as an exponential isotropic variance model (iexpv) for the 16 rows and 24 columns of the 384well qPCR plates (the C q values were spatially correlated; Appendix S4 and S6, Supporting information). Amongwell bias was in the range of AE0.24 C q and more than twice the among-plate bias range. We fitted the residuals, which reflected the technical among-well variation for the gene-by-individual-by-tissue combinations (i.e. the triplicates per gene) conditionally for each gene because the variance for the T gene was larger than for the C gene. The results of the model-testing procedure to select the above-described covariance structure are summarized in Appendix S6 (Supporting information). The model was specified as follows: C q ¼ Gene*Tissue*Temperature*Size þ Riverþ River : Gene:bin þ River : Tissueþ River : Tissue : Gene:bin þ Individual : Tissueþ Individual : Tissue : Gene:binþ Finally, we evaluated whether a relationship between RTL and river temperature might be directly confounded with a relationship between the average body size for each river and the river temperature. Hence, we fitted an LMM for the response of individual body size (ln-transformed) with the fixed river temperature covariate and a random river term.
Model fitting and outlier handling. All statistical analyses were conducted in R 3.0.2 (R Core Team 2013), and the modelling under residual maximum likelihood (REML) was conducted using ASREML-R v. 3.0 software (Butler et al. 2009). We first fitted models with a full-fixed structure. To identify erroneous data and to ensure the validity of model assumptions, random-term effects, including residuals, were evaluated for distribution assumptions, heteroscedasticity and presence of outlier effects. We defined outliers as studentized random effects obtained using the ASREML-R 'aom' function with Student's t-distribution at a false discovery rate <0.05 (FDR). For qPCR data, identifying and removing outliers is a recommended practice because many technical errors can occur during processing and measuring of samples (Burns et al. 2005;Schefe et al. 2006). Generally, qPCR outliers are defined as divergent measurements relative to the distribution of the majority of measurements, but 'divergent measurements' definitions often vary among studies (Burns et al. 2005). Our chosen approach identified outliers as conditional (i.e. after adjusting measurements for all other gene-specific effects in the model such as river, sample, tissue, plate and well) studentized measurements that deviate from the overall theoretical t-distribution of all conditional studentized measurements with a 0.05 probability of falsely identifying outliers. Thus, our approach accounted for (i) the change in the shape of the theoretical underlying distribution with sample size (the t-distribution converges with the normal distribution only under large sample sizes but has longer tails otherwise) and (ii) the increase in the probability of observing divergent measurements relative to the overall distribution with increasing sample size. In total, we removed 1.2% (60 measurements, including 13 data on an entire faulty well) of the total technical measurement values for T m of all samples (4830), 1.3% (six measurements) of the measurements for the threshold quantification values (C q ) of the samples that we replicated across plates (462) and 1.7% (79 measurements) of the C q measurements of all samples (4769). A re-analysis of the data with all outliers included resulted in the same inferences (Appendix S6, Supporting information). After removal of technical outliers and under the data-supported covariance matrix, we tested the fixed terms. We defined significance for fixed and random covariance terms as P < 0.05 and for random variance terms as P < 0.1. We evaluated random terms and covariance structures between nested models using likelihood-ratio tests (LRTs) and fixed terms using F-tests with denominator degrees of freedom approximated according to Kenward & Roger (1997). We removed nonsignificant fixed or random terms unless they represented the study design. We present all in-text estimates with standard errors that we approximated for correlations using the delta method.

Average summer temperature
We first quantified average water temperatures for the 10 rivers because testing for an association between RTL and water temperature would be meaningful only if water temperatures differed among rivers. The largest observed daily temperature differences among the 10 rivers were 10°C and occurred between mid-July and mid-August (Fig. 2a). The model-predicted average temperature across the investigated period differed up to 6.0°C among rivers (F 9/30.8 = 75.6, P < 0.001; Fig. 2b). Furthermore, water temperature exhibited different temporal curves among rivers (River:Time 2 ; F 9/ 44.8 = 5.8, P < 0.001), indicating that rivers warmed or cooled differently early and late in the year (Fig. 2a) but season-turning points were equal among rivers (River: Time; F 9/37 = 0.69, P = 0.715). Air temperature during the analysed period was on average 1°C warmer in 2013 (when we sampled) relative to 2014 (when we recorded water temperatures; Appendix S1, Supporting information).
Telomere length Next, we tested for associations between RTL and the model-predicted average water temperature, the body size, their interaction and for the differences in these associations among the six tissues. Most higher order interactions were nonsignificant and were removed from the model. In particular, the interaction effects of body size with temperature on RTL did not differ among or across tissues (Gene:Tissue:Temperature:Size; F 1/95.8 = 0.59, P = 0.706 and Gene:Temperature:Size; F 5/ 63.2 = 0.52, P = 0.475). Furthermore, the association of RTL with either temperature or body size did not vary among tissues (Gene:Tissue:Temperature and Gene:Tissue:Size). However, whereas the temperature-slope differences among tissues were clearly nonsignificant (Gene:Tissue:Temperature: F 5/34.5 = 1.64, P = 0.175), the size-slope differences among tissues were close to the significance threshold of 0.05 (Gene:Tissue:Size: F 5/ 98.8 = 2.15, P = 0.065), as was the average RTL among tissues (Gene:Tissue: F 5/39.4 = 2.17, P = 0.077). When we kept both the gene-by-tissue and the gene-by-tissue-by-size effects in the model to predict the nonsignificant size slopes for each tissue, the regression slope on the variance-standardized scale of RTL on size was steepest for muscle and gentlest for fin tissue (muscle: 0.18, blood: 0.15, liver: 0.13, skin: 0.13, kidney: 0.13, fin = 0.12).
Nonetheless, we found negative associations between the RTL and the average water temperature of the 10 rivers (Gene:Temperature) and the body size of the 112 individuals (Gene:Size; Table 1, Fig. 3a,b). When we used the modelled variance-standardized scale, the effect of water temperature on RTL was 0.226 AE 0.082 and that of body size was 0.157 AE 0.051. The backtransformed gene-by-temperature effect, adjusted to a common individual body size, was 0.123 AE 0.045 C q *°C À1 (this positive C q effect for the T gene suggests a negative RTL effect because RTL equals the C-T gene contrast). Thus, because we had standardized the C q to the log2 scale, fish in the warmest river were predicted to have 40% shorter telomeres than fish in the coldest river. The back-transformed gene-by-size effect, adjusted to a common water temperature, was 1.05 AE 0.34 C q *ln(mm) À1 , that is the largest studied fish (93.5 mm) was predicted to have 38% shorter telomeres (on the biological scale) than the smallest fish (48.5 mm). From the RTL model, we obtained high individual correlations between tissues for the fixed-effect conditioned RTL (range: r muscle,kidney = 0.924 AE 0.016 to r liver,kidney = 0.977 AE 0.006; Appendix S5 and S6, Supporting information).
We also tested whether a removal of the River and River:Gene terms affected the detection of the associations of RTL with temperature or body size. These terms provide the error terms for the significance testing of the RTL-temperature association and account for unspecified differences in the average RTL among rivers when estimating the RTL-size association. When the terms were removed, the significance of the RTL-temperature association greatly increased (F 1/113.0 = 36.76, P < 0.001) when tested based on individuals, whereas the significance of the RTL-size association disappeared (F 1/112.9 = 2.74, P = 0.10) when not controlled for riverspecific slope elevations. Results based on the additional model with body size as the response and water temperature as a continuous predictor indicated that the average fish size in the rivers and the river-specific water temperatures did not covary (F 1/8 = 1.52, P = 0.252; Fig. 3c).
Surprisingly, we also detected effects on C q common to both genes for tissues, temperature, size and interactions among the three terms (Table 1). This means that we probably detected effects resulting from the quantity of genomic DNA used in the qPCR assay. Although we had standardized total DNA quantity prior to the assay, these effects may underlie among-tissue differences in (a) (b) Fig. 2 Observed daily water temperature curves for the 10 studied rivers between June and end of September 2014 (a) and the corresponding predicted average water temperatures with 95% confidence intervals (b). Colours in (a) and (b) have been matched and scaled according to the temperatures in (b). df: degrees of freedom; ddf: denominator degrees of freedom; Gene: factor for the C-T gene contrast, equals RTL; Tissue: factor for six tissues of either blood, fin, kidney, liver, muscle or skin; Temperature: covariate for average temperature (scaled°C ) in each of the ten rivers; Size: covariate for individual fork length (scaled ln of mm). Bold values denote P < 0.05. genomic to mitochondrial DNA quantities (Leary et al. 1998) that may additionally vary with both body size and water temperature. These differences, however, did not affect the RTL estimates because they were standardized by a nuclear gene.

Discussion
Our empirical study of a poikilothermic wild salmonid fish demonstrates a negative association between multitissue RTL and both past river water temperature and individual body size for a common age. These results suggest that telomere length may be a useful biomarker of past environmental stress due to elevated temperatures and of internal stress due to rapid growth. The effect of water temperature on RTL was slightly larger than that of body size. Whereas RTL covaried with temperature among rivers, average body size among rivers did not covary with temperature among rivers, which indicated that direct confounding between water temperature and temperature-dependent growth was unlikely. Thus, we found no evidence for differences in temperature-induced growth among rivers that may have caused differential telomere attrition. Furthermore, the association of RTL with water temperature among rivers and the association of RTL with individual body size within rivers were estimated from a common multiple regression model (i.e. both as partial regressions) and at different hierarchical levels; we estimated the temperature trend across rivers at average body-size effects and the body-size trend within rivers (accounted for random river deviations from the temperature trend) at average temperature effects.

Telomere length and water temperature
The suggested link between telomere attrition and environmental stress across scientific disciplines (e.g. medical, ecological, evolutionary, and animal welfare research) is oxidative stress (reviewed by Houben et al. 2008;Monaghan 2014;Bateson 2016). Although oxidative stress levels across temperatures in natural settings remain to be quantified, strong biological evidence suggests that the summer peak temperatures in the warmest rivers did indeed bring about thermal stress (see below) and that this stress may have led to elevated oxidative stress. Oxidative stress has been shown to increase as a result of elevated water temperatures (reviewed by Lushchak 2011). However, oxidative stress has also been observed to increase at temperatures lower than the species-specific optimum (e.g. Vinagre et al. 2012), which may explain the shorter telomeres at lower temperatures found in mosquitofish (Rollings et al. 2014). These observations suggest that the relationship between telomere attrition rate and water temperature may be quadratic rather than linear.
Known biological temperature parameters for Salmo trutta indicate that the observed temperatures, which averaged up to 16.5°C during summer and peaked temporally at 23°C, greatly exceeded the preferred temperatures of 9-10°C for the studied age (reviewed by Elliott & Elliott 2010). On the basis of the recorded 2014 water temperatures and the presence of warmer air temperatures in 2013 relative to 2014, it is likely that peak temperatures in the warmest rivers in 2013 were even higher and, at least temporarily, were close to the incipient lethal temperature range (22-25°C; incurring 50% mortality after seven days; Elliott & Elliott 2010). Importantly, the recorded temperatures are likely to represent previously experienced temperatures because a potential spatial temperature escape appears to be unlikely in the investigated warmer rivers; the rivers are relatively small, water heats up in unshaded reservoirs upstream from our sampling sites (our unpublished results, see also Hayes et al. 2008), and we did not observe any deep pools where water may have cooled. In addition, although the congeneric Atlantic salmon (S. salar) seeks refuge from elevated temperatures in cooler pools, this behavioural response is limited to fish older than 1 year (Breau et al. 2007), and a similar age dependency can be expected for S. trutta. Recent anthropogenic influences on river temperature, such as removing canopy shading or building dams with top-layer discharge, have been found to result in water temperature increases comparable to the differences observed here (i.e. several degrees: Moore et al. 2005;Hayes et al. 2008). Similarly, the range of globally expected air temperature increase is 1-4°C until the end of this century (relative to 1986-2005; IPCC 2014). Moreover, extreme events, such as summer heat waves, are also expected to increase in frequency (Katz & Brown 1992). As a result, both present and future temperature increases, and their variation, are very likely to negatively affect salmonid populations through higher thermal stress.
Because we present field-study results, the detected association of RTL with water temperature might involve unmeasured confounding factors among rivers. We were able to substantially exclude confounding effects between water temperature and average growth in the rivers because body size and temperature did not covary and because the estimated growth effects within rivers were conditional on average temperature effects. Likewise, we estimated temperature effects among rivers conditional on average body size. However, the average age (defining zero age when first feeding is initiated), and not the size, among rivers may have differed by an unknown amount and may have led to a confounding between summer water temperature and average age among rivers. Unfortunately, we lack information about spawning periods (between autumn and early winter) in the year prior to sampling and about development time until the first feeding in spring. Therefore, the average age of the YOY in different rivers may have differed by several weeks (see below), and age, rather than growth or water temperature, may underlie RTL variation among rivers. Nonetheless, we believe that the probability is low that the age ranking for YOY among the ten rivers, which is a compound effect of differences in spawning time and developmental time until first feeding, was correlated with the average summer temperatures.
As an inevitable consequence of analysing field data, water temperature and geography are confounded in our study, and this may have introduced a variety of unknown additional environmental factors that covaried with both geography and water temperature. For example, the environmental stress exerted by other factors, such as toxins, heavy metals and other pollutants, may increase with water temperature (Heugens et al. 2008), and this presence of unknown environmental factors could make water temperature a potential amplifier but not the direct cause of the detected association, which was based on quantified cumulative total stress effects. Likewise, genetic effects on telomere-length quantification such as variation in interstitial telomeric sequences, which get co-amplified by qPCR, may have differed among the populations. Nonetheless, our inferences were not based on an explorative data-mining analysis of many environmental variables whose statistical significance might arise by chance when not controlling for multiple testing (Parker et al. 2016). Instead, we tested for an association between the RTL and a selected environmental variable (water temperature) under an a priori hypothesis, and this reduced the probability of water temperature being a falsely identified environmental variable that associates with RTL.
Alternatively, we may have detected possible longterm effects of the temperature regimes rather than effects of contemporary, cohort-specific stress. In the yeast Saccharomyces cerevisiae, a transgenerational elevated temperature regime (37 vs. 30°C) has been shown to cause the strongest telomere attrition among 13 environmental stressors (Romano et al. 2013). Nevertheless, because we did not assess RTL before and after summer (i.e. longitudinally), we were unable to differentiate between transgenerational effects of river temperature regimes and cohort-specific past thermal stress. The potential presence of confounding environmental effects and the uncertainty of underlying (transgenerational) mechanisms for the detected RTL emphasize the advantage of carrying out longitudinal studies that may help alleviate some of these methodological problems (Price et al. 2013). However, the necessity of lethal sampling for multiple-tissue sampling precluded a longitudinal study. In our study, we can assume that the ancestors of the sampled individuals experienced comparable temperature differences every year; S. trutta has been shown to exhibit strong natal homing (Keefer & Caudill 2013). Accordingly, the detected effects of water temperature differences on RTL may also hold even if there are underlying transgenerational, long-term effects. Thus, although it was not possible to disentangle transgenerational and cohort-specific effects in our cross-sectional field study, our results suggest that RTL may be a useful biomarker of past thermal environmental stress in poikilothermic animals.

Telomere length and body size
Even though telomerase is expressed across tissues in all life stages in poikilothermic animals, telomere attrition appears to occur under rapid growth. Rapid growth is exhibited during early stages (Hatakeyama et al. 2008(Hatakeyama et al. , 2016 and is accompanied by high oxidative stress (as reported for our study species by Carney Almroth et al. 2010), which may explain why we detected an association between RTL and within-cohort body size. It has previously been shown that telomere attrition may be confined to early stages when it is usually greatest in many organisms (reviewed by Bateson 2016) such as young rice fish (Hatakeyama et al. 2008(Hatakeyama et al. , 2016, and it probably applies to YOY brown trout studied here. It remains to be tested whether temperature effects are also more likely to be detected at early stages when telomerase activity may not keep up with telomere attrition. Our conclusions may be limited because we cannot exclude the possibility that body-size differences within rivers may also underlie small within-cohort age differences. The brown trout spawning period, which is generally between autumn and early winter, can last for several weeks (Elliott 1984), and therefore, age among sampled individuals within and among rivers may have differed by several weeks. Thus, age (i.e. the duration of past growth) rather than magnitude of past growth rate may have contributed to the detected relationship between RTL and body size, although differential growth rates among families (underlying both genetics and environment) between the first feeding in spring and the September sampling probably blurred or even eradicated age-related size differences.
Importantly, our study results show that an appropriate study design and the statistical account of environmental and growth effects are necessary to provide reliable results. In particular, when omitting the random river terms from our model, which adsorbs unknown environmental and genetic sources of RTL variation among rivers that are not accounted for by the included fixed effects, and regarding individuals as experimental units (an often falsely made assumption in [molecular] ecological studies ;Hurlbert 1984;Meirmans 2015), effects from past growth (or age) were not detectable, and thermal effect significance was greatly overestimated. Nonetheless, our results indicate that past elevated environmental water temperature may represent an additional and more important factor than past growth for predicting YOY telomere length in this cold-water fish.

Differences among tissues
The association between RTL and average past water temperature was statistically similar among the six tissues, thus suggesting that temperature stress in poikilothermic animals is largely nonspecific with respect to tissue. However, the association of RTL with cohortspecific body size appeared (though nonsignificantly) strongest in muscle and weakest in fin (even though both were preserved and processed similarly). Thus, it is possible that stress resulting from rapid growth may be stronger in muscle. This latter (nonsignificant) trend has practical implications because fins, which exhibited the weakest RTL-body-size association, appear to be the most common tissue used in ecological or evolutionary fish studies on telomere dynamics (e.g. N€ aslund et al. 2015; Pauliny et al. 2015). Our results indicate that fin clips, which can easily be sampled nonlethally, may still be useful in studies on thermal stress but that growth effects on telomere attrition may be strongest in muscle, which can also be sampled nonlethally but requires certain biopsy techniques (e.g. Van Meter 1995).
In summary, we detected an association between cross-sectional telomere length and both past experienced water temperature and past realized growth. Both relationships may indicate that telomere attrition may reflect a causal relationship between telomere attrition rate and temperature-or rapid-growth-dependent (but age-specific) past oxidative stress. Thus, telomere attrition appears to also occur in wild poikilothermic species that express lifelong telomerase across tissues. Finally, an RTL constancy across tissues appears to hold across temperature conditions but probably not across growth conditions (which tended to affect muscle more than other tissues). Given that many known short-and long-term anthropogenic activities might cause water temperature increases of several degrees, anthropogenic actions may negatively affect poikilothermic aquatic animals through increases in high-temperature-induced environmental stress. Future studies to quantify this thermal stress may benefit from longitudinal assessments of (growth-adjusted) telomere attrition rate. editor for greatly improving earlier versions of the manuscript. This study was supported by Estonian Ministry of Education and Research institutional research project funding to AV (IUT8-2), Academy of Finland grants to PI (135653, 272713) and AV (266321) and a German Research Foundation research fellowship to PVD (DE 2405/1-1). The authors declare no conflict of interest.

Supporting information
Additional supporting information may be found in the online version of this article.
Appendix S1 Water temperature reconstruction for the Preedi river.
Appendix S2 Dissociation curves for the two qPCR assays.
Appendix S3 Melting temperature outlier model.
Appendix S4 Well effects of the qPCR assay.
Appendix S5 RTL correlations between tissues.

Fig. S1
Observed summer air temperature curves for summer 2013 and 2014 (a); observed and predicted water temperature curves for Preedi River and observed water temperature curve for Vodja River for summer 2013 (b); and observed water temperature curve for Vodja River or predicted water temperature curve for Preedi River with 95% confidence intervals for summer 2014 (c).

Fig. S5
Mixed-model-predicted well effects for threshold quantification cycle [DC q ; (a)], estimated across 13 qPCR plates (384 wells per plate), which were measured on a 7900HT Fast Real-Time PCR system (Life Technologies).

Fig. S6
Pairwise relationships between the six tissues for relative telomere length (RTL; tissues on y vs. x indicated in each panel header).
Table S1 Results of the model testing via likelihood ratio tests to determine the covariance structure of the mixed model on C q values to quantify relative telomere length (RTL) and test specific fixed effects.

Table S2
Variance (s 2 ), correlation (r), or isotropic autocorrelation parameter (/) estimates with REML standard errors (in parentheses) for the random terms of the model for C q values to determine relative telomere length (RTL) and test specific effects. Table S3 ANOVA table for results of the full mixed model on Cq values for two genes to determine relative telomere length (RTL) and test specific effects among 10 rivers and 112 individuals.

Table S4
For comparison to what happens when the qPCR measurement data has not been purged of outliers, i.e. when amplification efficiencies are based on all data and all Cq data have been kept in the model, we re-run the RTL model and provide these additional results.

Table S5
For comparison to what happens when the qPCR measurement data has not been purged of outliers, i.e. when amplification efficiencies are based on all data and all C q data have been kept in the model, we re-run the RTL model and provide these additional results.