Do New World pitvipers “scale‐down” at high elevations? Macroecological patterns of scale characters and body size

Abstract Bergmann's rule describes the macroecological pattern of increasing body size in response to higher latitudes and elevations. This pattern is extensively documented in endothermic vertebrates, within and among species; however, studies involving ectotherms are less common and suggest no consistent pattern for amphibians and reptiles. Moreover, adaptive traits, such as epidermal features like scales, have not been widely examined in conjunction with Bergmann's rule, even though these traits affect physiological processes, such as thermoregulation, which are hypothesized as underlying mechanisms for the pattern. Here, we investigate how scale characters correlate with elevation among 122 New World pitviper species, representing 15 genera. We found a contra‐Bergmann's pattern, where body size is smaller at higher elevations. This pattern was mainly driven by the presence of small‐bodied clades at high elevations and large‐bodied clades at low elevations, emphasizing the importance of taxonomic scope in studying macroecological patterns. Within a subset of speciose clades, we found that only Crotalus demonstrated a significant negative relationship between body size and elevation, perhaps because of its wide elevational range. In addition, we found a positive correlation between scale counts and body size but no independent effect of elevation on scale numbers. Our study increases our knowledge of Bergmann's rule in reptiles by specifically examining characters of squamation and suggests a need to reexamine macroecological patterns for this group.

Squamate reptiles (i.e., lizards and snakes) do not exhibit consistent patterns of body size along environmental gradients, and a single mechanism may not explain this variability (Olalla-Tárraga et al., 2006;Watt et al., 2010). Beginning with one of the earliest macroecological studies on lizards, Bogert (1949) found that larger species inhabited warm, low elevation areas, while smaller species were found in cooler, higher elevations (i.e., a contra-Bergmann pattern), which was recently observed in Sceloporus (Oufiero, Gartner, Adolp, & Garland, 2011). Across several families of snakes, Lindsey (1966) observed a slight tendency of larger species inhabiting regions with lower temperatures. Later, Reed (2003) found little support for Bergmann's rule in either Elapidae or Viperidae. Finally, a review by Millien et al. (2006) highlighted that snakes were the vertebrate group with the lowest agreement with Bergmann's rule. This deviation from the normal pattern of Bergmann's rule among snakes might be due to their elongated bodies that impact heat exchange through relatively high surface area to volume ratios (Feldman & Meiri, 2014;Lillywhite, 1987;Lindsey, 1966;Olalla-Tárraga & Rodríguez, 2007).
Squamates in particular possess keratinized scales that may play important roles in water balance presumably driving strong correlations between elevation and scale numbers, such as patterns found in Sceloporus (Acevedo, 2009). Variation in climate could select for differences in the size, shape, number, color, and perhaps other features of scales. For example, large scales are related to latitude (Bogert, 1949;Oufiero et al., 2011) but also to larger species of lizards. In addition, scale number can vary intraspecifically along altitudinal gradients (e.g., Anolis marmoratus, Malhotra & Thrope, 1994).
Other environmental factors like precipitation may also influence scale characters. In lizards, larger and fewer scales are observed in hot, dry regions (Hellmich, 1951;Horton, 1972;Lister, 1976;Sanders, Malhotra, & Thorpe, 2004), while snakes have more scales in these types of habitats (Brown, Thorpe, & Baez, 1991;Klauber, 1941Klauber, , 1997Licht & Bennett, 1972;Malhotra & Thrope, 1994;Soulé & Kerfoot, 1972;Thorpe & Báez, 1987). One hypothesis is that heat and water balance is related to the amount of exposed interstitial skin, which is influenced by scale size or number (Pough et al., 2001;Sanders et al., 2004). Alternatively, environmental conditions such as temperature during development may influence scale characters (Osgood, 1978). Body size is also often positively correlated with the number of ventral scales (Klauber, 1945;Lindell, Forsman, & Merila, 1993) while the distribution of body size has also been shown to vary along environmental gradients (Bogert, 1949;Pincheira-Donoso et al., 2008). Additional studies are needed to understand the range of patterns and the complex interplay between environmental gradients, body size, and scale characteristics in squamate reptiles in order to test potential mechanisms.
An excellent model to study macroecological patterns in reptiles is New World pitvipers, an ecologically and morphologically diverse, broadly distributed group with limited dispersal for a vertebrate. As carnivores, studying these snakes also reduces the influence of varying trophic level (Reed, 2003). Furthermore, the lineage is monophyletic, meeting a key requirement of Bergmann's rule (Cruz et al., 2005;Reed, 2003) and has a well-resolved phylogeny supporting comparative analyses among clades (Castoe & Parkinson, 2006;Gutberlet & Harvey, 2004;Jadin, Smith, & Campbell, 2011;Kraus, Mink, & Brown, 1996). In this study, we use information on evolutionary relationships and previously published data on geographic distribution and morphological traits to assess correlations between scale counts, body size, and elevation across lineages of New World pitvipers. We investigate whether a pattern emerges that suggests environmental change across elevation gradients produces physiological constraints and, therefore, selects for body size and scale characters. We predict that increasing elevation will be associated with smaller body size. This extends previous work investigating snake body size and elevation while also investigating important morphological traits (Reed, 2003;Terribile et al., 2009).
In addition to their broad geographic and ecological distribution, New World pitvipers have extensive morphological diversity (see F I G U R E 1 Two distantly related species, Mixcoatlus browni (a, b) and Crotalus intermedius (c, d), showing similar body size (maximum TL equals 51.5 and 57 centimeters, respectively), color pattern, and scale numbers (e.g., few scales on the head). This contrasts the distinction between C. intermedius (c, d) and its closer relative C. molossus (e, f), which lives in lowland, arid areas and has a larger body size (max TL 133 cm), different color pattern, and a greater number of scales. Photographs were taken by RCJ (b, d, e), Eric N. Smith (a, c), and Jonathan A.

| Clade assessment
Essential to understanding the patterns of species distributions along environmental gradients requires incorporation of phylogenetic information in data analysis (Gaston et al., 2008;Harvey & Pagel, 1991). Importantly, our understanding of the diversity and evolutionary relationships of pitvipers has become quite robust over the past two decades with a strong congruence between gross morphology and mitochondrial genes (see review in Gutberlet & Harvey, 2004). Phylogenetic relationships within most of the New World pitviper genera are strongly supported, and individual species appear to be accurately assigned to their respective genera.
For example, there is strong support for relationships within the Lachesis (Bushmasters). This lack of phylogenetic resolution concerning how the genera are related constitutes a large knowledge gap that hinders phylogenetic comparative analyses. Therefore, we grouped the species of New World pitvipers within their respective genera as well as some strongly supported clades composed of several closely related genera (e.g., Porthidium group) for statistical analyses. We conducted our analysis on several of these smaller clades and within particular genera to examine whether or not patterns were apparent across different taxonomic scales (Meiri & Thomas, 2007).

| Body size, scale counts, and elevation data
We used the literature to obtain data on scale morphology, body size, and elevation for 122 of the 150 currently described species of New World pitvipers (Uetz, 2019). These data were generally taken from Campbell and Lamar (2004) but data from several recently described or revised taxa were obtained from other published sources (Table 1). For each species, we derived a single value for each scale character of interest, obtained from data within the geographic range of the taxa (Gaston et al., 2008). This method is unaffected by species richness and is conservative by only considering each species once as opposed to grid-based methods (Meiri & Thomas, 2007). Specifically, we used the maximum total length recorded as our measurement of body size as this is recommended for use in snakes because of their slender, elongated bodies, and average body size is unavailable in the literature (Terribile et al., 2009). For squamation, we recorded number of intersupraoculars, mid-dorsal scale rows, subcaudals, supralabials, and ventrals (as described in Klauber, 1997). Most of these values were found with a range of values (e.g., ventrals), and we therefore recorded the mode (as provided in Campbell & Lamar, 2004)

| Statistical analyses
To test for a relationship between maximum body size and elevation, we used a general linear mixed effects model in R with the package "nlme" (Pinheiro, Bates, DebRoy, Sarkar, & R Core Team, 2014; R Core Team, 2013). We included an initial random intercept and slope structure of genus nested within clade to account for nonindependence of each observation due to species relatedness. In other words, the model accounts for the fact that the body sizes of snakes within a genus (and clade) are likely more similar to each other than to a random species outside of this genus (or clade). We then used likelihood ratio testing to identify the best random structure for the model. We chose this nested random effect method of accounting for nonindependence among snake species instead of using a linear model with phylogenetically independent contrasts because we did not have a fully resolved molecular phylogeny for the entire group of snakes. For the analysis, we log 10 -transformed maximum body size and square-root transformed midpoint elevation in order to meet the assumptions of normally distributed residuals, and we therefore used a Gaussian family error distribution with identity link. We also conducted separate regression analyses for the three most speciose clades: the Porthidium group, the rattlesnakes, and the South American lanceheads.
We then determined how scale counts were related to body size and elevation. First, in order to collapse the highly correlated data on scale counts, we performed a principal components analysis (PCA) on the log-transformed (centered and scaled) scale count data. The first axis of the PCA explained 44% of the variation in the data set, and the loadings from this axis were used as a proxy for overall scale counts. This first axis was positively correlated with all scale characteristics except the modal supralabial counts, allowing us to easily interpret higher PCA scores as representing generally higher scale counts ( Figure 2). To statistically model the effects of elevation and body size on scale counts, we again used a general linear mixed effects model with the same initial random effects structure and data transformations as above. We also calculated the variance inflation factor (VIF = 1.12), which verified that the colinearity between elevation and maximum body size was not strong enough to bias the model (Fox, 2008).

| RE SULTS
We found that, when we look across all clades, maximum body size significantly declines with elevation in pitvipers, displaying a contra-Bergmann's rule pattern (overall effect of elevation: types, suggesting that scale counts generally increase with body size. This positive association was clear both among and within clades and genera (Figures 6 and 7). In this case, the model that included random slopes and intercepts of clade and genus was indistinguishable from one that only included random intercepts. This suggests that the relationship (i.e., slope) between scale count and body size is similar among all pitvipers, though slight differences may exist among and within clades (Figure 7).

| D ISCUSS I ON
Across 15 genera of New World pitvipers, the data support a negative relationship between elevation and body size and a positive relationship between scale counts and body size (Figures 3 and 6), although elevation does not seem to directly correlate with scale counts after accounting for body size. The negative relationship between body size and elevation was clearly supported across taxonomic levels. When we examined patterns among genera, we observed larger bodied clades (e.g., Lachesis) at low elevations and smaller bodied clades (e.g., Mixcoatlus) at high elevations. Among the more speciose clades examined separately, only the genus Crotalus, which spans the largest range of elevation, was found to have a F I G U R E 5 The relationship between maximum body size and elevation for each genus. Unique intercepts result from the significant random effect on the intercept, where genus is nested within clade. The gray lines represent the 95% confidence prediction intervals taking into account the uncertainty in the intercepts and the overall slope. Clades are represented by the number preceding the genus labels Research investigating how body size changes along environmental gradients in squamate reptiles has found no consistent pattern (Oufiero et al., 2011). Some studies identified negative correlations between body size and latitude and elevation (Ashton & Feldman, 2003), while others have found positive relationships with latitude and no relationships with elevational range (Cruz et al., 2005). Within some studies, effects of latitude were stronger when elevation was taken into account (Cruz et al., 2005;Olalla-Tárraga et al., 2006;Terribile et al., 2009). Our study contributes to this growing body of research, specifically for elevational gradients, by supporting a contra-Bergmann's pattern across 15 genera in the New World pitvipers. Our results are similar to those found by Reed (2003) but in direct contrast to findings of Olalla-Tárraga et al. (2006) where the largest species of snakes occurred at higher, colder elevations in Western North America. Different patterns across studies could also be based on the geographic ranges of species studied. For example, Terribile et al. (2009) found that vipers follow a contra-Bergmann's pattern when analyses were restricted to South America but not North America. Similarly, in Europe snakes increased in body size with decreasing latitude, but the pattern was inconsistent and more complex in North America (Olalla-Tárraga et al., 2006). Additional discrepancies among studies may be related to different phylogenetic levels displaying significant variation in patterns of body size (Cruz et al., 2005;Olalla-Tárraga et al., 2006;Terribile et al., 2009). For example, eastern clades of Crotalus viridis sensu lato followed Bergmann's rule, while western clades showed the opposite trend (Ashton, 2001;Cruz et al., 2005). Different patterns of body size distributions across clades may also be the result of different physiological mechanisms (Terribile et al., 2009). Some of the most common mechanisms hypothesized to contribute to Bergmann's rule include heat conservation (Ashton & Feldman, 2003;Cowles, 1945;Olalla-Tárraga et al., 2006), embryonic development (Angilleta, Niewiarowski, Dunham, Leaché, & Porter, 2004;Oufiero et al., 2011;Tousignant & Crews, 1995), prey availability (Ashton & Feldman, 2003), and sexual selection (Pincheira-Donoso et al., 2008). However, our analyses were focused on correlations and could not distinguish among these hypotheses, making this an important subject of further investigation.

F I G U R E 7
The relationship between scale counts and maximum body size. Unique intercepts and slopes result from the random effects, where genus is nested within clade. However, it should be noted that, for this model, there was no clear difference between a model with no random slope. The gray lines represent the 95% confidence prediction intervals taking into account the uncertainty in the intercepts and the overall slope. Clades are represented by the number preceding the genus labels Physiological mechanisms involving temperature may also contribute to variation in other physical features and display variation along elevational gradients (Bergmann, 1847;Cowles, 1945;Watt et al., 2010). As an important first step, we described the pattern of scale counts with respect to both body size and elevation because of the potential role scales play in physiological processes of reptiles. For example, small Sceloporus occurring at high elevations have numerous small scales, while larger bodied, low elevation species have fewer larger scales (Bogert, 1949;Oufiero et al., 2011). It is hypothesized that, for lizards, having smaller, more numerous scales in colder environments facilitates heat retention, while the larger scales seen in warmer climates function as heat shields (Oufiero et al., 2011;Regal, 1975 Sanders et al., 2004;Soulé & Kerfoot, 1972). This is supported by patterns of Sceloporus with fewer scale rows found with decreasing aridity (Oufiero et al., 2011) and higher numbers of ventral scales in snakes from drier, open habitats (Martínez-Feíra et al., 2009). Another key aspect to scales related to variation in habitat is locomotor performance (Kelley, Arnold, & Gladstone, 1997). Snakes may have a closer functional association between scale counts and locomotor ability due to the elongation of the body and loss of limbs (Kerfoot, 1970 with stronger natural selective pressure on scale counts and body size than climate variation associated strictly with elevation alone. Our results suggest important potential mechanisms for future investigation but should be interpreted in light of the limitations of data collection and study design. First, our approach limits environmental gradients to a single point in geographical or environmental space and ignores the interaction of these factors (Blackburn & Ruggiero, 2001;Gaston et al., 2008;Olalla-Tárraga et al., 2010;Ruggiero & Lawton, 1998). In addition, although elevation as a proxy is extremely common in this literature, this is an assumption that might not show strict associations in all cases (Oufiero et al., 2011).
We also did not include any information about latitude in our analysis, which further affects climate. Future research could investigate the patterns of pitviper body size and scale traits with climatic factors directly, such as georeferenced snake specimens that can be directly linked to a local climate. Second, by analyzing single data points for species level and above, we cannot address how population-level variation in both body size and scale numbers influences the associations with elevation or climate (Martínez-Feíra et al., 2009). Third, our general approach is based on availability of data from the literature, and some sample sizes from which the values were taken may not accurately represent the taxa studied (Gaston et al., 2008). Finally, the marginal R 2 value was low (.04) for the relationship between elevation and body size, suggesting that there are additional factors driving body size distributions beyond elevational ranges.
Overall, our results contribute to the generality of documented patterns of body size clines in reptiles in the western hemisphere and our understanding of evolutionary and ecological mechanisms shaping reptile species richness, distribution, and assemblages (Adams & Church, 2008;Ashton & Feldman, 2003;Cruz et al., 2005;Gaston et al., 2008;Olalla-Tárraga et al., 2006). By examining scale traits in addition to body size, our study also has important implications for taxonomy. Among widely distributed taxa over diverse habitats, confusion of taxa may exist because most species descriptions depend on morphological traits (Sanders et al., 2004). For example, similar selective pressures can result in similar morphological traits of closely related species of viper (Martínez-Feíra et al., 2009). Our research also has implications for understanding the role of climate change on pitvipers due to potential shifts in elevational ranges and decreasing body size (Cruz et al., 2005;O'Brien, Fox, Planque, & Casey, 2000;Portner, 2001). Better understanding of the relationship between climate and morphology may support the development of new hypotheses pertaining to the response of these taxa to climate change and enhance conservation efforts (Oufiero et al., 2011;Sinervo et al., 2009).

ACK N OWLED G M ENTS
Only now after significant investment in general systematic and natural history studies of biodiversity are we able to investigate and understand broader patterns of ecology and evolution. We thank those involved in funding and working in natural history collections that make these insights possible. We thank Jonathan A. Campbell and Eric N. Smith of the Amphibian and Reptile Diversity Research Center at the University of Texas at Arlington for access to specimens and photographs.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
RCJ conceived the study, while SAO and JRM improved on those ideas; RCJ collected the data; JRM, SAO, and RCJ analyzed the data; All authors contributed to the writing of the paper and have approved of the final version.