Elevation and body size drive convergent variation in thermo-insulative feather structure of Himalayan birds

Globally, high elevation habitats have been independently colonized by taxa separated by millions of years of evolution. Mountains thus represent excellent systems to study how distantly related species adapt to the same environmental challenges. Cold temperatures influence the elevational distribution of birds along montane gradients. Yet the eco-physiological adaptations that may explain this pattern, such as variation in insulative feather structure across high elevation and low elevation species has not been quantified. We used a comparative approach to understand if elevation, evolutionary history and body size drive variation in thermo-insulative feather traits across 1715 specimens of 249 Himalayan passerines. Controlling for phylogenetic relationships between species, we found that the proportion of the feather’s plumulaceous (downy) section increased with elevation. Body size also had a predictable effect on thermo-insulative variables with small birds having relatively longer feathers and thus a more insulative plumage than large birds. We show that an increase in the proportion of the feather’s downy section at colder temperatures is an evolutionarily widespread response across temperate and tropical taxa, and overall, smaller-bodied birds tend to have longer and more insulative feathers. Our results reveal convergent patterns in feather structure variation as a response to cold temperatures across species separated by millions of years of evolution.


Introduction
The evolution of traits that allow taxa to cope with abiotic constraints may be driven by an organism's environment, the region of its evolutionary origin, or a combination of both (Miller et al. 2013, Khaliq et al. 2014. The relative importance of these effects can be examined by investigating species that show independent evolution of adaptations to similar abiotic constraints. Mountains show predictable changes in environmental variables along the elevational gradient and montane communities are composed of genera that have independently colonized, and evolved modifications to cope with the environmental challenges of high elevation habitats (Cadena et al. 2012, Quintero andJetz 2018). Mountains thus offer powerful systems to understand whether species that 681 are separated by millions of years of independent evolution, but currently share the same environmental constraints, have evolved convergent or divergent adaptations towards these abiotic challenges.
Environmental temperature is a significant abiotic constraint for birds (Swanson and Liknes 2006) and temperature regimes are considered an important factor directly or indirectly driving species elevational distributions on mountains (Jankowski et al. 2013). On both tropical and temperate mountains, birds have shifted their distributions upslope due to warming temperatures (Tingley et al. 2009, Freeman and Freeman 2014, Freeman et al. 2018. However, for individual species, the magnitude of the effect of temperature on distribution is often mediated through more complex relationships between temperature, habitat and interspecific competitors (Tingley et al. 2012, Elsen et al. 2016, Hanz et al. 2019. To disentangle the interactive effects of temperature and other biotic variables (Srinivasan et al. 2018) on elevational distribution of birds, a crucial first step is to uncover the eco-physiological adaptations that allow high elevation taxa to survive temperature extremes. Such a mechanistic understanding can reveal how temperature may affect not only current but also future elevational distribution of montane birds.
Avian thermoregulation is chiefly a combination of insulation by plumage, and heat production brought about by enhanced respiration through increasing the basal metabolic rate (Swanson and Liknes 2006). Research on the ability of various species for metabolic rate variation has demonstrated that tropical birds may not have the capacity to increase their basal metabolic rate significantly (Swanson 2001, Londoño et al. 2015, but may be able to cope with low temperatures by having cold tolerance traits like broad thermoneutral zones (the range of temperatures where endotherms can maintain internal body temperature with minimal metabolic regulation), and lower conductance (body heat lost to the environment) (Londoño et al. 2017, Pollock et al. 2019. Despite these findings, there's scant evidence for widespread macro-ecological patterns that reduce conductance in tropical birds, such as an increase in body size along elevational gradients (Bergmann's rule) (Freeman 2017).
Thermal conductance in birds is greatly reduced by the insulation provided by plumage (Scholander 1955). Differences in metabolic capacities for thermoregulation among temperate and tropical taxa may lead to the prediction that highland tropical species may be better insulated than sympatric temperate species to compensate for this difference in metabolic thermoregulatory capacity. However, adaptations for thermoregulation in temperate taxa are likely a combination of metabolic and plumage-driven changes and thus whether tropical taxa are better insulated than temperate species needs to be tested in a community where both lineages live in sympatry. Yet, there are very few multi-species studies exploring the association between environmental temperature and the variation in feather structures in terrestrial birds (Pap et al. 2017, Osváth et al. 2018. Furthermore, there are no studies documenting this variation across a mountain assemblage although such differences have been suggested previously (Bangs 1899). Knowledge of thermoinsulation provided by plumage thus represents a critical gap in our understanding of thermoregulation in birds in general and montane birds in particular.
Bird contour (body) feathers cover most of a bird's body (Stettenheim 2000, Gamero et al. 2015 and are made up of proximal plumulaceous henceforth downy, and distal pennaceous sections (Osváth et al. 2018). The downy section of the feather is important for thermo-insulation (Stettenheim 2000, Dove et al. 2007. Previous research has shown that, within species, feather structure is variable and driven by environmental conditions (Broggi et al. 2005, Vagasi et al. 2012. Overall, in comparison to passerines living in warm habitats, species living in cold regions have longer feathers, feathers with proportionately longer downy sections, and lower downy barb density (Wolf andWalsberg 2000, Pap et al. 2017). Longer, and thus more overlapping feathers increase the depth -the distance between the outside air and the bird's skin -of the plumage. Greater proportional length of the downy section increases thermo-insulative capacity of an individual feather. Insulative capacity is further enhanced by lower downy barb density that creates larger air pockets to trap body heat (Stettenheim 2000, de Zwaan et al. 2017). However, the nature of the relationship between barb density and thermo-insulation is likely quadratic, first increasing and then decreasing thermo-insulative capacity beyond a threshold of lower barb density, and has not been fully explored. Low downy barb density is also predicted to be associated with high pennaceous barb density, such that sparse downy barbs trap more air while dense pennaceous barbs work as a seal reducing heat loss to the environment (Butler et al. 2008).
The Himalayas offer an excellent system to understand how thermo-insulative properties of bird feathers may vary by 1) elevation, with a predictable decrease in temperature along the gradient, and 2) clade origin, the temperate or tropical evolutionary history of the taxon. There is decrease of about 25°C in mean annual temperature between 0 and 5000 m in the Himalayas (Supporting information). Higher elevations are more seasonal (difference between highest and lowest temperatures) and there is a strong precipitation gradient in the Himalayas with low elevations having higher precipitation than higher elevations (Price et al. 2011). These mountains support a hyper-diverse bird community found along the world's largest elevational gradient (Price et al. 2014). The high Himalayas have been independently colonized by hundreds of avian taxa separated by millions of years of evolutionary history (Price et al. 2011). White et al. (2019) showed that freezing temperatures may go beyond limiting individual species, but also shape the phylogenetic composition of entire avian communities in the Himalayas. In these mountains, cold temperatures explain quantitative differences in community composition and species assemblages above and below the freezing temperature line (more tropical warm-adapted species below the freezing line and more temperate or cold-adapted taxa above the freezing line). Thus, Himalayan bird communities are an admixture of both temperate and tropical clades (Päckert et al. 2012(Päckert et al. , 2020. We explored the relationship between insulative feather structures, elevation, clade origin and body size in a suite of Himalayan birds. Across 249 species (Supporting information) of Himalayan passerines, we predicted that 1) proportion of the downy section and relative feather length would increase with elevation (Butler et al. 2008, de Zwaan et al. 2017. Given that small animals have higher conductance than large ones because of a larger body surface area to volume ratio (Scholander et al. 1950), we predicted that 2) small bodied birds would have more insulative feather structures than larger species. Finally, given that temperate species have higher metabolic flexibility for thermoregulation (Wiens et al. 2010), we predicted that 3) tropical taxa should show a stronger association between insulative feather structures and elevation (but see Pollock et al. 2019). However, it is also plausible that the ability of temperate species to survive cold temperatures is due to a combination metabolic flexibility and increased thermo-insulation in which case, all taxa, irrespective of clade origin, would show changes in feather traits with elevation. Investigating these differences can help uncover the mechanisms that drive current differences in species elevational distribution.

Material and methods
We examined passerine bird specimens (n = 1715, 249 species) in the Smithsonian National Museum of Natural History (USNM, the Supporting information lists the specimen numbers used in the study) collections with known elevation of collection from the Himalayas (Fig. 1). We restricted our analysis to species known to breed within the greater Himalayas excluding any regions of the Tibetan plateau and winter migrants to the Himalayas (Price et al. 2014).

Feather structure variables
Following the non-invasive methodology outlined in Butler et al. (2008), we took photographs of individual feathers in the breast region isolated on the specimen using a light and stereo microscope with camera. Feathers from this region of the body have been used in several other studies investigating the thermo-insulative structure of feathers (Broggi et al. 2005, Osváth et al. 2018. Although there is variation in the structure of feathers in other parts of the body (e.g. the head), previous studies have shown that there is no significant difference in the thermo-insulative structures of feathers on the ventral and dorsal sides of the bird (Pap et al. 2020), regions where contour feathers cover the greatest proportion of body. Thus, we restricted our data collection to only ventral feathers of the specimen in our dataset. We measured the following feather structure variables from these photographs at the USNM using Leica Application Suite, LAS V2.6 image analysis software with image calibration for measuring length in millimeters. 1) Length of plumulaceous (downy) section of the feather: the downy section was designated as the region with non-interlocking fluffy textured barbs. 2) Length of pennaceous section: the region of the feather with inter-locking barbs. 3) Total length of feather: the sum of the downy and pennaceous sections. 4) Downy barb density: number of downy barbs per millimeter measured in the middle portion of the downy section of the feather. 5) Pennaceous barb density: number of pennaceous barbs per millimeter measured in the middle portion of the pennaceous section of the feather. 6) Downy barb length: mean length of four downy barbs in the middle section of the downy section that were not visibly damaged. 7) Pennaceous barb length: mean length of four pennaceous barbs in the proximal section of the feather that were not visibly damaged.
Based on these measurements we calculated 1) proportion of downy section of the feather, computed as the length of the downy section divided by the total length of the feather and 2) proportion of the pennaceous section of the feather: length of the pennaceous section of the feather divided by the total length of the feather. We first assigned body mass to each species (n = 249) as the median body mass provided in the Birds of the World (Billerman et al. 2020) and calculated relative feather length to body size (using body mass as a proxy for body size) by dividing the total feather length by the species body mass. To ensure repeatability of our measurements within an individual specimen, out of the 1715 specimens, we measured the same variables on two adjacent contour feathers of 479 specimens with repeatability across every species. We calculated an interclass correlation coefficient ICC for each variable using the package ICC (Wolak and Wolak 2015). We restricted our sampling to taxa with body size less than 200 g due to size constraints on the feathers we could accurately measure given our microscope apparatus.

Clade origin -temperate or tropical
Latitudinal mid-point of the range has often been used as a proxy for region of evolutionary origin for the species (Weir andSchluter 2007, Page andShanker 2020). We also used latitudinal mid-point of the global range as a proxy for region of evolutionary origin (clade origin) for each species. Clade origin was calculated by leveraging Bird Life International species distribution maps for all species in our dataset (BirdLife International and NatureServe 2015, Supporting information). For species that are latitudinal migrants we only used regions where the species breeds for the calculation. For species that have both a migratory and resident populations we chose the maximum and minimum breeding latitude as the greater or lesser of the two respectively. We classified a species as having a tropical clade origin if the latitudinal mid-point of distribution was between 0° and 23.4° North or South while a species with a latitudinal mid-point of distribution more than 23.4° North or South was classified as temperate (Supporting information).

Statistical analyses
We assembled a dataset that included the following variables for each specimen and each feather measured within a specimen 1) proportion of downy section, 2) proportion of pennaceous section, 3) downy barb density, 4) pennaceous barb density, 5) mean downy barb length, 6) mean pennaceous barb length, 7) log body mass for species, proxy for body size, 8) elevation of specimen collection (min = 75 m and maximum 4924 m) and 9) clade origin temperate or tropical based on latitudinal mid-point of species distribution.
To decide which dependent variables to retain in our analysis, we explored correlations among feather structure variables and only retained those in the final analysis which were not strongly and significantly correlated with others. These included three feather structure variables, 1) proportion of downy section ( Fig. 2A), 2) downy barb density and 3) total feather length relative to body size. However, we did not include downy barb density in the analysis since it had low repeatability in our dataset and opposing predictions exist about its relationship with cold weather in the literature (Pap et al. 2017, Butler et al. 2008). Since all feather structure variables are collected from the same feather, these measurements cannot be considered independent datapoints. Thus, to understand the relative importance of elevation, body size and clade origin on feather structure we implemented a phylogenetic multivariate generalized least squares (pmgls) model using the 'mvgls' function in the MvMORPH package (Clavel et al. 2015) on mean values of all variables for each species (n = 249). Since the number of species far exceeded the number of traits in our dataset, we used a log-likelihood error structure in our pmgls model as suggested by (Clavel et al. 2019). Using proportion of downy region and relative total feather length as simultaneous dependent variables as in a multivariate multiple regression, we tested the effect of the following predictor parameters 1) mean elevation, 2) median body size (median body mass), 3) clade origin (temperate or tropical) 4) an interaction between mean elevation and clade origin to test whether feathers of tropical and temperate taxa vary differently in relation to elevation and 5) an interaction between median body size and clade origin to test whether the effect of clade origin is significantly different on the variation in relative feather length. We calculated p-values for the effect of predictor variables using the 'manova.gls' function in the MvMORPH package. In the pmgls model, we controlled for phylogenetic relationships between species by including a phylogenetic correlation matrix derived from a subset of the phylogeny of all Himalayan birds provided in Price et al. (2014).
Proportion of the downy segment and relative feather length are predicted to increase with increasing elevation and decreasing temperatures. However, the nature of this relationship may vary among taxa within our dataset. Additionally, relative feather length may also simultaneously vary with body size. To examine the relationship between feather thermo-insulative structures and elevation in our dataset in more detail, we tested for correlation between 1) proportion of downy segment and 2) relative length of the feather with elevation at the level of the genus. Although the sample sizes and elevations sampled for each species varied considerably in our dataset, through this analysis we wanted to investigate those taxa whose feathers vary significantly with elevation. All analyses were done in the software package R ver. 4.0.2 (<www.r-project.org>).

Repeatability of measurements and correlation within feather structure variables
Our dataset included 152 species with a temperate and 97 species with a tropical clade origin. Mean elevation for a species was positively correlated with clade origin such that specimens of species with a temperate clade origin were recorded at higher elevations than specimens of tropical taxa (r = 0.60, df = 247, p < 0.01, Fig. 2B). However, there was no relationship between body size and elevation across the specimens in our dataset (r = 0.001, p = 0.9, Fig. 2C). Feather structure variables were highly repeatable within an individual and interclass correlation coefficients were comparable to other feather structure studies (Gamero et al. 2015, Osváth et al. 2018, length of downy section ICC: 0.99; length of pennaceous section ICC: 0.98; downy barb length ICC: 0.99; pennaceous barb length ICC: 0.99; downy barb density ICC: 0.58; pennaceous barb density ICC: 0.52. Downy barb density and pennaceous barb density have been shown to have similar levels of repeatability in other studies as well (Pap et al. 2017). Several feather structure variables were significantly correlated to each other in our large comparative dataset of terrestrial passerine birds (Fig. 2D-F); downy and pennaceous barb length (r = 0.96, p < 0.01), downy barb length and total feather length (r = 0.81, p < 0.01) and, contrary to prediction in the literature, downy and pennaceous barb density (r = 0.46, p < 0.01). These correlations helped us narrow down our list of dependent variables to the two variables 1) proportion of the downy section and 2) relative feather length used in the next section as they were not correlated to one another.

Associations between feather structure and ecological variables
Thermo-insulative feather structure variables were significantly associated with elevation and body size in Himalayan passerines. Our phylogenetic multivariate generalized least squares analysis revealed that, thermo-insulative feather structures had moderate phylogenetic signal λ = 0.51. We found that elevation was significantly positively associated with both proportion of downy segment and relative feather length (p = 0.02) while body size had a significant negative relationship (p < 0.01) with both these thermo-insulative variables (Table 1). Our results suggest that overall, high elevation birds have more downy and relatively longer feathers that are predicted to increase thermo-insulation in colder environments (Fig. 3A-B). However, this relationship with elevation is confounded by a negative association with body size such that small birds in general have more insulative feathers than large birds (Fig. 3C-D) and this relationship persists at high elevations. We did not find a significant effect of clade origin on the variation in thermo-insulative feather structure in Himalayan birds. Nor did we find a significant association between an interaction of clade origin and elevation or body size.
In addition to these general patterns, at the level of the genus, species from 17 out of the 107 genera (comprising 53 out of 249 species) showed statistically significant correlations between proportion of downy section and elevation. 14 out of these 17 genera showed the expected positive Table 1. Model estimates from a phylogenetic multivariate generalized least squares model that tested the relationship between thermoinsulative feather structures of Himalayan birds and various predictor variables. The p-values were calculated by conducting a multivariate analysis of variance (MANOVA) test on the generalized least squares object (Clavel et al. 2019  correlations with elevation (Fig. 4A). The direction of the correlations between relative feather length and elevation were mixed. Species from 22 genera showed significant correlation between relative feather length and elevation (comprising 68 out of 249 species). Eleven out of these 22 genera showed the expected positive association with elevation ( Fig. 4B) while the other 11 genera showed a significant negative relationship with shorter feathers at higher elevations.

Discussion
Our results demonstrate that across a diverse montane assemblage, species showed convergent variation in feather structure predicted to increase thermo-insulation. Previous studies of variation in feather structure have focused on temperate taxa (Gamero et al. 2015, Koskenpato et al. 2016, Osváth et al. 2018). Our study is the first to document feather adaptations to cold environments across a suite of species with both temperate and tropical evolutionary origins. These findings suggest that an increase in the proportion of the downy section may be a widespread response to a decrease in environmental temperature throughout passerines (Fig. 3A, Pap et al. 2017). Montane taxa showed 1) an increase in the proportion of down, and a small but significant increase in relative feather length, with an increase in elevation and 2) a significant decrease in the proportion of the downy region and relative feather length with an increase in body size. Birds show intraspecific variation in feather structure associated with environmental variables (Andreasson et al. 2020), a pattern demonstrated along latitudinal (Broggi et al. 2005) and elevational axes (de Zwaan et al. 2017). Pap et al. (2017Pap et al. ( , 2020 in a phylogenetic comparative approach, showed that feather structure variables known to vary within species, also vary across species, in association with environmental temperatures in Eurasian temperate birds. High elevation species are known to have lower conductance than low elevation congeners (Londoño et al. 2017). Our results demonstrate that such important eco-physiological differences between high and low elevation birds may be in part explained by the variation in insulative feather structures such as the proportion of the downy segment and relative feather length that we found in Himalayan birds.
Since we used museum specimens for our study, our non-destructive sampling method did not allow us to estimate feather density which has been shown to vary with temperature (Cheek et al. 2017, Osváth et al. 2018. We did, however, find evidence for increase in plumage depth through an increase in feather length, which also increases insulation. Downy and pennaceous barb length, and density, were significantly correlated across both tropical and temperate taxa (Fig. 2D-E). This finding was contrary to the prediction in the literature (Butler et al. 2008) that birds may increase insulation and reduce conductance by having low downy barb density that traps air pockets but high pennaceous barb density that seals these air pockets from escaping. Instead, the positive relationship between downy and pennaceous barb density suggests that birds with sparse downy barbs may reduce heat loss not by having dense pennaceous barbs, but by having deeper plumage. Downy barb length and total feather length were also significantly correlated (Fig. 2F) which, in conjunction with other results (Fig. 2E), demonstrated that birds with relatively long feathers also have relatively long downy barbs that can interweave extensively with barbs on adjacent feathers, creating an insulating air pocket. Such relationships among feather structure variables may be important for understanding the general role of feathers in insulating birds from cold temperatures.
Large-bodied animals are better at retaining body heat than small-bodied animals (Scholander et al. 1950). This predicts that small-bodied birds should show a stronger response to cold temperatures, in this case higher elevation, than large birds which may be better able to withstand cooler temperatures. Although evidence for within-species variation in body size along the elevational gradient to offset constraints of cold-stress is scant for birds (Freeman 2017, Sander and Chamberlain 2020), we found that thermo-insulative feather structure was significantly driven by body size. The weak association between relative feather length and elevation (Fig. 3B) suggests that even at low elevations, small birds have disproportionately longer feathers than larger birds. Both proportion of downy section and relative feather length were significantly negatively affected by body size. This suggests that large montane birds 1) may not rely on feather modifications for thermoregulation but use other eco-physiological routes, 2) can withstand cold temperatures disproportionately better than smaller taxa because of their smaller ratio of body volume to surface area or 3) use behaviors such as short-term elevational migration to escape harsh weather (Barçante et al. 2017), or a combination of the above to survive cold weather. Our results thus show a more general pattern in avian thermoregulation where feathers may play a bigger role in the thermo-insulation of small birds than larger taxa.
High elevation Himalayan bird communities are composed of more temperate, cold-tolerant species than tropical species (Elsen et al. 2016, White et al. 2019, Päckert et al. 2020). This pattern of evolutionary history influencing elevational distribution on mountains is not unique to birds or to the Himalayas (Fjeldså et al. 2012), and has been show in other taxa including mammals (Patterson et al. 1998), plants (Griffiths et al. 2021) and butterflies (Acharya and Vijayan 2015). In fact, cold-tolerant, high-latitude plants being distributed in high elevation tropical habitats has been known for almost two centuries (Von Humboldt 1849, Wallace 1869). This turnover in community composition to more cold-tolerant species at high elevations could be explained by the ability of temperate species to more efficiently thermoregulate at cold temperatures than tropical lineages (Wiersma et al. 2007, Wang et al. 2016. While temperate avian taxa may have more metabolic flexibility (Londoño et al. 2015) than tropical species, we did not find evidence that temperate species also show a stronger association between feather thermoinsulative structure and elevation than their high elevation tropical counterparts. Insulation provided by feathers may thus play an important role in allowing tropical taxa to have broad thermoneutral zones and live at high latitudes and elevations (Khaliq et al. 2015). The physical environment, and not evolutionary relationships, driving adaptations to cope with other abiotic constraints has also been shown for environmental hypoxia, another important environmental constraint in Himalayan birds (Barve et al. 2016).
Thermoregulatory ability is just one of many factors, biotic and abiotic, and the interactions between the two, that ultimately drive the elevational distribution of a montane bird (McCain 2009, Cadena et al. 2012, Jankowski et al. 2013, Barve and Dhondt 2017. Avian thermoregulation in itself is complex, ranging from behavioral (Barçante et al. 2017) to molecular adaptations (Stager and Cheviron 2020) for surviving cold stress. Even within the genera we sampled, there was large variation in the relationship between thermoinsulative feather structure variables and elevation. While several genera from both tropical (e.g. Pellorneum, Acridotheres, Arachnothera) and temperate (e.g. Parus, Luscinia, Anthus) lineages showed positive associations between elevation and feather structure, some genera (e.g. Hemixos, Pnoepgyga, Liocichla, Yuhina) actually showed trends opposite to those predicted to increase insulation along the elevational gradient ( Fig. 4A-B). The lack of a pattern in many species may be influenced by small sample sizes and limited elevation sampling for several genera. A lot remains to be investigated about how montane birds use their feathers to survive cold temperatures. Experimental studies quantifying differences in the insulative ability of feathers with varying amounts of down, and investigating both intraspecific and intrageneric variation, and broader comparisons between birds from other high mountains systems such as the Afromontane birds, Andean birds and the Himalayas may be fruitful next steps.
Himalayan birds show a trend where only a small subset of the tropical taxa found in the warm and humid south-east of the mountain range extend their distribution to higher elevations and the more seasonal north-west (Price et al. 2011, Srinivasan et al. 2014). However, the eco-physiological traits that may drive this pattern have not yet been investigated. Feathers play a large role in reducing the loss of body heat to the environment and are thus an integral part of avian thermoregulation (Stettenheim 2000, Dove et al. 2007, Du Plessis et al. 2012. We find that the tropical species that have been successful in expanding to high elevations have feather modifications convergent with those of temperate species. In our analysis, genera with high species diversity in the tropics, but that extended their distribution into the temperate and high-elevation Himalayas, such as Trochalopteron, Pomatorhinus and Pericrocrotus showed potentially adaptive variation in feather structure (Fig. 3A-B). This suggests that despite thousands of species independently colonizing montane habitats around the world, there may be convergent modifications in morphology such as an increase in the proportion of down and relative feather length arising repeatedly as species move to colder habitats. Our findings also underscore the importance of studying metabolic rate variation and feather structure variation in tandem to get a more complete understanding of cold tolerance in birds. Understanding thermoregulation by montane taxa is now more critical than ever because not only are mountains warming but extreme warm and cold events are increasing in frequency and severity around the world (Rangwala andMiller 2012, Wang et al. 2016).

Data availability statement
All relevant data are available as Supporting information.