The role of atmospheric humidity in the evolution of endotherms’ thermoregulatory performance remains largely unexplored, despite elevated atmospheric humidity being known to impede evaporative cooling capacity. Using a phylogenetically informed comparative framework, we tested the hypothesis that pronounced hyperthermia tolerance among birds occupying humid lowlands evolved to reduce the impact of humidity-impeded scope for evaporative heat dissipation by comparing heat tolerance limits (HTL; maximum tolerable air temperature), maximum body temperatures (T_bmax) and associated thermoregulatory variables in humid (19.2 g H₂O m^{− 3}) versus dry (1.1 g H₂O m^{− 3}) air among 30 species from three climatically distinct sites (arid, mesic montane and humid lowland). Humidity-associated decreases in evaporative water loss and resting metabolic rate were 27 - 38% and 21 - 27%, respectively, and did not differ significantly between climatic sites. Decreases in heat tolerance limits were significantly larger among arid-zone (mean ± SD = 3.13 ± 1.12 °C) and montane species (2.44 ± 1.0 °C) compared to lowland species (1.23 ± 1.34 °C), with more pronounced hyperthermia among lowland (T_bmax = 46.26 ± 0.48°C) and montane birds (T_bmax = 46.19 ± 0.92°C) compared to arid-zone species (45.23 ± 0.24°C). Our findings reveal a functional link between facultative hyperthermia and humidity-related constraints on evaporative cooling, providing novel insights into how hygric and thermal environments interact to constrain avian performance during hot weather. Moreover, the macrophysiological patterns we report provide further support for the concept of a continuum from thermal specialization to thermal generalization among endotherms, with adaptive variation in body temperature correlated with prevailing climatic conditions.

Materials and Methods

Study areas

We collected data at three climatically distinct areas (hot arid, mesic montane and humid lowland) between latitudes of 27.90° and 30.04° S in South Africa [see Table 1 for general information and climatic data (Smit et al. 2011, Fick and Hijmans 2017)].

Study species

We collected physiological data under either dry (~1 g H₂O m^-3) or humid conditions (~19 g H₂O m^-3). Overall, our analysis includes data from 627 individuals (humid, n = 307; dry, n = 320) representing 30 species, 15 families and three orders – Passeriformes, Piciformes and Coraciiformes. Humidity and dry protocol data were collected during the austral summer between January – February (2021) [arid-zone species (see Appendix S1: Table S3, S4 and S5a and S5b)] and September 2021 – February (2022) [montane species (see Appendix S1: Table S9 and S10) and lowland species (see Appendix S1: Table S6, S7 and S8)]. Where available, data for species responses under dry air conditions at our montane and lowland sites were obtained from Freeman et al. (2022) (see Appendix S1: Table S1 and S2).

Air and body temperature measurements

A temperature-sensitive passive integrated transponder (PIT) tag (Biotherm 13, Biomark, Boise, ID, USA) was injected into the peritoneal cavity of each bird prior to the commencement of experimentation to measure T_b. Data from the PIT tags were acquired using a reader and transceiver system (HPR +, Biomark, Boise ID, USA). During experimentation, T_air within the metabolic chamber was measured using a thermistor probe (TC-100, Sable Systems, Las Vegas, NV, USA) inserted through a small hole in the side of the chamber and sealed by a rubber grommet.

Experimental protocol

We measured T_b, evaporative water loss and resting metabolic rate using both the dry and humid protocols. Measurements typically lasted 2 – 4 h and began with a bird placed in a chamber at T_air = 28 °C and given at least 1 h to habituate. For the dry protocol, T_air setpoints beginning from T_air = 28 °C were initially increased incrementally by 4 °C to T_air = 40 °C and thereafter increased incrementally by 2 °C until birds reached their thermal endpoints, following Freeman et al. (2022).  For the humid protocol, following habituation at 28 °C, T_air was gradually increased to 34 °C, above which setpoints were increased incrementally by 2 °C until birds reached their thermal endpoints. Although the initial T_air setpoints differed between dry (28 °C) and humid (34 °C) protocols, the initial habituation periods at T_air = 28 °C prior to the commencement of data collection were identical and rates of heating similar between protocols. For these reasons, we do not think these differences between protocols had any effect on observed patterns of thermoregulation, particularly at T_air approaching the upper limits of thermoregulation. Transitions between successive T_air setpoints took 10–15 min. At each setpoint T_air, birds were exposed to stable T_air and humidity for a minimum of 15-20 minutes until concentrations of CO₂ and H₂O were stable for at least 5 min. We used the stepped respirometry protocol involving brief (15-20 min) exposure to each T_air setpoint which has been shown to yield patterns of evaporative water loss, metabolic rate and T_b similar to those using a steady-state protocol where birds experience each T_air setpoint for several hours (Short et al. 2022).

Gas exchange measurements

Evaporative water loss and carbon dioxide production ( ) were measured using an open flow-through respirometry system, with our set-up identical to that described by Freeman et al. (2020, 2022) and described in full in Appendix S1: Section S1 of the supplementary material.

During dry air measurements, flow rates were adjusted to minimise water vapour pressure within the metabolic chamber (mean chamber humidity across sites = 1.07 ± 0.84 g H₂O m^{− 3}) and varied between 3 L min^{− 1} and 24 L min^{− 1}. Humidity measurements were made following methods similar to the dry protocol, with modifications to the respirometry setup permitting the manipulation of in-chamber humidity (see Appendix S1: Figure S1 and Appendix S1: Section S1 for detailed description). During our humidity measurements, mean chamber absolute humidity was 19.21 ± 1.20 g H₂O m^{− 3} and varied by < 1.5 g H₂O m^{− 3} between sites.

Data analyses

Sample sizes (n) for dry or humidity treatments were generally n = 10 individuals per species, but lower for a few (Appendix S1: Table S1). All data were analysed in the R 4.0.5 (R Core Team, 2020) environment, using R Studio 1.1.463 (RStudio, Inc.). For each species, respective inflection T_air values above which T_b, evaporative water loss, the ratio of evaporative heat loss to metabolic heat production (EHL/MHP) and metabolic rate increased rapidly were identified using the R package segmented.lme (Muggeo 2016), with individual identity included as a random effect to account for measurements at multiple T_air values per individual to avoid pseudoreplication. Each response variable including T_b, evaporative water loss, and metabolic rate were analysed above and below inflection points separately using linear mixed-effect models in the R package nlme (Pinheiro J, Bates D, DebRoy S, Sarkar D 2015). Slopes for the relationships of thermoregulatory response variables were estimated as functions of T_air.

We used the “dredge” function in the MuMIn package to undertake model selection (Barton 2019). The standardised model used for within-species analysis included T_air (or T_air−T_b), M_b, the T_air: M_b interaction and Bird ID (individual) as a random factor. The model with the highest rank among competing models was selected using the AIC_c(Akaike information criterion values corrected for small sample size) as well as the Akaike weights (Burnham and Anderson 2002). If competing models were within ΔAIC_c< 2, we selected the most parsimonious model. Significance was assessed at α < 0.05 and values are presented as mean ± SD.

Interspecific analyses

We accounted for the effect of phylogenetic non-independence in observed patterns of our thermoregulatory response variables by constructing a maximum-likelihood tree including all study species using Mesquite (Maddison and Maddison 2014). Making use of the Hackett phylogeny as a back-bone (Hackett et al. 2008), we downloaded 100 phylogenies from www.birdtree.org (Jetz et al. 2012). Determining the necessary branch-length transformations was achieved by comparing an Ornstein-Uhlenbeck model (Martins and Hansen 1997) with a Brownian motion model of trait evolution (Grafen 1989) using AIC values. Lower AIC values were attained for the Brownian motion model and it was therefore retained. We used Pagel’s λ (Pagel 1999) to test for phylogenetic signal in the residual error of our PGLS (phylogenetic least square regression models) while simultaneously estimating regression parameters (Revell 2010),and rescaled our models using the estimates of λ. When testing for λ we included mean M_b for species to account for the known allometric scaling of physiological traits such as basal metabolic rate (McNab 2002, McKechnie and Wolf 2004) and heat tolerance limit (van Jaarsveld et al. 2021). We detected a significant phylogenetic signal for heat tolerance limit (λ = 0.50), metabolic rate (λ = 0.95), evaporative water loss (λ = 0.85) and EHL_MHP (λ = 0.47) across climatic study sites, and we therefore present results from our PGLS analysis and phylogenetically informed post hoc tests (PhylANOVA). The results of conventional analysis (i.e., phylogenetic non-independence not controlled for) are available in the supplementary material (Appendix S1: Table S11). Conventional analyses are mostly consistent with our findings following phylogenetic correction.

The R package caper (Orme et al. 2012) along with the “pgls” function was used to conduct phylogenetic regression analyses. To detect patterns and quantify differences in heat tolerance limit between climatic study areas as well as determine which physiological variables were predictors of heat tolerance limit patterns, we constructed a phylogenetically informed linear mixed effects model (Appendix S1: Table S12). The MuMIn package and “dredge” function was again used to detect which model selection procedure best explained observed patterns of heat tolerance limit under humid conditions (Barton 2019). Model selection was conducted using AICc values and weights. In addition to model selection, we ran analyses to detect auto-correlation among predictor variables (Appendix S1: Table S13, Durbin Watson test) and assessed the normality of residuals using a Shapiro-Wilk test. Model 328 (Heat tolerance limit ∼ Climate + MaxEHL/MHP + EvapScope + Tbmax + MaxEHL/MHP:climate; Appendix S1: Table S12) was selected.

The anova.pgls function in the R-package caper (Orme et al. 2012) was applied to our model output to determine the significance of predictor variables and assess whether the response variable differed significantly among study localities (Appendix S1: Table S11). We conducted post hoc multiple comparison tests taking into account phylogenetic relatedness using the PhylANOVA function in the R package phytools (Revell 2012) to obtain pairwise differences in both response and predictor variables between climatic areas. The PhylANOVA function conducts a simulation-based phylogenetic ANOVA and performs all post hoc comparisons of means among groups providing a p-value by phylogenetic simulation (Garland et al. 1993).