Genetic architecture and heritability of early‐life telomere length in a wild passerine

Early‐life telomere length (TL) is associated with fitness in a range of organisms. Little is known about the genetic basis of variation in TL in wild animal populations, but to understand the evolutionary and ecological significance of TL it is important to quantify the relative importance of genetic and environmental variation in TL. In this study, we measured TL in 2746 house sparrow nestlings sampled across 20 years and used an animal model to show that there is a small heritable component of early‐life TL (h2 = 0.04). Variation in TL among individuals was mainly driven by environmental (annual) variance, but also brood and parental effects. Parent‐offspring regressions showed a large maternal inheritance component in TL ( hmaternal2 = 0.44), but no paternal inheritance. We did not find evidence for a negative genetic correlation underlying the observed negative phenotypic correlation between TL and structural body size. Thus, TL may evolve independently of body size and the negative phenotypic correlation is likely to be caused by nongenetic environmental effects. We further used genome‐wide association analysis to identify genomic regions associated with TL variation. We identified several putative genes underlying TL variation; these have been inferred to be involved in oxidative stress, cellular growth, skeletal development, cell differentiation and tumorigenesis in other species. Together, our results show that TL has a low heritability and is a polygenic trait strongly affected by environmental conditions in a free‐living bird.

Several studies have investigated the potential of telomere dynamics (i.e., individual differences in TL and telomere loss rate) in mediating life-history trade-offs both across (Dantzer & Fletcher, 2015;Pepke & Eisenberg, 2020) and within relatively long-lived species (Monaghan, 2010;Spurgin et al., 2018). However, despite being an ecologically important trait in many species , knowledge about the genetic architecture of TL and its adaptive potential in wild populations remains scarce .
Quantifying the additive genetic variance of a trait is required to understand mechanisms driving adaptive evolution, that is, the response to selection on a trait (Ellegren & Sheldon, 2008;Kruuk et al., 2008;Lande, 1979). However, the magnitude of the heritability and mode of inheritance of TL is not well-known in populations of wild animals, and few general patterns have been described (Bauch et al., 2019;Horn et al., 2011). Utilizing longterm pedigree data, individual variation in early-life TL can be decomposed into various genetic and environmental sources of variation through a type of mixed-effect model ("animal model"), which takes all relationships from the pedigree into account (Kruuk, 2004;Wilson et al., 2010). Estimates of TL heritabilities from studies using animal models (reviewed in ) have varied considerably across wild bird populations, from h 2 = 0 (n = 177, in whitethroated dippers, Cinclus cinclus, Becker et al., 2015) to h 2 = 0.74 (n = 715, in western jackdaws, Coloeus monedula, Bauch et al., 2021).
While most studies are characterized by relatively small sample sizes, recent long-term studies on Seychelles warblers (Acrocephalus sechellensis, n = 1317, h 2 = 0.03-0.08, Sparks et al., 2021) and common terns (Sterna hirundo, n = 387, h 2 = 0.46-0.63, Vedder et al., 2021) also revealed contrasting estimates of TL heritabilities. Epidemiological studies of humans have documented consistently high TL heritabilities, ranging from h 2 = 0.34-0.82 . In humans, some studies reported strong paternal inheritance (e.g., Njajou et al., 2007) or maternal inheritance (e.g., Broer et al., 2013) or that there were no differences in parental mode of inheritance (e.g., Eisenberg, 2014). In birds, several studies have documented maternal effects on offspring telomere dynamics (Asghar et al., 2015;Heidinger et al., 2016;Horn et al., 2011;Reichert et al., 2015), or effects of parental age at conception on offspring TL (Eisenberg & Kuzawa, 2018;Marasco et al., 2019;Noguera José et al., 2018). Reichert et al. (2015) found a significant correlation between mother-offspring TL measured at 10 days of age in king penguins (Aptenodytes patagonicus), but not when TL was measured at later ages (>70 days). This may be because post-natal telomere loss rate is strongly influenced by individual environmental circumstances (Chatelain et al., 2020;Wilbourn et al., 2018) and does not always correlate strongly with chronological age (Boonekamp et al., 2013(Boonekamp et al., , 2014. Faster growth in early life is associated with reduced longevity (Metcalfe & Monaghan, 2003) and TL may be involved in mediating the trade-off between growth rate and lifespan (Salmón et al., 2021;Young, 2018). Accordingly, a negative phenotypic correlation between TL and body size or growth rate has been documented within several species (Monaghan & Ozanne, 2018, but see Boonekamp et al., 2021. Telomeres are known to shorten during growth (Ringsby et al., 2015), but a negative phenotypic correlation may also indicate the existence of a negative genetic correlation (Roff, 1995;Roff & Fairbairn, 2012). Froy et al. (2021) reported a modest negative genetic correlation (r A = -0.2) between bodyweight and TL in feral Soay sheep (Ovis aries). Furthermore, we have previously shown that artificial directional selection on body size in wild house sparrows (Passer domesticus) affected TL in the opposite direction . This suggests that there is a genetic correlation between the two traits. Thus, quantifying the genetic correlation between TL and body size enables us to determine whether the two traits can evolve independently of each other or if the pattern of selection on both traits is needed for predicting evolutionary responses (Kruuk et al., 2008).
Telomere length is a complex phenotypic trait (Aviv, 2012;Hansen et al., 2016) expected to be polygenic, that is, affected by small effects of many genes Hill, 2010). Accordingly, numerous genome-wide association studies (GWAS), which tests correlative associations of single-nucleotide polymorphisms (SNPs) with specific traits, have identified several loci correlated with TL in humans that map to genes involved in telomere and telomerase maintenance, DNA damage repair, cancer biology, and several nucleotide metabolism pathways (e.g., Andrew et al., 2006;Codd et al., 2010Codd et al., , 2013Coutts et al., 2019;Deelen et al., 2013;Delgado et al., 2018;Jones et al., 2012;Levy et al., 2010;Li et al., 2020;Liu et al., 2014;Mangino et al., 2012Mangino et al., , 2015Mirabello et al., 2010;Nersisyan et al., 2019;Ojha et al., 2016;Soerensen et al., 2012;Vasa-Nicotera et al., 2005;Zeiger et al., 2018). None of the GWA studies in humans specifically tested the marker associations of early-life TL, which pose a challenge to the interpretation of the results, as TL shortens through life in humans (Blackburn et al., 2015) and genes may have different impacts at various life stages (Weng et al., 2016). Furthermore, large sample sizes and dense sampling of genetic loci is needed to ensure high power in GWA studies (Mackay et al., 2009) and resolve any pleiotropic effects (Prescott et al., 2011). The genes influencing TL in humans that were identified through GWAS only explain a small proportion of the interindividual variation in TL (<2%, Aviv, 2012;Codd et al., 2013;Fyhrquist et al., 2013). One GWAS on TL of a nonhuman species (dairy cattle, Bos taurus) was recently performed (Ilska-Warner et al., 2019) supporting the polygenic nature of early-life TL. However, domesticated species in captivity may display TL dynamics that are not representative of natural populations (Eisenberg, 2011;Pepke & Eisenberg, 2021). There are to the best of our knowledge no previous GWAS on TL performed in natural populations.
In this study, we aimed to provide novel insights into the genetic architecture of TL and the evolutionary mechanisms by which natural selection can alter telomere dynamics using data from a passerine bird. We obtained a single measure of TL in individuals (n = 2746) born within 20 cohorts in two natural insular populations of wild house sparrows at a similar age (c. 10 days), in addition to individuals at the same age in two insular populations that underwent artificial selection on body size for four consecutive years (n = 569, Kvalnes et al., 2017;. First, we estimated the phenotypic correlations between TL and tarsus length (as a proxy for body size, Araya-Ajoy et al., 2019;Senar & Pascual, 1997) in house sparrow nestlings. Second, we tested for effects of parental age on offspring TL. Third, we estimated heritability, environmental variances, and parental effects on early-life TL, and test for genetic correlations between TL, tarsus length, and body condition in the natural populations (primary analyses). Nestling body condition (body mass corrected for structural body size, Schulte-Hostedde et al., 2005) is included here to account for the component of body size that is not explained by tarsus length, which could be due to variation in the mass of other tissues or fat reserves (Peig & Green, 2010). We then used similar analyses in the artificially selected populations to validate our results from the primary analyses. Finally, we used highdensity genome-wide single nucleotide polymorphism (SNP) genotype data (Lundregan et al., 2018) in a GWAS to identify genetic regions and potential candidate genes underlying variation in earlylife TL within wild house sparrows (up to n = 383).

| Study populations and data collection
The study was performed in four insular house sparrow populations off the coast of northern Norway ( Figure S1.1 in Appendix S1).
The study periods differed between the populations with data from   and were analysed separately in a set of secondary analyses as replications of the primary analyses. All four islands are characterized by heathland, mountains, and sparse forest. The sparrows live closely associated with humans and within the study area they are found mainly on dairy farms (Hestmannøy, Vega and Leka), where they have access to food and shelter all year, or in gardens and residential areas (Traena), where they may be more exposed to weather conditions (Araya-Ajoy et al., 2019). Natural nests inside barns or artificial nest boxes were visited at least every ninth day during the breeding season (May-August) to sample fledglings (5-14 days old, with a median of 10 days). All individuals were ringed using a unique combination of a metal ring and three plastic colour rings. Fledged juvenile sparrows and unmarked adults were captured using mist nets from May to October. These procedures ensured that approximately 90% of all adult birds were marked on all islands during the study period (Jensen et al., 2008;Kvalnes et al., 2017). We measured tarsometatarsus (tarsus) length using digital slide calipers to nearest 0.01 mm and body mass to nearest 0.1 g with a Pesola spring balance (see Appendix S1). Morphological measurements were taken by different fieldworkers. All fieldworkers were carefully trained to consistently use the same measurement technique of THR or, in some cases, another experienced fieldworker (Kvalnes et al., 2017). For 234 out of 2746 nestlings, no nestling morphological measurements were made. Following Schulte-Hostedde et al. (2005) nestling body condition was calculated as the residuals of a linear regression of mass on tarsus length (both log 10 -transformed). To avoid collinearity in models where both nestling age and tarsus length were included as covariates, we age-corrected tarsus length by using the residuals from a regression of tarsus length on age and age squared (to account for the diminishing increase in tarsus length with age).
One blood sample (25 μl) was collected from each fledgling, which was stored in 96% ethanol at room temperature in the field and subsequently at -20°C in the laboratory until DNA extraction.

| Molecular sexing and pedigree construction
DNA extraction is described in Appendix S1. Sex of most fledglings (n = 2641) was determined using amplification of the CHD-gene located on the avian sex chromosomes as described in Griffiths et al. (1998). 21 individuals were sexed exclusively based on their phenotype as adults and 84 nestlings could not be sexed. The pedigree construction is detailed in previous studies (Billing et al., 2012;Jensen et al., 2003Jensen et al., , 2008Rønning et al., 2016). Briefly, we used individual genotypes on 13 polymorphic microsatellite markers scored using the genemapper 4.0 software (Applied Biosystems) to assign parentage in cervus 3.0 (Kalinowski et al., 2007). Nestlings within the same clutch were assumed to have the same mother. Nestlings with missing parents (unassigned: n = 662 with missing mother and n = 700 with missing father) were assigned dummy parents, assuming that nestlings within the same clutch were full siblings and thus had the same (unassigned) parents. The dummy parents were included in the pedigree as founders. We calculated individual inbreeding coefficients (F) based on the microsatellite pedigree using the R package pedigree (Coster, 2012). Pedigrees were ordered using the R package MasterBayes (Hadfield et al., 2006) and pruned to only contain informative individuals. The pruned pedigrees included 4118 individuals (3093 maternities and 3130 paternities) in the natural populations, and 1057 individuals in artificially selected populations.
Maximum pedigree depth was 13 generations, the number of equivalent complete generations (the sum of the proportion of known ancestors across all generations, Wellmann, 2021) was 1.510, and mean pairwise relatedness was 0.003.

| Telomere length measurements
Relative erythrocyte telomere lengths (TL) of 2746 nestlings from Hestmannøy and Traena (sample sizes are detailed in Table S1.1 in Appendix S1) were successfully measured using the real-time quantitative polymerase chain reaction (qPCR) amplification method by Cawthon (2002) with modifications by Criscuolo et al. (2009). Primer sequences, PCR assay setup and thermal profiles followed  and are detailed in Appendix S1. Briefly, this method measures the ratio of telomere sequence relative to the amount of a nonvariable gene (GAPDH) and a reference sample. The reference sample consisted of pooled DNA from six individuals, which was also included as a two-fold serial dilution (40-2.5 ng/well) on all plates to produce a standard curve, in addition to a nontarget control sample (all in triplicates). Samples were randomized and run on 2 × 125 96well plates (telomere and GAPDH assays, respectively). The qPCR data was analysed using the qbase software (Hellemans et al., 2007), which computes relative TL as the ratio (T/S) of the telomere repeat copy number (T) to a single copy gene number (S) similar to Cawthon (2002). In qBASE the T/S ratio is calculated as calibrated normalized relative quantities (CNRQ) that control for differences in amplification efficiency between plates and for inter-run variation by including three inter-run calibrators from the standard curve. All individual plate efficiencies were within 100 ± 10% (mean telomere assay efficiency was 97.5 ± 3.9%, and 97.6 ± 4.2% for GAPDH assays).
The average of the reference sample cycle thresholds (Ct) across all plates were 10.54 ± 0.03 SD and 21.53 ± 0.02 SD for telomere and GAPDH assays, respectively. Thus, while reproducibility of TL measurements within the reference sample of the same DNA sample extract is high, we performed DNA re-extraction of the same blood samples for 25 individuals to test TL consistency across DNA extractions (Appendix S1). The re-extractions were run on different plates and the TL estimates of these samples remained highly correlated (R 2 = 0.75, Figure S1.3 in Appendix S1). For these individuals, the average of the TL measurements was used in subsequent analyses. All reactions for the primary analyses (from the populations on Hestmannøy and Traena) were performed by the same person (MLP).
MLP and WB generated the secondary data set (n = 569 on 2 × 21 plates, from the populations on Leka and Vega) as described in . The primary and secondary data sets used different reference samples and are therefore not combined in the analyses.

| The correlation between tarsus length and telomere length
We first tested the phenotypic correlation between TL and tarsus length (as a proxy for body size) within 2462 house sparrow nestlings from Hestmannøy and Traena. TL (response variable) was log 10 -transformed and linear mixed-effects models (LMMs) were fitted with a Gaussian error distribution (R package lme4, Bates et al., 2015). Sex differences in TL are known for house sparrows . Thus, models included sex (continuous) fledgling age at sampling, hatch day (numbered day of year mean centred across years), and island identity as fixed effects. We fitted random intercepts for brood identity, year, and qPCR plate identity to account for the non-independence of nestlings from the same brood, year and plate. Because our study populations are known to be affected by inbreeding depression , we included the inbreeding coefficient (F, continuous) as a fixed effect (Reid & Keller, 2010). We then compared models with and without (age-standardized) tarsus length using Akaike's information criterion corrected for small sample sizes (AICc, Akaike, 1973;Hurvich & Tsai, 1989), and Akaike weights (w) and evidence ratios (ER) to determine the relative fit of models given the data (Burnham & Anderson, 2002). Models were validated visually by diagnostic plots and model parameters are from models refitted with restricted maximum likelihood (REML). Estimates and 95% confidence intervals (CI) are reported.
We applied within-subject centring (van de Pol & Wright, 2009) to separate within-parental age effects (e.g., senescence) from between-parental age effects (e.g., selective disappearance), by including both the mean parental age at conception and the deviation from the mean parental age for each parent as fixed effects in two LMMs (for fathers and mothers, respectively) explaining variation in offspring TL (log 10 -transformed). Both models included island identity and sampling age as fixed effects, and random intercepts for year, qPCR plate identity, and either maternal identity or paternal identity.

| Heritabilities and genetic correlation of telomere length, tarsus length, and body condition
We used a multivariate Bayesian animal model (Hadfield, 2019;Kruuk, 2004)  were fitted with a Gaussian error distribution using the R package MCMCglmm (Hadfield, 2010). Models included sex, fledgling age at sampling (associated only with TL and condition), island identity, and inbreeding coefficient (F) as fixed effects (Wilson, 2008), which were fitted such that different regression slopes were estimated for each trait (Hadfield, 2019). To estimate variance components, random intercepts were included for individual identity linked to the pedigree ("animal", V A ), brood identity (V B ) nested under mother identity, father (V F ) and mother identity (V M ), and birth year (cohort effects, V Y ). Parental effects include those influences on offspring TL that are repeatable across the lifetime of the mother or father (Kruuk & Hadfield, 2007), while brood identity accounts for other common environmental effects (McAdam et al., 2014). House sparrows are multibrooded laying up to three clutches in a season and may breed in multiple years, with an average of 3.6 ± 1.3 SD fledglings per brood in this study. Furthermore, to account for variance associated with measurement error we included qPCR plate identity (V O , associated only with TL, see e.g., Froy et al., 2021;Sparks et al., 2021). Random effects were generally specified with 3 × 3 covariance matrices to estimate the variances and covariances between the effects for each trait. We used inverse-Wishart priors for random effects and residual variances in the multivariate model (V = I 3 and nu =3, Hadfield, 2019). We reran analyses with other relevant priors (parameter expanded) to verify that results were not too sensitive to the choice of prior. The MCMC chain was run for 2,000,000 iterations, sampling every 500 iterations after a burnin of 5% (100,000 iterations). Mixing and stationarity of the MCMC chain was checked visually and using Heidelberger and Welch's convergence test (Heidelberger & Welch, 1983) implemented in the "coda" package (Plummer et al., 2006). All autocorrelation values were <0.1 and effective sample sizes were >3000. The narrowsense heritability was calculated as the posterior mode of the proportion of phenotypic variance (V P ) explained by additive genetic variance : We also estimated heritabilities excluding V O from the total phenotypic variance since it does not represent biological variance (de Villemereuil et al., 2018). Estimates are provided as their posterior mode with 95% highest posterior density intervals (HPD). All analyses were performed in R version 3.6.3 (R Core Team, 2020).
We also ran univariate models of TL, tarsus length and body condition including the same fixed and random effects as in the multivariate model (Appendix S2). For comparison with previous studies (e.g., Asghar et al., 2015), we tested whether maternal TL and/or paternal TL predicted offspring TL using two LMMs (parent-offspring regressions, Appendix S2). Parental heritabilities ( h 2 maternal and h 2 paternal ) can be estimated from parent-offspring regressions as the slope multiplied by two (one sex contributes half of the genes to their offspring).
We used the R package pedantics  to show that, based on parent-offspring regression, the pruned pedigree of the natural populations had ≥80% power to detect heritabilities ≥0.21 (see Figure S1.2 and Appendix S1). Furthermore, we estimated maternal (V DAM ) and paternal (V SIRE ) genetic effects (e.g., Wolf & Wade, 2016) in a multivariate animal model by fitting random intercepts for maternal and paternal identity linked to the pedigree to quantify these effects while accounting for the environmental variances specified above (Appendix S2). Maternal and paternal heritabilities were calculated as: (Wilson et al., 2005). To test for sex-specific heritabilities (e.g., Jensen et al., 2003;Olsson et al., 2011), we ran a bivariate animal model of TL in females and males as two different phenotypic traits with a genetic correlation between them (Appendix S2).

| SNP genotype data and association analyses
Nestlings that survived to adulthood (recruited) on Hestmannøy and Traena were genotyped on a high-density 200K SNP array (detailed in Lundregan et al., 2018) with median distances between SNPs shorter than 5000 bp. SNPs were originally identified from whole-genome resequencing of 33 individual house sparrows which were mapped to the house sparrow reference genome (Elgvin et al., 2017). DNA was extracted as described in Hagen et al. (2013), separately from telomere analyses. Data preparation and quality checks were performed using the GenABEL package (GenABEL project developers, 2013). We removed SNPs or individuals for which there was more than 5% missing data, the minor allele frequency (MAF) was less than 1%, or pairwise identity-bystate (IBS) was more than 95%. After quality control, the genomic relationship matrix (GRM) was computed based on 180,650 (180,666) autosomal markers in 373 (383) (Rönnegård et al., 2016): The first model included age-standardized tarsus length as a covariate, and the second model did not. Both models included sex, age, hatch day (mean centred), F, and island identity as fixed effects, and brood identity, year, qPCR plate, and the GRM fitted as random effects. We estimated the proportion of phenotypic variance explained by each SNP as: h 2 SNP = 2pq 2 V P , where p and q are the allele frequencies and β is the estimated allele substitution effect (Falconer & Mackay, 1996). Finally, we determined if SNPs significantly associated with TL were within 100 kb of any gene within the annotated house sparrow genome, because this is the distance that linkage disequilibrium decays to background levels in this species (Elgvin et al., 2017;Hagen et al., 2020). Gene ontology (GO) searches were performed using the gene ontology annotation (GOA) database (Binns et al., 2009;Huntley et al., 2015) to obtain an overview of biological processes and molecular functions known to be influenced by the genes.

| The correlation between tarsus length and telomere length
The model explaining variation in TL that included tarsus length was ranked higher than the model without tarsus length (∆AICc = 2.5, w 1 = 0.78, ER 1 = w 1 /w 2 = 3.55). There was a negative association between tarsus length and TL (β tarsus length = -0.004 ± 0.002, CI = [-0.007, -0.000], n = 2462, Figure 1 and Table 1), such that larger nestlings generally had slightly shorter early-life telomeres. Thus, an increase in (age-corrected) tarsus length of 1 mm was associated with a decrease in TL of 0.8%.

| Parental age effects on offspring telomere length
There was no evidence for associations between offspring TL and   Table S2.2 in Appendix S2). We found no evidence of differences in sex-specific heritabilities of TL (Table S2.3 in Appendix S2).

| Heritabilities and genetic correlations of telomere length, tarsus length, and body condition
In the analyses of the artificially selected populations (Leka and  Table S2.4 in Appendix S2).

| GWA analyses
When controlling for the phenotypic effect of tarsus length on TL, nine SNPs showed evidence for an association with early-life TL (Table 3, Figure 3), with a Bonferroni (1935) corrected threshold of p ≤ 2.77 × 10 −7 at the genome-wide p-value threshold (i.e., the nominal p = .05 divided by 180,650 markers) and a genomic inflation F I G U R E 1 The negative association between age-corrected tarsus length and telomere length (log 10 -transformed) in 2462 house sparrow nestlings with a regression line from a LMM shown in Table 1. The 95% confidence interval (grey) reflects only the fixed effects factor of λ = 1.0489 ± 0.0002 ( Figure S2.3 in Appendix S2). Using the annotated house sparrow genome, a total of 22 genes on five chromosomes were found to be located within proximity (± 100 kb) of six of the top SNPs (Table 4). Four SNPs that showed weak evidence for an association with TL (nominal .05 < p < .10) are also shown in Table 3. Among three of these SNPs we identified three genes within ±100 kb on three chromosomes (Table S2.5 in Appendix S2).
SNPa429690 is located on chromosome 2 within the Aquaporin-1 (AQP1) gene, which encodes the AQP1 water channel membrane protein. The AQP1 protein is abundant in erythrocytes (where TL is measured) and important in regulating body water transport and balance (Nielsen et al., 2002), but also in a range of other physiological functions including cell migration, wound healing, fat metabolism and oxidative stress (Saadoun et al., 2005;Verkman et al., 2014). The same SNP is located 39 kb from the growth hormone-releasing hormone receptor (GHRHR), which controls body growth (Mullis, 2005), and has been associated with telomerase activity (Banks et al., 2010), lifespan (Soerensen et al., 2012) and the progression of several types of cancer (Chu et al., 2016;Schally et al., 2018;Villanova et al., 2019).

Humans with overexpression of growth hormones and consequently
insulin-like growth factor 1 (IGF-1) have shorter telomeres (Aulinas et al., 2013;Deelen et al., 2013;Matsumoto et al., 2015;Monaghan & Ozanne, 2018). SNPa17235 was close (11 kb) to FRMD4B (FERM domain-containing protein 4B), which is involved in epithelial cell polarity that is important in tissue morphogenesis (Ikenouchi & Umeda, 2010). This SNP was also near other genes related to cell proliferation (UBA3 and TMF1), skeletal muscle organization (LMOD3) and oxidative stress (ARL6IP5, see Table 4). SNPa450086 was 76 kb from OXR1 (oxidation resistance protein 1) that regulates expression of several antioxidant enzymes (Volkert et al., 2000). SNPa108592 was in the vicinity (43-84 kb) of several genes on chromosome 15 linked to cell proliferation, ubiquitination and immune response (Table 4). SNPi16410 was closest to SHCBP1 (70 kb) and CDCA4 (76 kb), which are both involved in cell proliferation and probably apoptosis (Asano et al., 2014;Wang et al., 2008;Xu et al., 2018;Zou et al., 2019). SHCBP1 is upregulated by growth factor stimulation (Schmandt et al., 1999). CDCA4 is likely involved in the regulation of hematopoietic stem cells from where erythrocytes (reflecting TL) are derived (Abdullah et al., 2001). Expression of the SCN4A gene (68 kb from SNPa491204) has previously been correlated with TL in human stem cells (Wang et al., 2017). SNPa491204 was 49 kb from the growth hormone gene GH (which is linked to TL as mentioned above, see also Pauliny et al., 2015) and WNT9B (40 kb) of the Wnt/β-catenin signalling pathway, which is modulated by telomerase (Park et al., 2009). In Appendix S2 we mention interesting genes found beyond the ±100 kb limits of the top SNPs.
When not controlling for the effect of tarsus length on TL, the same top SNPs were identified as in the analysis above where tarsus length was included (Table S2.6 in Appendix S2). In addition, SNPa208275 was associated with TL and found 47 kb from FGFR2 encoding a tyrosine-protein kinase that is a receptor for fibroblast growth factors that regulates several aspects of cell proliferation and bone morphogenesis (Table S2.7 in Appendix S2, Katoh, 2009).

| DISCUSS ION
The evolutionary response to selection on TL depends on the additive genetic variance of TL and the strength and sign of any genetic correlations with other traits under selection (Lande & Arnold, 1983). Dugdale and Richardson (2018) criticized past quantitative genetic studies of TL on the main grounds that (1) they applied basic regression analyses that did not consider environmental effects impacting TL and as a consequence of that, additive genetic effects may have been overestimated in previous studies; (2) TL changes and (3) collecting TL data from more than 3300 individuals across four populations, which represent a considerably larger sample size than those of previous wild animal studies.
We found that around 4% of the variation in early-life TL in house sparrows at the end of the nestling growth period was determined by additive genetic variation. The relatively small additive genetic variance and large year variance in early-life TL appears to be in accordance with the effects of relative growth and weather conditions on TL in similar sparrow populations .
The lack of repeated individual TL sampling in this study may prevent us from fully separating between permanent environmental effects and the common environmental effects (brood effects and parental effects, Wilson et al., 2010). However, with several offspring measures for each brood, mother, and father, most of any permanent environmental variance would be included in the residual variance (Kruuk & Hadfield, 2007). In addition, recent longitudinal studies have found negligible permanent environmental effects on TL van Lieshout et al., 2021;Seeker et al., 2018;Sparks et al., 2021;Vedder et al., 2021).
Similarly small but significant heritabilities of TL have been reported using animal models for example, nestling collared flycatch- in which TL correlates with several weather variables. These studies also documented considerable year or cohort effects on TL (Foley et al., 2020;Sparks et al., 2021) similar to studies finding no heritability of TL in white-throated dippers (Becker et al., 2015) and European badgers (Meles meles, van Lieshout et al., 2021). In TA B L E 3 Single nucleotide polymorphisms (SNPs) with evidence (italics) or weak evidence for an association with early-life telomere length in house sparrows (n = 373

F I G U R E 3
Manhattan plot showing genomic location plotted against -log 10 (p-value) of the GWA analysis results for early-life telomere length in house sparrows (n = 373). The dotted line indicates the genome-wide significance threshold (corresponding to p < .05 divided by the number of tests n = 180,650 SNPs) used to determine the top SNPs listed in Table 3 TA B L E 4 Genes found within ±100 kb of SNPs in Table 3 with evidence for an association with early-life telomere length house sparrows Chr.  Gong and Yeh (1999) and Osaka et al. (1998) 12 TMF1: TATA element modulatory factor (Homo sapiens)

67,507
Cell growth, immune response, androgen receptor coactivator Garcia et al. (1992) and Perry et al. (2004) 12  et al., 2011). The heritability of TL in house sparrows is comparable to that of many life-history traits and considerably lower than many morphological traits (e.g., Mousseau & Roff, 1987;Visscher et al., 2008), which may suggest that TL is under strong selection in the wild (Voillemot et al., 2012)  A considerable proportion of the phenotypic variance in TL could be attributed to brood and parental effects ( Figure 2). However, we did not find evidence that parental effects were transmitted through a parental age at conception effect ( Figure S2.1 in Appendix S2).
Paternal age effects, which has been observed in several other species (Eisenberg & Kuzawa, 2018), may not manifest in these house sparrows because the mean age at reproduction in this study was low (around 2 years). Parent-offspring regressions ( Figure S2.2 in Appendix S2) suggested a stronger component of maternal heritability (h 2 maternal = 0.44) rather than paternal heritability of TL. Maternal heritability estimates from parent-offspring regressions includes both direct additive genetic, maternal genetic and maternal environmental effects (Wilson et al., 2005), and we found a lower maternal heritability (h 2 maternal = 0.078) when using the animal models (Table  S2.2 in Appendix S2). Maternal inheritance of TL has been found in several bird species (Asghar et al., 2015;Becker et al., 2015;Horn et al., 2011;Reichert et al., 2015) and in some studies on humans , where this has been attributed to an X-linked gene (Nawrot et al., 2004) or implied genomic imprinting . In our study, we did not detect sex differences in TL heritability (Table S2.3 in Appendix S2). However, we would expect higher heritability for the sex in which TL is less strongly associated with fitness given similar genetic architectures (Merilä & Sheldon, 1999;Roff, 2012). Such an association with fitness was found by Heidinger et al. (2021), where early-life TL was positively related to lifetime reproductive success in house sparrows, but only in females.
Maternal effects on offspring TL are expected to be strongest in Chr.  Kassmann et al. (2008), Tsujino et al. (2003). See also Wang et al. (2017) Bergstein et al. (1997) and Bourhis et al. (2010). See also Park et al. (2009) Note: Chromosome number, distance (in bp) between SNP and gene, general molecular or biological function or relevance to telomere biology are indicated with references. The list is sorted first by SNP p-value and then by gene distance.

TA B L E 4 (Continued)
early-life (Wolf et al., 1998) and could act through for example, yolkdeposited components in the egg Noguera et al., 2020;Stier, Hsu, et al., 2020) or post-laying through maternal care behaviour (e.g., incubation and feeding rate, Viblanc et al., 2020). Our results suggest that such effects may have a genetic basis that will respond to selection: For heritable traits like TL, maternal inheritance of offspring TL may be expected to increase the expected rate of adaptive evolution of TL above what would be expected from the heritability alone (Hadfield, 2012;Lande & Kirkpatrick, 1990;Räsänen & Kruuk, 2007;Wolf et al., 1998).
There was evidence for additive genetic variance in the tarsus length of sparrow nestlings, but the heritability estimate (h 2 = 0.080, Table 2) was considerably smaller than those of adult house sparrows in a larger sample of populations in the same area (Araya-Ajoy et al., 2019;Jensen et al., 2008) and other bird species . However, there was a large brood effect on nestling tarsus length suggesting common environmental effects within broods (e.g., Potti & Merino, 1994). For instance, variation in clutch size, seasonal differences in food availability, weather conditions (Ringsby et al., 2002), and provisioning rates by parents (Ringsby et al., 2009) may induce intraclutch competition and variation in the degree to which nestlings are able to achieve their adult tarsus lengths at fledging (Metcalfe & Monaghan, 2001;Naef-Daenzer & Keller, 1999). Furthermore, measurement error is probably higher for the incompletely ossified nestling tarsi, which are covered by a soft fleshy skin tissue that contributes to the measured length.
Individuals with shorter tarsi (a proxy for structural size, Araya-Ajoy et al., 2019) were found to have longer telomeres, although the effect of tarsus length on TL was small and there was considerable variation in TL for a given size (Figure 1). This confirms previous observations of a prevailing negative phenotypic correlation between body size and TL within house sparrows Ringsby et al., 2015) and other species (Monaghan & Ozanne, 2018). We did not find any evidence of a significant negative genetic correlation between TL and tarsus length (Table 2). Instead, the negative phenotypic association between TL and tarsus length may have no genetic basis but is shaped by common environmental effects that affect both traits in opposite directions (e.g., Kruuk et al., 2008) including processes related to the incomplete replication of chromosome ends during cell division and increased oxidative stress during growth (e.g., Monaghan & Ozanne, 2018). The lack of a genetic correlation between TL, tarsus length or body condition could also be attributed to selection acting simultaneously on some correlated, unmeasured trait . Both with and without controlling for the effect of tarsus length on TL, our GWAS on TL identified several genes involved in skeletal development, cellular growth and differentiation that may regulate body growth or size (Table 4, Tables S2.5 and S2.7 in Appendix S2), which could, however, suggest some genetic basis of the negative correlation between TL and size. For instance, several growth factors were downregulated in telomerase deficient mouse bone marrow stromal stem cells (Saeed & Iqtedar, 2015) suggesting that short telomeres or telomere loss could also be a constraint on proliferation potential. Thus, because several of the genes that may regulate TL during early development appear to also be involved in cell proliferation or morphogenesis, such genes may have co-evolved.
None of the genes highlighted in our analysis have previously been linked to TL in GWA studies (reviewed in the introduction). Table 4 does not provide an exhaustive list of potential biological processes or molecular functions associated with variation in TL, and with little a priori information on the identified SNPs, we are cautious in interpreting these apparent associations as causal (Pavlidis et al., 2012). Furthermore, the low heritability and polygenic nature of TL make it challenging to identify putative causal genes, which consequently only explain a small part of the total phenotypic variance in TL, and our limited sample size (n = 383) is likely to upwardly bias effect sizes and SNP heritabilities due to the Beavis effect (Slate, 2013). Our GWAS on TL was limited to a subset of recruiting individuals, which may affect power to detect associations between SNPs and TL if the genotype or phenotype affects recruitment probability.  found no association between TL and first-year survival in house sparrows, but that recruits had longer tarsi.
Several of the identified candidate genes (Table 4, Tables S2.5 and   S2.7 in Appendix S2) are involved in cell proliferation and apoptosis during which TL and telomerase activity invariably play an important role (Greider, 1998;Masutomi et al., 2003). Telomerase activity has not been investigated in house sparrows. However, for example the RHOF gene product functions cooperatively with CDC42 and Rac to organize the actin cytoskeleton (Ellis & Mellor, 2000), and the latter complex participates in the control of telomerase activity in human cancer cells (Yeh et al., 2005). CDC42 is activated by FGD4 (Chen et al., 2004), which was found within a major locus affecting TL in humans (Vasa-Nicotera et al., 2005). SNPa108592 was found near several genes involved in cell proliferation, differentiation, immune response, and ubiquitination (Table 4). Ubiquitination regulates several shelterin components and telomerase activity (Peuscher & Jacobs, 2012;Yalçin et al., 2017). The closest gene, ORAI1 (43 kb), the keeper of the gates of calcium ions (Homer, 1924), is crucial for lymphocyte activation and immune response (Feske et al., 2006).
Although not linked to ORAI1 mutations, calcium ion levels can modulate telomerase activity (reviewed in Farfariello et al., 2015).
We identified a particularly interesting gene associated with TL, AQP1. The AQP1 channel not only conducts water across cell membranes, but also hydrogen peroxide, a major reactive oxygen species Tamma et al., 2018), and nitric oxide (Herrera et al., 2006), which is an important regulator of oxidative stress (Pierini & Bryan, 2015) and a weak oxidant itself (Radi, 2018). Furthermore, increased availability of nitric oxide may activate telomerase and thereby prevent replicative senescence (in endothelial cells, Vasa et al., 2000). Enhanced oxidative stress associated with endothelial cell senescence may also be mediated by AQP1-regulated nitric oxide flow (Chen et al., 2020;Tamma et al., 2018). In AQP1 knocked-out erythrocytes (where TL was measured) cell lifespan was shortened (Mathai et al., 1996) and angiogenesis is inhibited in AQP1 knockedout chicken embryos (Camerino et al., 2006) and mice (Saadoun et al., 2005). Telomeres are particularly sensitive to ROS and shorten due to oxidative stress during growth (Reichert & Stier, 2017;von Zglinicki, 2002). For example, Kim et al. (2011) found a negative genetic correlation between growth and resistance to oxidative stress in yellow-legged gull (Larus michahellis) chicks, which could be mediated by TL (see also Smith et al., 2016). Another candidate gene, OXR1, 76 kb from SNPa450086, has a well-described antioxidant function (Oliver et al., 2011;Volkert et al., 2000) and is upregulated in senescent human cells . Knockdown of OXR1 increases ROS production and ultimately induces apoptosis (Oliver et al., 2011;Zhang et al., 2018), which could be due to telomere crisis.
Over-expression of AQP1 has been associated with several types of cancer (Verkman et al., 2008), suppression of apoptosis (Yamazato et al., 2018) and may play an important role in tumour biology (Saadoun et al., 2005;Tomita et al., 2017). Other candidate genes including GHRHR, SHCBP1 (Tao et al., 2013), GH (Boguszewski & Boguszewski, 2019), and OXR1 (Yang et al., 2015) are also involved in tumorigenesis. Cancer prevalence is not well-studied in wildlife (Pesavento et al., 2018), but tumours have been documented in house sparrows (Møller et al., 2017). Genes affecting both TL and cancer risk (Jones et al., 2012;Tacutu et al., 2011) could underlie the antagonistic pleiotropy of trade-offs between long telomeres in early-life (with potential benefits to growth, reproduction, and other oxidative stress inducing processes) and later-life cancer mortality (Tian et al., 2018). Cancer is often viewed as a senescence-related pathology (Lemaître et al., 2020). However, the absence of cancer in early-life should not lead us to conclude that a somatic and potentially fitness-related cost is not paid to maintain that status (Thomas et al., 2018).
We have shown that TL is a heritable, polygenic trait with considerable environmental variation and a maternal inheritance component in a wild passerine. It is, however, important that future studies attempt to confirm the putative candidate genes identified here as associated with TL in other wild populations. Even though the additive genetic component was small, selection on variation in TL may produce evolutionary change in TL over time in wild populations.
The large component of variation in early-life TL caused by annual environmental stochasticity suggests that this will generate heterogeneity in TL among cohorts. Although we did not find a negative genetic correlation underlying the negative phenotypic correlation between TL and body size, we may hypothesize that selection for larger nestling size, which may enhance survival until recruitment (Ringsby et al., 1998), will be associated with selection for shorter early-life TL due to nongenetic mechanisms, which can ultimately influence lifespan or reproductive success.

ACK N OWLED G EM ENTS
We thank all field workers who collected data for this study and the island inhabitants that made it possible, laboratory technician for helpful comments on an earlier version of this manuscript, and "Forsknings-og undervisningsfondet i Trondheim" for financial support (to MLP). This work was funded by the Research Council of Norway (274930) and through its Centres of Excellence scheme (223257).
The study was carried out in accordance with permits from the Norwegian Animal Research Authority (FOTS id 11904) and the Ringing Centre at Stavanger Museum, Norway.

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