Evolution of Pelage Luminance in Squirrels (Sciuridae)

Introduction

Gloger’s rule is an ecogeographic rule describing the tendency of animals living in warm and wet environments to be darker (Gloger, 1833; Rensch, 1929). This rule has been supported by several quantitative tests from a wide variety of taxa including birds, mammals, reptiles, amphibians, and insects (Delhey, 2019; Delhey et al., 2019; Stanchak and Santana, 2019). There are multiple mechanisms that have been proposed to explain the color patterns that conform to Gloger’s rule (Delhey, 2019) with the most common being that darker coloration is selected for in environments with lower light conditions to facilitate greater concealment from predators. Dark environments are often created by greater vegetative production, which is generally driven by higher temperatures and environmental moisture. Moreover, a high number of studies have found support for the UV protection hypothesis (Chaplin, 2004) by finding latitude, a strong correlate of solar radiation, to be significantly associated with pelage luminance (e.g., Stoner et al., 2003a, b, Caro, 2005; Santana et al., 2012). Thus, biologists that study Gloger’s rule often investigate how well multiple environmental variables (i.e., temperature, environmental moisture, primary productivity, solar radiation, and their surrogates) predict the brightness of animal color. Another possible mechanism is selection for pelage color to match that of the background substrate (Nachman et al., 2003; Vignieri et al., 2010), the color of which may be predicted by the variables mentioned above.

The rodent family Sciuridae is a fascinating group to study color-pattern evolution because, unlike most other rodent families, squirrels have evolved extremely diverse pelages that include many types of patterns (e.g., stripes, spots, and patches) and hues (e.g., red, brown, tan, and black), and a range of luminances from black to very light (Figure 1) (Thorington et al., 2012). This phenotypic diversity is likely related to the evolution of squirrels in many different terrestrial biomes worldwide (e.g., desert, tundra, and rainforest), as well as their evolution of a diverse range of lifestyles (terrestrial, arboreal, and fossorial), life histories (diurnal, nocturnal, social, and solitary), and body sizes (16 g up to 1200 g) (Thorington et al., 2012). For example, many species of desert squirrels have pale-colored fur that is likely related to living in habitats with lightly colored soils that contain minimal amounts of dark organic content.

Figure 1. (A) Average dorsal color of all species for which we have photo and distribution data superimposed onto that species range map. (B) Dorsal regions of select squirrel species showing the distribution of pelage luminance across the Sciurid family. The representative species from left to right are as follows: Callosciurus finlaysonii, Sciurus aestuans, Sciurus vulgaris, Sundasciurus samarensis, Callosciurus nigrovittatus, Sciurus yucatanensis, Notocitellus adocetus, Petaurillus kinlochii, Urocitellus townsendii, Callosciurus finlaysonii.

Patterns of Gloger’s rule have been investigated in many mammal groups, but a large majority of these studies have been focused at the intraspecific level (Delhey, 2019). Studies of Gloger’s rule in squirrels are limited, but are also focused on intraspecific variation. A study in fox squirrels (Sciurus niger) with melanistic polymorphisms has identified elevation and humidity as positively correlated with pelage luminance, but also showed that precipitation and temperature were negatively correlated with pelage luminance (Kiltie, 1989). They also found fire frequency was positively correlated with darker pelage. A similar study of the African giant squirrel (Protoxerus stangeri) found a negative relationship between pelage luminance and mean annual precipitation and a positive correlation between pelage luminance and mean annual temperature.

Studies of Gloger’s rule in primates offer potentially intriguing insight to selective pressures that may drive color evolution in Sciurids because both mammalian groups have a wide diversity in color patterns and are similar in many other characteristics, such as possessing a mostly diurnal behavior, a range of both terrestrial and arboreal lifestyles, and a global distribution. A comparative study of Gloger’s rule in primates showed that pelage luminance was negatively correlated with actual evapo-transpiration (AET), which is a measure of the amount of water leaving an environment (Kamilar and Bradley, 2011). As AET is a composite measure of many environmental processes, this result implied a negative correlation between pelage luminance and precipitation, temperature, humidity, and primary productivity.

In this study, we used a comparative phylogenetic approach to investigate environmental predictors of dorsal-pelage luminance in Sciurids. We hypothesized that, in accordance with many other animal systems that conform to predictions by Gloger’s rule, pelage luminance will be significantly associated with variables related to environmental moisture, temperature, and light conditions. We also hypothesized that pelage luminance will be significantly associated with other environmental variables that are more directly associated with background matching (soil carbon content, soil density, and frequency of wildfires).

Materials and Methods

Phylogenetic Inference

To conduct phylogenetic comparative analyses, we constructed a phylogeny of 137 Sciuridae species across 51 genera. We downloaded all gene sequences from GenBank representing seven mitochondrial loci (12S, 16S, Cyt-b, COI, COII, COIII, mitochondrial control region) and four nuclear loci (C-myc, GHR, IRBP, RAG1). We then aligned all sequences using MUSCLE (Edgar, 2004) implemented in the EMBL web app (Madeira et al., 2019) and then manually edited egregious alignments. We then used Partition Finder (Lanfear et al., 2012) to find the best substitution model based on the Akaike Information Criterion. Finally, we inferred a phylogeny using Bayesian inference in MrBayes (Huelsenbeck and Ronquist, 2001) using Markov Chain Monte Carlo (MCMC) sampling. Two different runs (each with one cold and three heated chains) were analyzed for 3 × 10⁷ generations (with trees sampled every 100 generations). Branch lengths were than calibrated according to the dates used by Zelditch et al. (2015).

Luminance Data From Museum Study Skins

To obtain color measurements of the 137 squirrels used in our phylogenetic analysis, we photographed the dorsal profile of museum study skins from the Division of Mammals at the Smithsonian National Museum of Natural History. We used a Sony Alpha 5000 to take digital photographs of specimens against a gray matte background. Two light-boxes containing ESDDI 85 W 5500K Day Light Florescent bulbs were positioned at opposing ends of the specimen stage in an interior, windowless room. All images were processed the MacOS image editing software Paintbrush 2.1 to exclude any region of the image depicting the specimen stage. We then processed masked images through a custom python script that converted the image data to CIE LAB colorspace. Average luminance was then calculated from all unmasked pixels of the image, thus capturing information from the entire dorsal region of the study skin. Higher luminance values indicate lighter pelage color.

Environmental Data

Our study largely used environmental variables that have been shown to be important predictors in Gloger’s rule in other animal systems where studies have typically shown that animals are darker where environmental conditions are warmer, wetter, and darker (Delhey, 2019). The variables directly related to Gloger’s rule that were used in our study can be grouped into three categories: (1) temperature (mean annual temperature), (2) environmental moisture (mean annual precipitation and humidity), and (3) light conditions (insolation, canopy closure, and net primary productivity). Given that more direct background matching with surrounding substrates can also be an important selective pressure on pelage color in rodents (Caro, 2005; Vignieri et al., 2010), we also included a fourth category of variables associated with the background environment (soil carbon content, soil density, and fire frequency). Soil density was included because sandy soils, which are generally lightly colored, have higher bulk densities than fine silts and clays because they have larger, but fewer, pore spaces (Ayres et al., 1973). Soil carbon content was included because soils with higher carbon content have higher rates of decomposition and thereby produce darker substrates (Schulze et al., 1993). Fire frequency was included because higher frequency of fires was shown in fox squirrels (Kiltie, 1989) to be an important predictor of darker pelage.

The environmental data for the nine predictor variables (mean annual temperature, mean annual precipitation, humidity, insolation, canopy closure, net primary productivity, soil density, soil carbon, and fire frequency) for all 137 squirrel species were acquired from multiple databases. Both mean annual precipitation and mean annual temperature values were acquired across each species range from the PanTHERIA database, which provides point estimates at species level resolution for a range of environmental and life history characteristics of mammals (Jones et al., 2009) (PanTHERIA metadata available at http://esapubs.org/archive/ecol/E090/184/metadata.htm). Canopy closure raster data, which serves as a proxy measure for environmental light conditions, were acquired from high-resolution global maps of forest cover (Hansen et al., 2013). This dataset provides world-wide canopy closure data at 1 m by 1 m resolution. Due to memory constraints we reduced this dataset to a 5 m by 5 m resolution using the “aggregate” function in the raster package in R (R Core Team, 2013; Hijmans and van Etten, 2014). Both soil carbon content and soil density raster data at a 5-km resolution was acquired from the ISRIC Soil Grids web app (Batjes et al., 2019). Humidity (column water vapor amount), net primary productivity, fire frequency, and solar radiation were obtained from the NASA MODIS satellite hosted on NASA Earth Observations¹. These data represent time averaged values for each month for approximately the past 20 years. We downloaded these data at an 11 km by 11 km resolution and averaged the monthly values into a single raster file. All raster data was clipped with extant range maps of all 137 species downloaded from the IUCN (2019) using the raster and rgdal packages in R (R Core Team, 2013; Hijmans and van Etten, 2014; Bivand et al., 2019) to obtain species level averages.

Model Selection

Following the recommendations of Cooper et al. (2016), we used an MCMC based approach to fit both a Brownian motion and Ornstein Uhlenbeck model of evolution to our data. This was implemented by the “fitContinuousMCMC” function in the R package geiger (Pennell et al., 2014). This was run for 1 × 10⁶ generations with posterior likelihood values sampled every 100 generations.

Phylogenetic Generalized Least Squares

All data for our phylogenetic generalized least square (PGLS; Martins and Hansen, 1997) analyses were centered by subtracting the mean and scaled by dividing by the standard deviation. We then ran pairwise phylogenetic regressions using phylolm with the model specified as “OUfixedRoot” between each environmental variable and the phenotype with 1000 bootstrap replicates. Next, we ran 36 pairwise phylogenetic regressions using the same approach between each of the nine environmental variable and each of the eight other environmental variable to assess multicollinearity among the data. Variables were considered collinear if their correlation coefficient was greater than 0.7. Finally, we constructed 12 regression models representing all possible combinations of non-collinear variables (Table 1). As above, this was implemented with the function “phylolm” from the R package phylolm (Ho and Ane, 2014). Model weights were then calculated according to Symonds and Moussalli (2011). Partial R² values were calculated using “R2.resid” function of the R package rr2 (Ives and Li, 2018).

Table 1. The 12 regression models tested by PGLS.

Results

Phylogenetic Tree

Partition Finder selected the GTR + I + G substitution model for nearly all partitions, expecting the wobble positions of GHR and CMYC being better modeled by K80 + G and HKY + G, respectively. Our Bayesian phylogenetic inference of Sciuridae revealed 57 nodes with significant (>0.95) posterior probabilities (Figure 2).

Figure 2. Phylogeny representing 137 species within Sciuridae. Branch lengths are proportional to the expected number of changes per site. Dark circles represent nodes with >95% posterior probability.

Model Selection for Comparative Phylogenetic Analysis

Both runs of the MCMC converged as evidenced by the stationarity of the trace plots. The Ornstein Uhlenbeck model (Log likelihood −883) had a consistently higher likelihood than Brownian motion (Log likelihood −918) after burn-in and was therefore used for all comparative analyses.

Phylogenetic Generalized Least Squares

The univariate regression analysis revealed that seven out of nine environmental variables exhibit a significant correlation with the phenotype (Table 2 and Figure 3). Pelage luminance is negatively correlated with soil carbon content, humidity, primary productivity, canopy closure, temperature, and precipitation. Soil density is the only variable to have a significantly positive correlation with pelage luminance. Solar radiation and fire frequency exhibit no significant correlation with pelage luminance.

Table 2. Results of univariate regression analysis.

Figure 3. Plots showing the relationship between pelage luminance and each of the nine environmental variables included in this analysis. The solid line represents the line of best fit, while the dashed lines represent 95% confidence intervals from 1000 bootstrap replicates.

The regressions among environmental variables revealed that soil density is collinear with soil carbon content and insolation. Soil carbon content is also collinear with insolation. Water vapor is collinear with primary productivity, canopy closure, temperature and precipitation. Primary productivity is collinear with canopy closure and precipitation. Canopy closure is also collinear with precipitation.

The top three multiple regression models form the credible set, the set whose cumulative model weight exceeds 0.95 (Symonds and Moussalli, 2011) and we reject the remaining models. Model one finds significant negative effects between pelage luminance and fire frequency, primary productivity, and mean annual temperature as well as a significant positive correlation between the phenotype and solar radiation. Model three differs from model one only in the replacement of primary productivity with canopy closure, where the correlation is likewise significantly negative. Model two includes significant negative correlations between pelage luminance and soil carbon content and humidity, while the correlation with fire frequency is non-significant (Table 3). It is worth noting that mean annual precipitation, commonly found to be a significant correlate of pelage luminance, does not appear in the credible set.

Table 3. Coefficients of top three models from Table 1.

Discussion

Our comparative study of dorsal pelage luminance across Sciuridae based on multiple regression found support for the predictions of Gloger’s rule and other background matching expectations. The dorsal pelage of squirrel species become darker in environments with greater mean annual temperature, greater environmental moisture (humidity), lower light conditions (lower levels of solar radiation, greater primary productivity, and greater canopy closure), and darker background environments (greater soil carbon content and greater fire frequency). Our univariate regression analyses further supported these relationships. These findings are consistent with most predictions of Gloger’s rule that have been observed in mammals, including carnivores (da Silva et al., 2016; Caro et al., 2017), artiodactyls (Stoner et al., 2003b), Palearctic shrews (Stanchak and Santana, 2019), primates (Kamilar and Bradley, 2011; Santana et al., 2012), rodents (Lai et al., 2008), and marsupials (Nigenda-Morales et al., 2018; Cerezer et al., 2020). These empirical studies have commonly used temperature, precipitation, latitude, and vegetation to explain patterns of color variation that are associated with climatic factors or surrogates of those factors. Our study extended the general predictions of Gloger’s rule by also revealing that environmental variables related to background matching to surrounding substrates also contribute to pelage luminance.

Our results largely agreed with other quantitative investigations of Gloger’s rule in other Sciurids that were focused on intraspecific variation in Scirus niger (Kiltie, 1989) and Protoxerus stangeri (Amtmann, 1965). The former investigated several populations within S. niger and found support for darker squirrels being associated with greater temperature, moisture, and fire frequency. The mechanism proposed by Kiltie (1989) for the role of fire frequency is that wildfires leave behind a dark charred substrate, which may provide greater concealment for darker squirrels. However, they also stated that the lasting effect of this mechanism may be tenuous because of the rapid regrowth of vegetation following fires usually obscures the dark background of charred wood. Perhaps the mechanism is not fire frequency, but other environmental characteristics that are associated with habitats that burn frequently. More research is needed to determine if there are other probable mechanistic explanations for the association of fire frequency with pelage luminance and whether this trend is generalizable beyond Sciuridae. Amtmann (1965) analyzed samples from across the range of P. stangeri and similarly found a negative relationship between mean annual precipitation and pelage luminance. However, contrary to our study, P. stangeri was found to exhibit a positive correlation between pelage luminance and mean annual temperature.

Primates and squirrels present an interesting opportunity to study parallel evolution in animal color patterns. Primates, like Sciurids, possess exceptional diversity in color patterns (Bradley and Mundy, 2008; Santana et al., 2012) and share many other attributes, including wide-spread diurnality, many species exhibiting an arboreal lifestyle, having similar predators (carnivores, raptors, and snakes), and having global distributions that span a range of biomes. Our results are consistent with prior work on primates (Kamilar and Bradley, 2011) that found pelage luminance across primates corresponded negatively with AET. AET measures the amount of water leaving an environment and is affected by temperature, precipitation, humidity and primary productivity. This result implies that primates are darker in areas with higher precipitation, higher temperature, higher humidity and higher primary productivity. The one major difference between primates and Sciurids with regards to color pattern evolution is that many primate species have evolved striking sexual dimorphisms in pelage coloration that are related to social signaling and sexual condition (Setchell, 2005; Dubuc et al., 2009; Marty et al., 2009). Squirrels do not exhibit sexual dimorphism in pelage coloration and thus the selective pressures on pelage coloration are likely to be strictly due to ecological and environmental factors.

The range of color patterns and hues in the family Sciuridae are incredibly broad and diverse (Thorington et al., 2012). This diversity includes species with patterns that are patchy (e.g., Sciurus variegatoides), striped (e.g., Tamias striatus), spotted (e.g., Xerospermophilus spilosoma), mottled (e.g., Otospermophilus beecheyi), and uniform (e.g., Microsciurus mimulus), as well as species with colors that are very light (e.g., Callosciurus finlaysonii), very dark (e.g., Callosciurus prevostii), orangish-red (e.g., Sciurus vulgaris), reddish-maroon (e.g., Ratufa indica), reddish-brown (e.g., Tamiasciurus douglasii), and gray (e.g., Sciurus griseus). Furthermore, color patterns can vary tremendously just within a single animal (e.g., tail colors are often distinct from body colors). The evolutionary forces shaping this diversity of color patterns and hues across species and on different body parts are not well studied in Sciurids, but may include selection to enhance conspecific and heterospecific signaling, despite camouflage also likely being another important driver of pelage color evolution. Squirrels often vocalize and exhibit tail-flagging when they detect potential predators (Owings and Coss, 1977; McRae and Green, 2014). Thus, it is possible that striking color patterns on tails or even the main body may improve visual communication with predators during precarious encounters. We suspect that multiple types of selective pressures are shaping the evolution of fur color in squirrels and our study sheds some light that precipitation-related and background matching pressures may be shaping one dimension (pelage luminance) of these color patterns.

One of the general statistical challenges in studies of Gloger’s rule is that many of the hypothesized explanatory variables are inherently correlated with each other and may be directly or indirectly related to animal-color patterns that conform to Gloger’s rule. Ideally, studies can illuminate which variables are directly driving variation in pelage luminance and which variables may be mediated through others. Unfortunately, univariate or multiple regression approaches do not offer the ability to distinguish between direct and mediated effects when questions are ecologically complex. As an example, such approaches cannot distinguish between a model in which temperature and moisture drive primary productivity, which in turn drives pelage luminance; or one in which temperature and moisture drive primary productivity and pelage luminance simultaneously. One possible approach to differentiating between these two hypotheses is through statistical methods like path analysis, which has recently been extended for use with comparative phylogenetic data (von Hardenberg and Gonzalez-Voyer, 2013; van der Bijl, 2018). However, one limitation of path analysis is that it requires the hypotheses to be tested to explicitly specify all relationships among included variables. Given the complicated interconnectedness of variables typically included in Gloger’s rule studies, this specificity is difficult to attain without further methodological development. We recommend that future research explores how path analysis can be effectively applied to questions concerning Gloger’s rule.

Data Availability Statement

The datasets generated for this study can be found in the Dryad repository 10.5061/dryad.tqjq2bvv8.

Author Contributions

AS and AC contributed to the conception and design of the study, and draft and revision of the manuscript. AS was responsible for the data acquisition, analysis, and data interpretation. Both authors approved the publication.

Funding

AS was supported by the Mayer’s Undergraduate Research Summer Fellowship and the Undergraduate Research Scholarship through the Ohio State University. AC was supported by startup funds through the Ohio State University.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank D. Lunde and the Smithsonian National Museum of Natural History for access to specimens. We thank R. Zenil-Ferguson and L. Kubatko for their invaluable advice on phylogenetic comparative methods as well as K. Delhey and S. Nigenda-Morales for their helpful comments. We also thank T. Piko for help with data collection and figure formatting.

Footnotes

1. ^ https://neo.sci.gsfc.nasa.gov

References