Leukocyte profiles are associated with longevity and survival, but not migratory effort: A comparative analysis of shorebirds

2.1 General field procedures

We captured shorebirds during autumn migration (August–September) at two stopover sites in central and northern Poland: (1) Jeziorsko dam reservoir (51°47′N, 18°40′E), which is one of the major inland stopover sites for shorebirds and other waterfowl in Poland; (2) the mouth of the Vistula River on the Baltic coast (54°21′N, 18°57′E), where large numbers of shorebirds concentrate on sandflats during autumn. The two study sites are located nearly 300 km apart, but there were anecdotal observations of several shorebird individuals using both locations during their autumn migration (unpubl. data), indicating that the sites lie within the same flyway corridor. No large geographical or ecological barriers for migrating shorebirds are located in the vicinity of the study sites. Birds were sampled during autumn migratory period for three main reasons: (1) it allowed us to test for a link between H/L ratio and migratory effort; (2) we were able to sample several species of shorebirds at one geographical location, where all individuals were subject to relatively uniform environmental and ecological conditions at the moment of sampling; (3) baseline stress levels during migration are typically characteristic for most stages of annual cycle, except for breeding, when concentrations of stress hormones tend to increase (Romero, 2002).

Fieldwork took place between 2013 and 2016 (2013, 2015 and 2016 at the Jeziorsko reservoir; 2015 at the mouth of the Vistula River). All birds were captured in the selective catching devices for shorebirds, the walk-in (funnel) traps. Foraging (probing) shorebirds are directed along the special fences and through a set of funnel-shaped entrances to the central capture chamber, from which they are usually unable to escape (Figure S1). Traps are made of wire frames and a cover of thick fishing or wire net. Birds inside the trap often seem not to be aware of being caught, they can move freely (within the borders of the capture chamber) and may continue feeding or sleep until people approach (Busse & Meissner, 2015; P. Minias, pers. observ.). Thus, it seems safe to assume that, in contrast to mist netting, using walk-in traps is unlikely to trigger stress response in shorebirds at the moment of capture. Walk-in traps were controlled every 2 hr and birds were transported in special storing containers to the ringing station, where they were ringed and aged by plumage (Prater, Marchant, & Vourinen, 1977). Sex was only determined for the ruff Philomachus pugnax, which shows marked sexual size dimorphism and, thus, males and females can be identified based on non-overlapping ranges of wing length (Karlionova, Pinchuk, Meissner, & Verkuil, 2007). Since we captured mostly young individuals (in the first year of life), only first-year birds were included in the analyses of leukocyte profiles and all the descriptive statistics are for first-year individuals only. In the study we considered only those species for which at least five first-year individuals were captured and sampled for blood (Table 1).

2.2 Leukocyte profiles

In total, blood samples were collected for 413 individuals from 19 shorebird species (Table 1). Two hundred fifty-four individuals from 17 species were sampled at the Jeziorsko reservoir, while 161 individuals from 12 species were sampled at the mouth of Vistula River (Table S1). This roughly represented 5%–15% of all shorebirds captured each year at each site. Average time between capture and blood sampling at the two sites (Jeziorsko and Vistula River) was similar due to the application of uniform field methods. Transport (from walk-in traps) and handling time of sampled birds was always kept below 1 hr. Although leukocyte profiles has been reported to change between 30 and 60 min since capture in the great tit Parus major (Cīrule, Krama, Vrublevska, Rantala, & Krams, 2012), no changes in heterophils and lymphocyte counts have been recorded within the first hour since capture in other passerines (Davis, 2005) and in shorebirds (Buehler et al., 2008; D'Amico et al., 2017). Thus, shorebird leukocyte profiles were unlikely to be affected by handling stress under our sampling protocol. The ulnar vein of each bird was punctured with a disposable needle and a drop of blood was transferred to a slide, where a blood smear of one cell layer was made. Blood smears were air-dried and stored in darkness until analysis. All smears were stained using the May-Grünewald–Giemsa method and scanned at 1,000× magnification under a light microscope. A random sample of 100 leukocytes from each blood smear was counted and classified into five cell types: heterophils, lymphocytes, eosinophils, basophils and monocytes. The H/L ratio was calculated by dividing the number of heterophils by the number of lymphocytes. To reduce variability, all blood smears were assessed by one of the authors (R.W.). Twenty-five randomly chosen smears were assessed twice to estimate repeatability of H/L ratio measurements. As indicated by the intra-class correlation coefficient (ICC), measurement repeatability was high (ICC = .86, p < .001). An outlier analysis using the criteria of >3 SDs identified 10 outliers in the dataset. All outliers (H/L > 4.10) were removed from the data, as they were likely to represent individuals under strong acute stress. After removing outliers, H/L ratio distribution was still right skewed (skewness = 1.66) and, thus, was squared rooted prior to the analyses at individual level. Median values were used as the species-specific measures of H/L ratio.

To test whether leukocyte profiles were species specific, we assessed intraspecific repeatabilities in heterophil and lymphocyte counts, as well as in H/L ratios using two methods. First, we used the MCMCglmm package (Hadfield, 2010) developed for r statistical environment (R Development Core Team 2013) to calculate ICCs by dividing within-species variance by total variance. However, since estimates of repeatability can be unreliable for unbalanced sample sizes across species (Pap et al., 2015), we also estimated ICCs for five randomly chosen individuals per species, which was consistent with our minimal sample size. For this purpose, we conducted 1,000 tests per trait and averaged all ICC values and corresponding confidence intervals to obtain the final estimates of repeatability using the Irr (Gamer, Lemon, Fellows, & Singh, 2012) r package.

2.3 Migratory traits

Migratory effort was assessed with fat load accumulated by each species during migration and with total migration distance. Species with high fat loads tend to migrate using long-distance non-stop jumps, often covering thousands of kilometres between successive stopover sites, whereas low fat loads are characteristic for species that migrate using short-distance hops (Warnock, 2010). Thus, accumulation of large fat reserves is often associated with greater energetic challenge during migration, and possibly with larger migratory stress. Irrespective of the migration strategy adopted by each species (jumps vs. hops), longer total migration distance requires not only higher energy input, but is also associated with more profound physiological and morphological adjustments (Piersma, 1998). Both fat load and total migration distance have been demonstrated to correlate with several morphological and physiological adaptations to migration in shorebirds, such as haemoglobin concentration (Minias, Kaczmarek, Włodarczyk, & Janiszewski, 2013) and wing shape (Minias, Meissner, et al., 2015).

In order to estimate fat load, the size of subcutaneous fat depots was assessed in the furculum and in the axilla according to an eight-point scale developed for shorebirds, which has been demonstrated to linearly reflect an increment in fat mass (Meissner, 2009). Body mass of each captured individual was also recorded with a digital balance to the nearest 1 g. Fat scores were transformed into the units of mass using regression equations of body mass against fat scores and expressed as percent of lean body mass (LBM) in order to account for interspecific differences in size (see details in Appendix S1). LBMs were estimated with the intercepts from the same regression equations (Table S2).

Total migration distance was assessed for each species at the two levels: (1) the focal population and (2) the entire distribution range (see details in Appendix S1). Since both measures were strongly correlated (r = .81, n = 19, p < .001), only the population-level estimates were used in the comparative analyses to avoid multicollinearity of independent variables (henceforth referred to as the total migration distance). As migration distance clearly depends on the breeding latitude, any relationship between total migration distance and H/L ratio could be mediated by harsher climatic conditions experienced at higher latitudes. To control for this possible bias, mean breeding latitude for each population was calculated as an average from the northernmost and the southernmost distribution limit of breeding range and included in the comparative models.

2.4 Life span and survival

Maximum life span was used as a basic measure of longevity, since it is widely acknowledged to best reflect the ageing rate in vertebrates (de Magalhães, Costa, & Church, 2007). Data on maximum life span were collected using the AnAge database (de Magalhães & Costa, 2009) integrated within the web portal of the Human Ageing Genomic Resources (Tacutu et al., 2013). The database features the highest reported values of maximum life span for a wide spectrum of animal species, including nearly 1,200 birds, and has been commonly used for comparative analysis of vertebrate longevity (e.g. González-Lagos, Sol, & Reader, 2010; Williams & Shattuck, 2015).

Longevity records, as compiled in the AnAge database, may provide reliable information on life span only if adjusted for sampling effort, as probability of recording an extremely old individual increases with sample size (Møller, 2007). Sample sizes used to estimate maximum life span were similar (medium) for all shorebirds included in this study, except for the dunlin Calidris alpina and the red knot Calidris canutus, which had smaller sample sizes. However, since these two species showed two highest values among all longevity records extracted for this study (Table S3), we assumed that they were unlikely to be markedly underestimated due to low sample size. Also, all longevity records came from wild shorebird populations and the data quality was classified as acceptable. Taking all these into account, we concluded that the complied values of maximum life span are likely to reliably reflect interspecific variation in shorebird longevity.

Data on annual survival rate were collected from literature for 15 out of 19 shorebird species (Table S3). Survival and mortality are considered context dependent and may vary considerably between populations from the same species (Sol, Székely, Liker, & Lefebvre, 2007). While multiple survival estimates for one of our focal species, the common ringed plover Charadrius hiaticula, showed great variability (Table S3), repeatability of survival rate for all other species with multiple data (n = 7) was moderately high (ICC = .52). Thus, multiple measures of survival rate were averaged for each species prior to analysis.

Comparative analyses of life span and survival may be biased due to interspecific variation in body size, as large species generally have greater longevity and survive better than smaller species (Møller, 2007). To control for allometry, we extracted residuals of a log–log regression of maximum life span against LBM (R² = 0.23, F_1,17 = 5.00, p = .039). These residuals were used as an independent variable in the comparative analyses and are henceforth referred to as “residual life span.” Log–log regression of survival rate against LBM was non-significant (F_1,13 = 1.92, p = .18) and, thus, we decided to use raw survival data in the comparative analyses. There was a significant positive correlation between survival and maximum life span (r = .59, n = 15, p = .020), while correlation between survival and residual life span approached significance (r = .46, n = 15, p = .081).

2.5 Comparative analysis

Because of shared evolutionary history, species do not provide independent points for analysis. Hence, when conducting analyses of trait correlations across species, phylogeny must be taken into account (Felsenstein, 1985; Grafen, 1989). We ran phylogenetically informed analyses at two different levels: species specific and individual. Species-specific analyses (median H/L ratios as the dependent variable) were conducted using phylogenetic generalized least squares (PGLS) regression models, which were fitted using the ape (Paradis, Claude, & Strimmer, 2004) and nlme (Pinheiro, Bates, DebRoy, & Sarkar, 2014) r packages. To incorporate within-species variation in leukocyte profiles and cross-species differences in sampling effort, we also used Bayesian phylogenetic mixed models (BPMM; Hadfield & Nakagawa, 2010), as implemented in MCMCGlmm (Hadfield, 2010) and mulTree (Guillerme & Healy, 2014) packages in r. In all BPMM we used uninformative priors (variance set to 1 and belief parameter set to 0.002) for both fixed and random effects. In each analysis, we ran two chains with 200,000 iterations, thinning value of 50 and a burn-in period of 150,000, resulting in 1,000 samples per model. Gelman–Ruben statistic (Gelman & Rubin, 1992) was used to check for the consistency of estimates, indicating that the two independent chains converged in each analysis (potential scale reduction values <1.1).

In all comparative analyses, a specific model of evolution has to be assumed. Basic model of evolution is consistent with the Brownian motion process (BM), which describes evolutionary dynamics of a trait evolving in a purely neutral way. Thus, BM takes no account in natural selection and may not adequately describe change in adaptive characters (Butler & King, 2004). Specifically, BM may not be an appropriate model when selection persists through time (evolutionary changes on successive branches are correlated) and when different lineages are subject to the same selective regime (Felsenstein, 1985). Also, features of BM are at odds with the notion of stabilizing selection, which assumes that variance in mean phenotypic traits remains bounded, while in BM-based models variance grows with time in an unbounded fashion (Butler & King, 2004). To weaken the assumption of the BM model, the phylogenetic scaling parameter λ can be incorporated as a branch-length transformation in the regression analysis (Pagel, 1999). In general, λ varies between zero and one, where zero indicates phylogenetic independence (all branches collapse to zero), and a value of one indicates that trait evolution corresponds to the BM model (the internal branch length of the phylogeny remains untransformed; Freckleton, Harvey, & Pagel, 2002). Alternatively, the Ornstein–Uhlenbeck (OU) process with multiple evolutionary optima can be used, as it has been suggested to address adaptive hypotheses more appropriately (Hansen, 1997). While the OU model is considered the simplest mathematical approach to model evolutionary process with selection, it also includes BM as a special case, which means that it may provide a variety of distributions that are consistent with phenotypic evolution under both selection and drift (Butler & King, 2004). In our PGLS models, the OU model performed similar to the BM model adjusted for λ, as indicated by the relative fit of the two models (the Akaike’s Information Criterion: ΔAIC < 2) and, thus, we used both evolutionary approaches (henceforth referred to as PGLS_OU and PGLS_BM) in the analyses. BPMM analyses were run, by default, using the BM model of evolution. In all BM-based models, the phylogenetic scaling parameter λ was set to its maximum likelihood estimate. λ for H/L ratio and absolute counts of heterophils and lymphocytes were estimated using phytools (Revell, 2012) r package.

Due to differences in sample size (number of species), PGLS and BPMM were constructed separately for residual life span and annual survival, which were entered as covariates. Migratory fat load, total migration distance, breeding latitude and LBM were included as the remaining covariates in each model. Median and individual H/L ratios were entered as the dependent variables in PGLS and BPMM respectively. All analyses were also repeated for the absolute counts of heterophils and lymphocytes. To obtain more parsimonious reduced models, we removed non-significant (p > .10) predictors from the initial full models. Phylogeny was reconstructed based on the most recent complete avian time-calibrated phylogeny (Jetz, Thomas, Joy, Hartmann, & Mooers, 2012) with a backbone tree developed by Ericson et al. (2006). To account for phylogenetic uncertainty, all models PGLS and BPMM were run for 1,000 and 100 alternative trees respectively. All phylogenies were downloaded from the BirdTree database (http://www.birdtree.org, Jetz et al., 2012). p Values were inferred based on z-score tests.

We also assessed the role of evolutionary history in explaining the current-day variation in H/L ratio of shorebirds. For this purpose, we used two methods that followed Sol et al. (2007): (1) nested ANOVA model was used to identify the taxonomic level associated with major diversification in H/L ratio; (2) spatial autocorrelation statistic Moran's I (Gittleman & Kot, 1990) was used to assess phylogenetic autocorrelation in H/L ratio. Phylogenetic correlograms of normalized Moran's I (I/I_max) were used to assess the strength of autocorrelation H/L ratio at different taxonomic levels (genus, family and order). The analysis of Moran's I was conducted using the ape package in r.