Hummingbird blood traits track oxygen availability across space and time

Predictable trait variation across environmental gradients suggests that adaptive evolution repeatedly finds the same solution to a challenge. Trait-environment associations can reflect long-term, genetic evolution across phylogenies, or short-term, plastic responses of individuals. When phylogenetic and population-level patterns match, it implies consistency between the long timescale of adaptation and the short timescale of acclimatization. Alternatively, genetic adaptation can find solutions that ‘break the rules’ of trait-environment covariation. For example, blood-hemoglobin concentration ([Hb]) increases at high elevations in animals, but genetic adaptations in some populations have been shown to augment tissue oxygenation while curtailing hemoglobin production, altering this predictable [Hb]-elevation association. Here, we tested whether species adaptation to elevation generally alters trait-environment relationships for blood. To do this, we measured blood traits of 1,217 individuals representing 77 species of Andean hummingbirds, across a 4,600 m elevational gradient. We used hierarchical Bayesian modeling to estimate blood trait responses to elevation, environmental temperature, precipitation, individual and species characteristics, and phylogeny. Strikingly, the effects of elevation on [Hb] and hematocrit (Hct) were nearly identical for individuals and species, implying that rules of elevational blood variation are set by physics of gas exchange in the hummingbird respiratory system and are unchanged by species adaptation. However, when we looked at mechanisms of [Hb] adjustment—by changes in red blood cell size or number—we did find a signal of species adaptation: To adjust [Hb], species at low and high elevations, respectively, tended to adjust cell size, whereas species at mid-elevations tended to adjust cell number. Despite scale-independent elevational variation in [Hb] and Hct, the species-specific balance of red blood cell size versus number appears to have been affected by adaptations that distinguish hummingbird species living at moderate versus high elevations.


Introduction
We used the individual-level dataset to model within species (i.e., interspecific) variation 157 using multilevel Bayesian models. For each of six blood parameters, we estimated the effects of 158 elevation, body mass, precipitation, air temperature, species-typical Hb β A 13-β A 83 genotype, as 159 well as within-species variation in elevational position (an individual's sampled elevation 160 relative to its species elevational range), body mass, air temperature, and precipitation. We 161 accounted for between-species differences by incorporating a species grouping variable. 162 Separately, we excluded sex and wing-loading as potentially important predictors (see Methods). 163 In total, we compared six models per set using LOOIC and WAIC (Tables S3-S4; Methods). In 164 addition to fixed effects, models included species as a random effect, and a phylogenetically 165 correlated random intercept (i.e., covariance matrix of phylogenetic distances among taxa). Thus, 166 our full model in each set reflects the hypothesis that all predictors, phylogeny, and species 167 identity influence blood characteristics. 168 Based on the top-performing model in each set, we found that for every 1,000-m increase 169 in elevation, [Hb] is expected to increase by 0.39 g/dl, Hct is expected to increase by 1.42 170 percentage points, TRBC is expected to increase by 310,000 cells per µl, and MCH is expected 171 to decrease by 0.48 pg. Elevation strongly predicted [Hb], Hct, TRBC, and MCH (Figs 2 A-F 172 and 3A-F). Elevational position was also a strong predictor of Hct, indicating that individuals 173 located higher in their species' elevational range tended to have higher Hct (Fig 3B). [Hb], 174 TRBC, and MCH also increased with increasing elevational position, although these 175 relationships were subtle (Fig 3A, C, E). 176 The finding that [Hb], Hct, and TRBC increase with elevation was consistent with 177 theoretical expectations for O2-carrying capacity of blood and with prior empirical findings in 178 birds [26,30,31,33,34], humans [41,42], mammals [11,43], and other taxa [44,45]. In prior work, 179 substantial variation in slopes of these parameters across elevation may have been attributable to 180 variation in temporal, taxonomic, and spatial scales of comparison, as well as constraints on the 181 time-course of evolutionary adaptation [22,42,46]. We expect increases in Hct at higher 182 elevations because an increase in the total mass of Hb enhances oxygenation under lower 183 ambient PO2 [18]; unless there is a proportional increase in plasma volume or MCHC, an 184 increase in Hb mass is expected to be accompanied by an increase in the proportion of blood 185 comprised of red blood cells. TRBC also increased with elevation and was, along with MCV, 186 one of two key drivers of [Hb] adjustment (Fig 2, Fig 4; [38]. MCV did not vary predictably in 187 our linear models ( Fig 2D). Although increasing within-cell [Hb], or MCHC, would seem to be 188 an alternative mechanism of increasing O2-carrying capacity, MCHC was less variable than other 189 blood parameters (Table S2), likely due to physical constraints [47,48]; accordingly, MCH and 190 MCV were strongly correlated (Table S5; Fig 4C). Contrary to recent empirical findings [30,31], 191 adjustments to MCHC appeared to be neither substantial nor predictable with respect to elevation 192 within or among species (Fig 2F). 193 We found that MCV and MCH decreased in individuals occurring in warmer 194 environments (Figs 3D, E), and that MCHC decreased in individuals occurring at warmer 195 localities within a species' range ( Fig 3F) (Tables S3-S4; Methods). Moderate effects of the species grouping variable suggest that 204 variation in blood traits is partly explained by aspects of species biology that evolve too fast to 205 have resulted in phylogenetic autocorrelation. Species-typical Hb β A 13-β A 83 genotype did not 206 predict [Hb] or any other blood characteristics (Fig 3A-F), providing no evidence that 207 evolutionary adjustments to O2-binding affinity affect the physical composition of blood or its 208 variation with elevation. 209 210

Among-species blood trait-environment variation 211
Using the species-level dataset, we evaluated the effects of elevation, body mass, species-typical 212 Hb β A 13-β A 83 genotype, and ambient temperature and precipitation on each of our 6 blood 213 response variables. We compared sets of models with LOOIC and WAIC (Table S3) Within-and among-species findings for elevational effects were similar (Fig 3). Based on 217 the top model in each set, for every 1,000-m increase in elevation, [Hb] is expected to increase 218 by 0.29 g/dl, and Hct is expected to increase by 1.55 percentage points. Elevation was a 219 consistently strong predictor of [Hb], Hct, and TRBC, all of which increased with increasing 220 elevation (Figs 2G-I and 3A-C). The importance of the quadratic effect of elevation for MCV 221 and MCH, and overall shape of the species-level data, suggest that O2 carrying capacity 222 adjustment differs for middle-elevation species versus those found at low and high elevations 223 ( Fig 2J-K; Fig 3D-E). Neither species-typical Hb β A 13-β A 83 genotype, body mass, 12 environmental temperature, nor precipitation explained variation in any species-mean blood 225 characteristic (Fig 3). 226

227
Mechanisms of blood-O2 carrying capacity adjustment 228 The 'banana' shape of the relationship between cell number versus size indicates that the highest 229 levels of blood-O2 carrying capacity per unit blood are achieved at high values of TRBC or 230 MCV, respectively. Our data support that the axis of variation in cell number versus size space is 231 tightly constrained (Fig 4A-B  Extending our analysis, we designed a set of models to test potential predictors of CeNS: 257 elevation, a quadratic elevation effect (elevation 2 ; to account for notable curvature in the data), 258 body mass, wing loading, species-typical Hb β A 13-β A 83 genotype, and phylogeny (see 259 Methods). Elevation, elevation 2 , and body mass strongly predicted CeNS (Fig 4D-F; Table S7). 260 We found a trend of decreasing CeNS with higher wing loading, but the 95% credible intervals 261 overlapped zero (Fig S2A-B). Species-typical Hb β A 13-β A 83 genotype did not predict CeNS 262 ( Fig S3), nor did we find phylogenetic signal in CeNS (Fig S4). 263 Examining the tradeoff between cell number and volume provides a potential explanation 264 for the striking pattern that middle elevation species tend to adjust [Hb] by different mechanisms 265 than low and high elevation species (Fig 4E). CeNS was uncorrelated with TRBC coefficients 266 (R 2 = 0.07, r = -0.27) but strongly correlated with MCV coefficients (R 2 = 0.51, r = -0.71), 267 indicating that variation in CeNS was primarily driven by the degree to which species modulate 268 cell size (Fig S2, Fig S5. We suggest that because low and high elevation species tend to possess 269 14 relatively small cells, they may be more able to rely on shifts in cell size, rather than cell number, 270 to increase [Hb] (Fig 2J; Fig 4). However, as cell size increased, species were more likely to 271 approach maximum size for red blood cells, beyond which they could incur a functional 272 detriment. As a result, species with larger red blood cells were more reliant on cell number to 273 adjust [Hb]. For middle elevation hummingbird species that have evolved under moderate (~15-274 30%) reductions in O2 availability, TRBC appears to be the primary axis for adjustment of 275 blood-O2 carrying capacity per unit volume of blood. Species that have evolved at the highest 276 elevations are expected to be the best adapted to hypobaric hypoxia, and it is therefore 277 remarkable that their mode of [Hb] adjustment (primarily by MCV) resembles that of 278 hummingbird species in the lowlands, the putative ancestral elevation for hummingbirds [37]. 279 280

Hummingbirds at extreme elevations 281
Do hummingbirds at extreme high elevations follow the same rules as hummingbirds 282 distributed at low and moderate elevations with respect to elevational blood variation? Our 283 models encompassed the entire elevational gradient and therefore could have missed important 284 phenomena such as threshold effects. To address this, we compared individual hummingbirds 285 sampled above and below 4,200 m elevation, a threshold that marks the upper ~10% of the 286 gradient. Hummingbirds sampled from these extreme high elevations had significantly smaller 287 cells and higher erythrocyte counts than others ( Fig S6). This finding suggests that high-altitude 288 specialist hummingbirds tend to have smaller cells, consistent with findings in other high altitude 289 vertebrates [11,29,44,52,53]; however, this pattern was only detectable under the extreme 290 conditions at the uppermost limits of the Andean elevational gradient. We proceeded to formally 291 test for elevational breakpoints using segmented regression (see Methods). Evidence for 292 objective blood trait thresholds along the elevational gradient was mixed, with breakpoint 293 estimates for individual traits ranging from ~1,400 to ~4,000 m ( Fig S7). Lack of consilience in 294 these estimates suggests that trait-specific thresholds are challenging to identify and may not be 295 generalizable across species and clades. 296 297

Phylogenetic signal in blood traits 298
Phylogenetic signal was low (λ = 0.0; K range = 0.14-0.21) across all blood traits for both 299 within-and among-species models. Blomberg's K and Pagel's λ results were consistent with 300 visual assessment of continuous trait maps, even after removing species with <3 individuals 301 sampled (Table S8; Fig S8). We additionally found no phylogenetic signal in CeNS ( Fig S4). 302 Non-existent to low phylogenetic signal in blood traits within and among hummingbird 303 species was surprising given demonstrated, phylogenetically conserved Hb genetic adaptations to 304 elevation [21,22] and in light of strong phylogenetic signal reported from recent comparative 305 analyses of blood traits that had broader taxonomic sampling [33,34]. However, it is likely that 306 diverse avian clades possess different blood trait optima; in previous comparative analyses that 307 included multiple families or orders, phylogenetic signal may have been driven primarily by 308 differences among clades. For example, Minias et al. [33] reported that small passerines had 309 consistently high [Hb] and Hct estimates, and that including the shorebird order, 310 Charadriiformes, which was characterized by high [Hb] but medium Hct values, strongly 311 influenced patterns of phylogenetic conservatism [33]. We tested for and found no interaction of 312 clade and elevation on the six studied blood traits, with the exception of the interaction of 313 elevation and the Hermit clade in the [Hb] model set ( Fig S9). Interestingly, hermits are a large, 314 early diverging hummingbird clade that contains diminishing species diversity with altitude (and 315 16 no species that exceeds 3,000 m), and its apparent increased sensitivity of [Hb] to elevation hints 316 at a functional explanation for this biogeographic tendency. The existence of clade-specific trait 317 optima [33] and strong species effects, demonstrated in this study, may at first seem difficult to 318 reconcile with our findings that there was no within-family phylogenetic signal in 319 hummingbirds; however, we view these findings as being consistent with steady evolution of 320 trait optima by a mean-reverting, Ornstein-Uhlenbeck process, as described by Minias et al. [33], 321 and marked by occasional major jumps that coincide with changes in body plan or respiratory 322 demands. 323 A strength of our study was that hummingbirds vary only subtly in aspects of their 324 biology that could plausibly confound our analyses, such as metabolic rate, energy budget, torpor 325 use, mode of locomotion, foraging strategy, or clutch size. Phylogenetic conservatism in these 326 idiosyncrasies--not phylogenetic affinities--were the primary drivers of variation in blood traits 336 affecting O2-carrying functions. In response to elevational changes in O2 partial pressure, 337 hummingbirds make predictable adjustments to the physical composition of blood. Remarkably, 338 after controlling for many potentially confounding variables, decreasing PO2 drove increases in 339 [Hb] and Hct that were nearly identical in slope within and among hummingbird species, 340 respectively. This implies that developmental and evolutionary optimization of blood traits to 341 ambient PO2 occurs in fundamentally similar ways, responding to the same set of physical rules 342 that govern gas exchange at the blood-gas barrier. In addition to O2 availability, species identity 343 itself explained 13-21% of blood trait variation, suggesting that fast-evolving, unmodeled 344 aspects of species biology underpin substantial hematological variation. Our data revealed a 345 striking constraint space that defines the cell number-size tradeoff for hummingbird erythrocytes; 346 within this space, individuals and species co-equally modulated cell number and size to adjust 347 collection, we pipetted a 10 µl sub-sample of blood-NaCl dilution into a hemocytometer and 382 allowed cells to settle for 1 min. After focusing the microscope, one or more photos were 383 19 obtained of cells in the hemocytometer grid. Later, cells were manually counted in the ImageJ 384 software [58] and TRBC was estimated following standard protocols [28,59]. 385 After drawing blood, birds were humanely killed by rapid cardiac compression. We 386 recorded collection elevation, latitude, longitude, sex, body mass, and date. Body mass was 387 typically recorded prior to blood sampling; however, body mass was corrected for estimated 388 volume of blood drawn for any individuals massed after blood sampling. Each collected bird was 389  We used principal components analysis (PCA) for temperature (Bio1-11) and precipitation 420 (Bio12-19) variables to create a composite measure of temperature and precipitation across the 421 gradient. PC1 of the temperature variables (hereafter "temperature") explained 73.4% of the 422 variation in composite temperature, and loadings suggested that this axis corresponded to 423 increasing temperature across sites (see GitHub). PC1 of the precipitation variables (hereafter 424 "precipitation") explained 79.8% of the variation in composite precipitation, and loadings 425 suggested that this axis corresponded to increasing precipitation across sites (see GitHub). 426 To account for multiple measurement effects (i.e., traits measured from multiple 427 individuals within species), we calculated between-and within-species differences following 428 Villemereuil et al. [65]. This approach uses within-group centering to separate each predictor 429 21 into two components: We first calculated species mean values (i.e., between-species variability) 430 and then subtracted the mean value of individual observations (i.e., within-species differences). 431 These measures of within-species variation for temperature, precipitation, and body mass were 432 included in models. Because hummingbird body masses may differ by sex, within-species 433 differences were calculated using mean body mass of males for males, mean body mass of 434 females for females, and species-level means for individuals of unknown sex (n=6). 435 We used wing area data from Skandalis  We initially constructed phylogenetic models (results not shown), and calculated phylogenetic 453 signal using the full and reduced models for each blood parameter, following the vignette and 454 recommendations of P. Bürkner (https://cran.r-project.org/web/packages/brms/ 455 vignettes/brms_phylogenetics.html), using the 'hypothesis' method. Estimates of phylogenetic 456 signal for all model sets (full and reduced in both within-and among species analyses) were low 457 (λ < 0.004); these estimates are supported by visualizations of continuous trait maps of mean 458 blood parameters (SI Appendix; Fig S8). Accordingly, we excluded phylogeny as an effect in 459 models and evaluated phylogenetic signal using Blomberg's K and Pagel's λ. Results from the 460 two tests were concordant and consistent across model subsets (Table S8). predictors whose 95% credible intervals (CIs) did not overlap zero in full models. Distribution 478 family differed by model set: Gaussian for [Hb], Hct, and TRBC models; skew-normal for MCV 479 and MCH models; and student's t-distribution for MCHC models. We ran four Markov chains 480 for 10,000 iterations with a burn-in period of 5,000, thinned every 10 steps. 481 We tested for interactions between clade and elevation by comparing three models: 1) the 482 top within-species model for each blood parameter; 2) the top model + (1|Clade); and 3) the top 483 model with an elev:Clade interaction. We found no meaningful interaction of elevation and clade 484 except for elev:CladeHermit in the [Hb] model, whose 95% credible intervals did not overlap 485 zero (Fig S9). In accordance with lack of phylogenetic signal in our data, the effect of clade Species mean values were calculated for all species with n=2 or more individuals sampled. We 493 split Patagona gigas into "north" and "south" haplotypes in light of evidence that the complex 494 may be comprised of two functionally and genetically distinct populations separated by strong 495 migratory divides [77], and we analyzed data only from individuals sampled during the breeding 496 months to ensure correct population assignment. We designed model sets for each of our six 497 blood traits. For each response, we analyzed distributions (qqPlots, Cook's Distance, leverage) 498 24 and removed outliers that strongly influenced fit. For each model, predictor variables included 499 species-typical Hb β A 13-β A 83 genotype, mean elevation, mean body mass (log-transformed; one 500 mean for all individuals due to balanced sampling between sexes), mean temperature, and mean 501 precipitation. All variables were standardized to a mean of zero and standard deviation of one. 502 For each response, we compared: 1) an intercept-only null model, 2) a model with all 503 predictors, and 3) a reduced-predictor model with only those predictors whose 95% CIs did not 504 overlap zero in full models (reduced model 2). Given the notable quadratic shape of species 505 mean data for TRBC, MCV, MCH, and MCHC (Fig 2G-L We then analyzed the factors that contribute to variation in CeNS. Our model response 534 variable was CeNS (0-1 index). Predictors, all species mean values, included elevation, a 535 quadratic elevation term "elevation 2 " to account for notable curvature in the data, body mass 536 (log-transformed), wing loading, and species-typical Hb β A 13-β A 83 genotype. We compared six 537 models: 1) an intercept-only null model, 2) a model with core predictors (elevation, body mass, 538 and wing loading), 3) a model with core predictors plus elevation 2 , 4) a model with core 539 predictors plus elevation 2 and species-typical Hb β A 13-β A 83 genotype, 5) a reduced version of 540 model 3, and 6) a reduced version of model 4. Reduced models included only those predictors 541 whose 95% CIs did not overlap zero in full models. In addition to model fit diagnostics, we 542 assessed correlations between variables by analyzing correlation matrices and by regressing 543 CeNS on each R 2 values from standardized species-specific CeNS models, unstandardized TRBC 544 and MCV beta estimates, and raw TRBC and MCV values (Fig S5). Distribution family differed 545 by model set: Gaussian for individual-level, species mean, and CeNS prediction models, and 546 student's t-distribution for species-specific models. For all CeNS models we ran four Markov 547 chains for 20,000 iterations with a burn-in period of 10,000, thinned every 10 steps. We assessed 548 phylogenetic signal with Blomberg's K, Pagel's λ, and visually with contmaps ( Fig S4). We tested for objective thresholds in each blood trait using piecewise linear regression of 567 elevation (standardized) regressed on each blood trait. In this procedure, we first fit a simple 568 linear regression model, iteratively searched for objective "breakpoints" along our sampled 569 elevational gradient, selected the breakpoint with the lowest residual error, and refit the linear 570 regression model with separate estimates for slopes and intercepts below and above the 571 breakpoint, respectively. We constrained the model to estimate breakpoints that were each 400 m 572 from the lowest and highest points in our data (i.e., from 400 m to 4,178 m in elevation), to avoid 573 biasing breakpoints towards small sample sizes and each low and high elevations, and anomalous 574 species at both ends of the distribution. We evaluated goodness of fit using fit metrics and by 575 visually assessing QQplots, Cook's D, leverage, and residual distributions. We compared models 576 using AIC values. All breakpoint models fit significantly better than simple linear regression 577 models except for the TRBC breakpoint model; this model did fit better in all respects than the 578 simple linear regression model, but the fit was not significant at the p=0.05 level. We used a one-579 way ANOVA to compare whether individuals below and above model-estimated breakpoints 580 differed in their blood characteristics (Fig S7). value from a single species. In all panels, solid trend lines correspond to relationships where 95% 854 posterior credible intervals in models did not overlap zero, while dashed trend lines correspond 855 to relationships where 89% credible intervals in models did not overlap zero (see Fig 3). 856 Credible intervals for both elevation (at the 95% level) and elevation 2 (at the 89% level) did not