Biologists aim to explain patterns of growth, reproduction, and ageing that characterize life histories, yet we are just beginning to understand the proximate mechanisms that generate this diversity. Existing research in this area has focused on telomeres but has generally overlooked the telomere’s most direct mediator, the shelterin protein complex. Shelterin proteins physically interact with the telomere to shape its shortening and repair. They also regulate metabolism and immune function, suggesting a potential role in life history variation in the wild. However, research on shelterin proteins is uncommon outside of biomolecular work.

Intraspecific analyses can play an important role in resolving these unknowns because they reveal subtle variation in life history within and among populations. Here, we assessed ecogeographic variation in shelterin protein abundance across eight populations of tree swallow (Tachycineta bicolor) with previously documented variation in environmental and life history traits. Using blood gene expression of four shelterin proteins in 12-day old nestlings, we tested the hypothesis that shelterin protein gene expression varies latitudinally and in relation to both telomere length and life history.

Shelterin protein gene expression differed among populations and tracked non-linear variation in latitude: nestlings from mid-latitudes expressed nearly double the shelterin mRNA on average than those at more northern and southern sites. However, telomere length was not significantly related to latitude.

We next assessed whether telomere length and shelterin protein gene expression correlate with 12-day old body mass and wing length, two proxies of nestling growth linked to future fecundity and survival. We found that body mass and wing length correlated more strongly (and significantly) with shelterin protein gene expression than with telomere length.

These results highlight telomere regulatory shelterin proteins as potential mediators of life history variation among populations. Together with existing research linking shelterin proteins and life history variation within populations, these ecogeographic patterns underscore the need for continued integration of ecology, evolution, and telomere biology, which together will advance understanding of the drivers of life history variation in nature. 

Study populations: Data were collected from 8 populations in the eastern United States, spanning nearly 10 degrees of latitude (Table 1, Fig 2A): Ithaca, New York (42.28°N, 76.29°W); Amherst, Massachusetts (42.22°N, 72.31°W); Linesville, Pennsylvania (41.65°N, 80.43°W); Bloomington, Indiana (39.17°N, 86.53°W); Lexington, Kentucky (38.11°N, 84.49°W); Knoxville, Tennessee (35.90°N, 83.96°W); Davidson, North Carolina (35.53°N, 80.88°W); and Santee, South Carolina (33.49°N, 80.36°W). These populations do not represent the entire breeding range of this species and in particular, do not extend to the northern edge in Canada and Alaska. All methods were approved by institutional IACUCs and conducted with appropriate state and federal permits.

Sampling of nestlings: Nest boxes were monitored for hatch dates, but in cases where hatch dates were missed (e.g., due to weather or COVID-related staffing shortages), hatch dates were estimated using existing growth curves (McCarty, 2001; Wolf et al., 2021) and accounted for in all statistical analyses. Data from multiple populations shows that the average peak of postnatal growth occurs around 6-days old (McCarty, 2001; Wolf et al., 2021). Growth then slows and plateaus near adult size by 12-days old, just as feather development accelerates. We targeted 12-day old nestlings because they have just completed the rapid period of postnatal growth. Many studies therefore use morphological data at this critical time period as a proxy of nestling growth (Gebhardt-Henrich & Richner, 1998; Haywood & Perrins, 1992; Magrath, 1991; Martin et al., 2018; McCarty, 2001). Population variation in growth rates occurs primarily after peak growth but does not map neatly onto latitude, at least not in the northern (historical) range where previous research has been focused (Ardia, 2006; McCarty, 2001).

We sampled nestlings at 12.03 ± 0.01-days old (hatch day = day 1, range = 10 – 14 days). We sampled ~30 nests per population (Table 1), though logistical constraints prevented collection of RNA in Kentucky. Upon arrival at each nest, we immediately collected whole blood from the brachial vein of 2-3 nestlings per nest (≤ 200 µl, below the maximum suggested volume based on body mass; Gaunt et al., 1997), and we avoided obvious runts with atypical growth. We collected blood in separate tubes for DNA and RNA analyses. We banded nestlings with a USGS band and weighed them to the nearest 0.1g. We also measured flattened wing length using a wing ruler. We stored blood on ice or dry ice in the field, and later stored it at -80°C.

Due to limited budgets, we made the decision a priori to conduct laboratory analyses for a single nestling per nest. When possible, we selected the nestling with the median mass. If the median-massed nestling was not bled or failed to produce a sufficient blood sample, we selected the nestling with the closest mass to the median. In nests with even brood sizes, we randomly selected one of the two nestlings with median mass for telomere and gene expression analyses. In all states except Indiana, telomere length and gene expression data come from the same individual.

 qPCR for Telomere length: We extracted DNA from whole blood (following Wolf et al., 2022) and used primers telc and telg (adapted from Cawthon, 2009) to quantify telomere length relative to the single copy gene GAPDH. Samples were run in triplicate, and mean values were used to calculate the T/S ratio of telomere repeat copy number (T) to our single gene copy number (S) using the formula: 2^-∆∆Ct, where ∆∆C_t = (C_t ^telomere – C_t ^GAPDH) _reference – (C_t ^telomere – C_t ^GAPDH) _sample. The intraclass correlation coefficient (ICC) for intraplate repeatability was 0.951 ± 0.03 (95% CI = 0.944, 0.957) for GAPDH C_t values and 0.926 ± 0.09 (95% CI = 0.916, 0.935) for telomere C_t values. The ICCs for interplate repeatability were 0.96 ± 0.03 (95% CI = 0.87, 0.98) for GAPDH Ct values, 0.89 ± 0.06 (95% CI = 0.73, 0.95) for telomere Ct values, and 0.79 ± 0.10 (95% CI: 0.54 - 0.90) for the T/S ratio (based on 2^-∆Ct values). Plates (n = 13) were balanced by population, sex, relative date of sampling within each population, and brood size.

Nestling Sexing Protocol: Nestlings were molecularly sexed using DNA following established methods (Griffiths et al., 1998; Wolf et al., 2022).

Shelterin Protein Primer Design: Shelterin proteins are relatively conserved across taxa (de Lange, 2018; Myler et al., 2021) and earlier work has identified at least four shelterin proteins in the chicken (De Rycker et al., 2003; Konrad et al., 1999; Tan et al., 2003; Wei & Price, 2004). Our shelterin protein primer sets were developed using the tree swallow transcriptome (accession #GSE126210; Bentz et al., 2019). TRF2 exhibits multiple variants in passerines, and a BLAST search confirmed that our primer set targets TRF2 in closely related barn swallows (Hirundo rustico). TPP1 and POT1 each have a single transcript in adult tree swallows that is highly expressed across tissues, and BLAST searches confirmed that our primer sets targeted TPP1 and POT1 genes, respectively, in multiple bird species. We also designed primers for RAP1 based on tree swallow transcripts of TRF2IP (TRF2-interacting protein), a common alias for RAP1. However, this study omits TRF1 due to negligible expression in nestling blood, and TIN2 because we could not confidently identify the passerine sequence for TIN2. Thus, altogether we quantified gene expression for four key components of the shelterin complex: TRF2, RAP1, TPP1, and POT1 (primer sequences in Table S1).

qPCR for Shelterin Protein Gene Expression: We extracted RNA using a phenol-chloroform-based Trizol method (Invitrogen, Carlsbad, CA) with PhaseLock tubes (5PRIME, #2302830). We synthesized cDNA using 1µg RNA and Superscript III reverse transcriptase (Invitrogen), treated with DNAase (Promega, Madison, WI) and RNase inhibitor (RNAsin N2111, Promega). For each gene of interest, we used the 2^-∆∆Ct method of quantification (Livak & Schmittgen, 2001), in which expression is normalized to the geometric mean Ct of two reference genes for each sample (Vandesompele et al., 2002), and relative to a calibrator sample on each plate. Reference genes correct for technical variation in cDNA quantity across samples, and as such, must (i) be highly expressed, (ii) exhibit low variability among samples, and (iii) show no significant variation among biological categories of interest. Our reference genes were PPIA (peptidylprolyl isomerase A; Virgin & Rosvall, 2018) and MRPS25 (Mitochondrial Ribosomal Protein S25; Woodruff et al., 2022). Preliminary work showed that New York samples exhibited markedly higher gene expression of these and a third reference gene (GAPDH). This violates assumption (iii) of the 2^-∆∆CT method, and we therefore had to omit New York gene expression data. The remaining six populations exhibited limited among-population variation in reference gene expression (non-significant state differences or ≤ 0.5 Ct of the study-wide average).

Samples were run in triplicate alongside No Template Controls (NTCs), using PerfeCta SYBR Green FastMix with low ROX (Quanta Biosciences, Gaithersburg MD) on 384-well plates using an ABI Quantstudio 5 machine with Quantstudio Design & Analysis software (v1.4.3, Thermo Fisher Scientific, Foster City, CA). Each well included 3µL of cDNA diluted 1:50 (or 3µL water, for NTCs) and primers diluted to 0.3µM in a total volume of 10µL. All reactions use the following thermal profile: 10 min at 95°, followed by 40 cycles of 30 s at 95°, 1 min at 60°, and 30 s at 70°, with a final dissociation phase (1 min at 95°, 30 s at 55°, and 30 s at 95°) that confirmed single-product specificity for all samples. All samples fell within the bounds of the standard curve and reaction efficiencies were within 100 ± 15%. Each gene was run on 1.5 plates, balanced by population. Intraclass correlation coefficients for triplicates were 0.996 ± 0.01 (95% CI = 0.995, 0.997) for PPIA C_t values, 0.994 ± 0.009 (95% CI = 0.993, 0.996) for MRPS25 C_t values, 0.975 ± 0.05 (95% CI = 0.967, 0.982) for POT1 C_t values, 0.940 ± 0.05 (95% CI = 0.923, 0.954) for TRF2 C_t values, 0.975 ± 0.05 (95% CI = 0.968, 0.981) for TRF2IP C_t values, and 0.996 ± 0.007 (95% CI = 0.995, 0.997) for TPP1 C_t values.