DOI: 10.5061/dryad.cc2fqz66z
title: Playa-Montmany, et al.: The thermoregulatory role of relative bill and leg surface areas in a Mediterranean population of Great tit

description: There is growing evidence on the role of legs and bill as ‘thermal
 windows’ in birds coping with heat stress. However, there is a lack of
 empirical work examining the relationship between the relative bill and/or
 leg surface areas and key thermoregulatory traits such as the limits of
 the thermoneutral zone (TNZ) or the cooling efficiency at high
 temperatures. Here, we explored this relationship in a Mediterranean
 population of Great tit (Parus major) facing increasing thermal stress in
 its environment. The lower and upper critical limits of the TNZ were found
 to be 17.7 ± 1.6°C and 34.5 ± 0.7°C, respectively, and the basal metabolic
 rate was 0.96 ± 0.12 ml O2 min-1 on average. The evaporative water loss
 (EWL) inflection point was established at 31.85 ± 0.27°C and was not
 significantly different from the value of the upper critical limit. No
 significant relationship was observed between the relative bill or tarsi
 size and TNZ critical limits, breadth, mass-independent VO2 or
 mass-independent EWL at any environmental temperature (from 10°C to 40°C).
 However, Great tit males (but not females) with larger tarsi areas (a
 proxy of leg surface area) showed higher cooling efficiencies at 40°C. We
 found no support for the hypothesis that the bill surface area plays a
 significant role as a thermal window in Great tits, but the leg surface
 areas may play a role in males’ physiological responses to high
 temperatures. On the one hand, we argue that the studied population
 occupies habitats with available microclimates and fresh water for
 drinking during summer, so active heat dissipation by EWL might be favored
 instead of dry heat loss through the bill surface. Conversely, male
 dominance behaviors could imply a greater dependence on cutaneous
 evaporative water loss through the upper leg surfaces as a consequence of
 higher exposure to harsh environmental conditions than faced by females.
description: All procedures were approved by the bioethical committee of the
 University of Extremadura, Spain (108/2016) and were conducted under the
 governmental license CN0032/18/ACA. <b>Capture
 and biometric measurements</b> The Great tit
 individuals examined belonged to a population located in the areas
 surrounding the city of Badajoz (SW Spain; 38° 56’ 7.85” N, 6° 56’
 33.129”). Classified as Csa according to the Köppen Climatic
 Classification (mean annual T<sub>a</sub>: 17.27 ± 0.05°C and
 summer mean maximum T<sub>a</sub>: 34.18 ± 0.07°C; data from
 1998 to 2018, State Meteorological Agency), this area has experienced a
 significant increase in the summer maximum T<sub>a</sub> and
 frequency and duration of heat waves over the last three decades (Acero et
 al. 2014, 2018). A total of 24 Great tits were collected as nestlings and
 hand-raised in the laboratory of the University of Extremadura during
 spring 2017. Individuals were maintained in artificial nests where they
 were fed every 2 h from 8:00 AM to 10 PM until fledging. When birds were
 completely independent, they were individually identified with an
 alphanumeric band and moved to outdoor aviaries (5 m × 2.5 m × 2 m each)
 equipped with small ponds with running water, natural vegetation and live
 prey where they stayed during several months before metabolic trials
 started (winter 2019). Taking advantage of a parallel study, we added data
 from 8 wild-living individuals (5 juveniles and 3 adults) caught in the
 same study area to optimize our sample size of individuals representing
 each sex. These Great tits were captured in the wild by mist-nets in the
 late afternoon, measured at night (winter 2019), and released early the
 next morning. The age and sex of the wild-living birds
 were determined according to their plumage characteristics (Svensson
 1992), and the sex of all individuals (16 males and 16 females) was later
 confirmed by CHD-based molecular sexing protocols (Griffiths et al. 1998).
 Hand-raised birds were released at their place of collection several weeks
 after respirometry measurements. Bill and tarsi surface
 areas were estimated individually following Greenberg et al. (2012).
 Briefly, we used an equation in which the bill area is approximated to an
 elliptical cone: where BW is the bill width, BD is the bill depth, and BL is the bill length (see figure in Svensson 1992). Measurements of the tarsus were used to estimate the tarsi surface area as a proxy of the leg surface area using the equation for an elliptical cylinder: where TW is the tarsus width, TD is the tarsus depth (both measured at the midpoint of the tarsus) and TL is the tarsus length. All bill and tarsus measurements were performed by the same person (JMAG) using a digital calliper (± 0.01 mm). We also measured the wing length (flattened and straightened) as a proxy of body size using a wing rule (± 0.5 mm) (Gosler et al. 1998). <b>Gas exchange measurements</b> We measured O<sub>2</sub> consumption (ml min<sup>-1</sup>) and EWL (mg hr<sup>-1</sup>) using an open flow-through respirometry system. Each individual was placed in a polypropylene metabolic chamber (232 × 165 × 162 mm; effective volume = 3.9 L), the floor of which was covered with a 1cm mineral oil layer to avoid evaporation from excreta. The chambers were equipped with a wire mesh platform located 3 cm above this oil layer to allow individuals to perch without touching the oil. All metabolic chambers were placed in a temperature-controlled cabinet (ICP, 750 Memmert GmbH, Schwabach, Germany), where the increasing or decreasing T<sub>a</sub> profiles (see below for details) were created automatically using control software. We introduced a calibrated thermistor probe (± 0.001°C) inside the metabolic chambers to monitor the T<sub>a</sub> during the metabolic trials. Exterior dry air (&lt;1 kPa WVP) was pumped from an air dryer compressor (MESTRA®) into a carboy (Lighton, 2008) and then directed to the metabolic chambers using mass flow controllers (MFS, Sable Systems International, Las Vegas, Nevada, USA). Flow rates of 1000 or 3000 ml min<sup>-1 </sup>(depending on data collection protocol; see details below) were used during metabolic trials. Excurrent airstreams from the chambers flowed through an eight-channel multiplexer (RM-8, Sable Systems International), which automatically alternated every 360 sec between metabolic chambers containing birds as well as an additional chamber left empty to obtain baseline values. The latter were obtained for 300 sec at the start of every trial and following two metabolic chamber measurements. We subsampled the downstream air at 200 ml min<sup>-1</sup> (SS3 subsampler, Sable Systems International) and pulled it sequentially through an H<sub>2</sub>O analyzer (RH300, Sable Systems), a Drierite® column and an O<sub>2</sub> analyzer (FC-10 Oxygen Analyzer, Sable Systems). The data were digitalized using an analog-to-digital converter (UI2 model, Sable Systems) and recorded with a sampling interval of 1 s using Expedata software (version 1.9.14, Sable Systems). Both analyzers were zeroed and spanned weekly using standard protocols (Lighton, 2008). <b>Data collection protocol</b> Gas exchange rates were measured through a wide range of T<sub>a</sub>s (10°C, 15°C, 20°C, 25°C, 30°C, 35°C, 37°C and 40°C) in a stepped manner in a maximum of six individuals at a time. The metabolic trials were performed at night (from 8:00 PM to 8:00 AM; the daily resting phase of Great tit) after the food was withheld from the birds for at least 2 h to ensure they were in a postabsorptive state (RER 0.70). To determine T<sub>lc</sub>, six birds at a time were exposed alternately to an increasing or a decreasing stepped T<sub>a</sub> profile ranging from 10°C to 30°C or vice versa. All individuals were exposed to each T<sub>a</sub> for a minimum of 65 min using a flow rate of 1000 ml min<sup>-1</sup>. For T<sub>uc</sub> determination, individuals were exposed to an increasing profile of T<sub>a</sub>s (35°C, 37°C and 40°C). We used a flow rate of to 3000 ml min<sup>-1 </sup>to ensure maintenance of low humidity levels (&lt;1 kPa WVP), which aided in keeping birds calm (Whitfield, et al. 2015), and only two birds were measured per trial; they were exposed to each T<sub>a</sub> for a minimum of 25 min. The first 65 min or 25 of the stepped T<sub>a</sub> profiles of each protocol were used to ensure that the individuals were acclimated to the metabolic chambers (i.e., stable VO<sub>2</sub> and EWL traces) after handling. To ensure captive individuals recovered from the stress of handling following T<sub>lc</sub> measurement, bird exposition to the highest T<sub>a</sub>s (35°C, 37°C and 40°C) was conducted after two weeks. Individual’s behavior within chambers was monitored directly by an observer (we did not record videos) using infrared cameras to ensure they remained calm during the metabolic measurements. All individuals were hydrated and weighed (± 0.1 g) before and after the metabolic measurements. The mean body mass (M<sub>b</sub>) of the birds was used in the analyses. <b>Data analysis</b> The VO<sub>2</sub> and EWL values at each T<sub>a</sub> were estimated as the lowest stable 2-min (see, for example, Boratyński et al. 2016) values using Eqs. 10.2 and 10.9 from Lighton (2008), respectively, with a custom macro designed in Expedata. We used a respiratory quotient of 0.70 (e.g., Kvist and Lindsröm 2001). To obtain the metabolic heat production (MHP), and converted the VO<sub>2</sub> values to metabolic rates (Watt, W) using an energy equivalent of 20 kJ 1<sup>-1</sup> O<sub>2</sub> (e.g., Caro and Visser, 2009). The drift of water and O<sub>2</sub> traces was corrected using the Catmull-Rom spline correction applied to baselines. The evaporative heat loss (EHL) was calculated assuming latent heat of vaporization values for water at 35°C, 37°C and ­40°C following Tracy et al. (1980). The evaporative cooling efficiency (EHL/MHP) was calculated at every T<sub>a</sub> above T<sub>uc</sub>. We used a generalized estimating equations (GEE) approach to simultaneously identify population limits of TNZ (T<sub>lc</sub> and T<sub>uc</sub>) of our Great tit population (n=24) using the ‘lme4’ package (Bates et al. 2014), the 'geepack' package (Halekoh et al. 2006) and a modified version of the 'segmented' package (Muggeo, 2009), in R 3.6.1. Then, we calculated T<sub>lc</sub> and T<sub>uc</sub> for each focal individual using the R packages ‘lme4’ and the modified version of ‘segmented’. We could not obtain TNZ values for wild-living individuals since they were only measured under one of both protocols. The VO<sub>2</sub> values were corrected by body mass using residuals from a regression between VO<sub>2</sub> and body mass (log-transformed values). The TNZ breadth was calculated as the T<sub>uc </sub>value<sub> </sub>minus the T<sub>lc </sub>value. Mean value of VO<sub>2</sub> within TNZ was considered to be the basal metabolic rate (BMR). The inflection point of EWL was also calculated using ‘lme4’ and ‘segmented’ packages in R. To obtain the relative appendage sizes (bill and tarsi index values), we computed the residuals of the regression of the bill or tarsi surface area on the wing length, as this is assumed to be the best proxy<i> </i>of body size in small-sized passerines, including Great tits (Gosler et al. 1998, Gardner et al. 2016). We log-transformed the variables to meet the assumptions of linearity, homoscedasticity and normality. The residuals were calculated separately for males and females due to sexual dimorphism in size. Great tit males had larger wing lengths (<i>t</i><sub>30</sub> = -2.27, <i>p </i>&lt; 0.05), higher M<sub>b</sub> values (<i>t</i><sub>30</sub> = -2.38, <i>p</i> &lt; 0.05) and bill surface areas (t<sub>30</sub> = -2.16, <i>p </i>&lt; 0.05) than females, but the sexes did not differ in tarsi surface area (<i>t</i><sub>30</sub> = -1.17, <i>p </i>= 0.25). We calculated mass-independent VO<sub>2</sub> and mass-independent EWL from the regressions of the VO<sub>2</sub> and EWL rates, respectively, on the mean M<sub>b</sub>. To ensure that the order of exposure of each individual to T<sub>a</sub> did not affect the metabolic measurements from 10 to 30°C, we performed a <i>t</i>-test to compare the mass-independent VO<sub>2</sub> and mass-independent EWL values between individuals measured at the decreasing or increasing stepped T<sub>a </sub>profiles. No significant differences were found in the analysis (all results <i>p </i>&gt; 0.09), so the order of T<sub>a</sub> exposure was not considered in the models. To test the effects of the bill and tarsi indices on physiological traits (T<sub>lc</sub>, T<sub>uc</sub>, TNZ breadth, the T<sub>a</sub> inflection point of EWL, EHL/MHP and the mass-independent VO<sub>2</sub> and mass-independent EWL at each T<sub>a</sub>), we built a series of generalized linear models (GLMs) that included physiological traits as response variables, sex (two levels) as a fixed factor, bill and tarsi indices as covariates, and the interactions between the bill index and sex, and between the tarsi index and sex, as fixed factors. In the case of EHL/MHP, M<sub>b</sub> was included as a covariate. Multicollinearity was tested by calculating the variance inflation factor (VIF) among all predictor variables using the ‘car’ package (Fox and Weisberg, 2019); we confirmed no collinearity problems (all VIF values &lt; 5; see Zuur et al. 2010). The model selection was based on the Akaike information criterion for small sample sizes (AIC<sub>c</sub>) to identify the top model(s) (models within 2 ΔAICc of the top model), and the AICc weights (<i>w</i><sub>i</sub>) were used to further distinguish among the top models (Burnham and Anderson 2002). We used the function ‘dredge’ from the R package MuMIn (Barton 2018) for this procedure. In cases where more than one model had ΔAICc &lt; 2 but <i>w</i><sub>i</sub> &lt; 0.9 (Burnham and Anderson, 2002), we performed model averaging (Grueber et al. 2011). A predictor was considered significant when the 95% confidence interval (CI) for the estimated coefficient did not overlap zero. We further calculated the relative importance weight (RIW) of each explanatory variable (see Table 2). Statistical analyses were conducted in SPSS Statistics 23 (SPSS Inc., Chicago, IL, USA) and R 4.0.3 (R Core Team, 2014), and figures were produced using the R package ‘ggplot2’ (Wickham et al. 2016). Values are shown as means ± SEs.