Ecological Drivers of Eastern Fox Squirrel Pelage Polymorphism

The color patterns of an animal’s pelage, feather, or skin serve a variety of adaptive functions; importantly, one function is concealment through background matching. In spatially and temporally heterogeneous environments, some species exhibit multiple distinct color patterns within a population (i.e., color polymorphism). The environmental drivers of color polymorphism are poorly understood. We used the polymorphic eastern fox squirrel (Sciurus niger ssp.; hereafter, fox squirrel) as a model species to investigate the role of environmental factors on pelage coloration. Building upon previous research that investigated the drivers of melanism, we measured fox squirrel pelage coloration across the visible light spectrum. Agouti-colored squirrels were positively associated with increased proportion of burned area in a fox squirrel dispersal buffer. Light-colored (less melanistic) squirrels were positively associated with increasing proportion of cropland in a fox squirrel dispersal buffer. We posit that agouti pelage is broadly adapted to a range of heterogeneous conditions created by fire. Conversely, croplands, once established, are relatively stable ecosystems which promote a consistently adaptive light-colored pelage morph. We suggest that in an increasingly human-dominated environment, spatially and temporally homogeneous processes, such as prescribed burning, may not sufficiently recreate environmental heterogeneity, which could result in lost pelage diversity.


INTRODUCTION
One of the most distinct forms of animal diversity is the variety of color patterns animals exhibit. The color patterns of an animal's pelage, feather, or skin serve a variety of adaptive functions such as concealment, communication, and thermoregulation (Caro, 2005;Ancillotto and Mori, 2017;Cuthill et al., 2017). Animal coloration patterns are an evolutionary response to selective forces that promote fitness-maximizing color patterns (Cloudsley-Thompson, 1999;Caro, 2005). In prey species, pelage coloration can reduce predation risk, often through concealment, thereby reducing an animal's likelihood of being identified by a predator (Krupa and Geluso, 2000;Boratyñski et al., 2014). In heterogeneous or changing environments, a prey species with a consistent color pattern may be maladapted (Kettlewell, 1955;Mather, 1955). Accordingly, color polymorphism is often associated with spatial and temporal environmental heterogeneity that promotes multiple strong selective forces on coloration (Ford, 1945;Mather, 1955;Merilaita et al., 1999;Bond and Kamil, 2006;Gray and McKinnon, 2007). Understanding the specific environmental factors that influence pelage coloration provides insight into the adaptations of polymorphic species and the ecological factors that promote the persistence of polymorphism.
Eastern fox squirrels (Sciurus niger ssp.; hereafter, fox squirrel) found in the southeastern United States are a model species for polymorphism research in heterogeneous environments (i.e., fire-prone pine [Pinus spp.] forests). The fox squirrel is the most variably colored mammal in North America (Hall, 1981), such that individuals are identifiable using their unique pelage patterns (Tye et al., 2015). In the southeastern United States, fox squirrels typically have gray or agouti dorsal pelage, with varying amounts of melanism (Moore, 1956). Agouti pelage has multiple colors on individual hair follicles, and often appears as a mottled brown color. Fox squirrels show greater variation in ventral coloration, which can be cream, red, beige, or black (Moore, 1956). Broadly, fox squirrels in Florida can be classified into one of six color morphs: gray, buff/tan, agouti with partial dorsal melanism, black-bellied, dorsal melanism, and melanistic (Figure 1; Florida Fish and Wildlife Conservation Commission, 2013). Fox squirrels commonly exhibit countershading, a pelage pattern where an animal is dorsally dark and ventrally light (Steele and Koprowski, 2001). Countershading may conceal animals from above and below, especially in well-lit environments (Caro, 2005), though evidence for countershading as an adaptive coloration is mixed (Kiltie, 1988;Ruxton et al., 2004;Yahner, 2012). Countershading is thought to be a more effective antipredator adaptation in marine environments, where predators swim above and below prey (Kiltie, 1988). Similarly, arboreal animals, such as fox squirrels, face predation threats from above and below (Koprowski, 1994), and may therefore benefit from countershading.
Previous research on fox squirrel color polymorphism in the southeastern United States suggests that pale-agouti dorsal pelage, defined as agouti pelage with no non-agouti black pelage, is most cryptic when static against a background of bark from multiple species of tree (Kiltie, 1992b) and in unburned FIGURE 1 | Exemplars of the six fox squirrel color morphs found in Florida. From left to right, gray, agouti with partial dorsal melanism, buff/tan, black-bellied, dorsal melanism, melanistic. environments (Kiltie, 1992a). Dark and intermediate dorsal coloration is more cryptic against burned backgrounds, and moving squirrels with partial dark dorsal pelage may be more difficult for avian predators to identify (Kiltie, 1992b). Melanistic fox squirrels are most concealed against fire-blackened substrates, but trees, other vegetation, and bare soil only remain blackened after a fire for approximately 2 weeks (Kiltie, 1992b). Two weeks is likely too brief for concealment to be the adaptive benefit that allows melanistic fox squirrels to persist (Kiltie, 1992b).
Research on the ecological drivers of fox squirrel coloration show a positive correlation between dorsal melanism and historic lightning-caused wildfire frequency and temperature, and a negative correlation with elevation (Kiltie, 1989). Kiltie (1992b) suggested that dark coloration persists because dark squirrels are more difficult for aerial predators to identify, not because dark pelage is more cryptic in a fire-dominated environment. Historic fire patterns are only one possible ecological driver of squirrel polymorphism. Other ecological factors such as landcover (Kettlewell, 1955;Sumasgutner et al., 2018), precipitation/moisture (da Silva et al., 2016, and soil type (Krupa and Geluso, 2000) drive polymorphism in other species, and may therefore shape fox squirrel polymorphism. Furthermore, previous work on fox squirrels focused on the ecological drivers of dorsal melanism (Kiltie, 1989(Kiltie, , 1992a, while explaining variation in fox squirrel pelage across the full color spectrum is likely to provide a more complete understanding of the evolutionary links between coloration and the environment. Color vision is important for hunting in some species of fox squirrel predators, including canids (Kasparson et al., 2013) and raptors (Håstad et al., 2005), suggesting that pelage darkness alone may not be enough to describe an animal's cryptic pelage adaptation.
The aim of this study was to identify environmental variables that drive fox squirrel pelage polymorphism. We hypothesized that fox squirrel pelage variation would correlate with the background color of the environment as an antipredator adaptation. We predicted that squirrels with agouti-colored pelage would positively correlate with fire frequency and burned area (Kiltie, 1989). Dark-colored pelage in rodents, including sciurids, is often associated with low-light environments, including closed-canopy forests and increased tree density (Wauters et al., 2004;Ancillotto and Mori, 2017). Thus, we predicted that squirrels with more melanistic pelage would positively correlate with closed-canopy forests where decreased light penetration through the canopy makes dark pelage more cryptic (Wauters et al., 2004;Ancillotto and Mori, 2017). Lastly, we predicted that countershading would negatively correlate with canopy cover as countershading is an antipredator adaptation in well-lit environments (Caro, 2005).

MATERIALS AND METHODS
Between 2012 and 2015, we collected roadkill fox squirrel specimens across Florida (Figure 2). To increase sample sizes and coverage throughout Florida, we enlisted the help of various wildlife agencies and citizen scientists to collect roadkill fox squirrel specimens. Although our sampling design was opportunistic, we managed to collect fox squirrels throughout most of their range within Florida. We recorded the location of each roadkill specimen using a handheld GPS. To prevent the pelage color from fading, we kept all specimens frozen until they were prepared as flat skins and dried. We stored all specimens in a dark cabinet at the Florida Museum of Natural History (FLMNH).
We used a digital color sensor (Nix TM Pro Color Sensor, Nix Sensor Ltd., Hamilton, ON, Canada) to quantify squirrel pelage color composition (Hodgen, 2016;Stiglitz et al., 2016;Holman et al., 2018). In a separate study, we found that the sensor gives highly correlated measures of fox squirrel pelage compared to traditional photography methods (Potash et al., 2020). The sensor measures the average color within a 15 mm diameter round aperture. We described and recorded our color readings using the Commission Internationale de L'Eclairages (CIELab) color space which uses three axes, L * , a * , and b * to define colors (Schanda, 2007). The L * axis measures brightness and ranges from 0 to 100, with higher values indicating a lighter color (i.e., L * = 0 is black and L * = 100 is white). Values on the a * and b * axes can be positive or negative, where at a * = 0 and b * = 0 corresponds to the color gray. The a * axis describes a continuous gradient where negative values indicate more of the color green, while positive values indicate more of the color red. The b * axis describes a continuous gradient where negative values indicate more of the color blue, while positive values indicate more of the color yellow.
We measured L * , a * , and b * values at five locations on each squirrel specimen. These locations were between the ears (head), centered between the front legs (shoulder), at the base of the tail (rump), midway between the shoulder and rump (dorsum), and laterally perpendicular to the dorsum along the edge of the specimen (venter; Figure 3). To increase the accuracy of our measurements, we averaged 7 repeated measurements at each location, removing and replacing the sensor at each location between each recording (Holman et al., 2018).

Environmental Data
We measured environmental characteristics that previous studies show influence pelage variation within a 3.7 km radius around each specimen . The 3.7 km radius (hereafter, buffer) is the average dispersal distance of southeastern fox squirrels (Wooding, 1997). We measured seven environmental variables, the proportion of burned area between 1992-2015 (Prop_Burned), the number of fires that occurred between 1992-2015 (Count_Fire), the mean percent canopy cover (Canopy), the mean elevation in meters (Mean_Elevation), the mean annual precipitation from 1981-2014 in cm (Precip), the dominant soil order (Soil), and proportion of landcover type.
We chose proportion of landcover as a metric to describe the available habitat to each squirrel, and to test how much variation in pelage could be ascribed to habitat. We calculated the proportion of landcover type in each squirrel buffer using the 2011 National Land Cover Database (NLCD; Homer et al., 2015). To reduce the number of landcover types, we binned the original classifications into 19 classes based on the Florida Natural Areas Inventory (Florida Natural Areas Inventory, 2010; Tye et al., 2017). To further reduce the number of landcover types in our analysis, we only included landcover classes where fox squirrels have been observed in Florida, which were Hardwood and Hammock, Pinelands, Crops, Coniferous Wetlands, Parks Cemeteries and Golf Courses, and Tree Plantations .
We selected Prop_Burned and Count_Fire as explanatory variables because previous studies have shown correlations between fox squirrel dorsal melanism and fire-blackened substrates (Kiltie, 1992b) as well as the number of historic wildfires (Kiltie, 1989). We calculated Prop_Burned and Count_Fire in each buffer using the spatial wildfire occurrence dataset from 1992-2015 (Short, 2017). The spatial wildfire occurrence dataset consists of ignition points for wildfires, each of which has an associated attribute for the final fire size. We created circular polygons around each ignition point such that the area of the circle equaled the fire size of the associated point and measured the amount of burned overlap in each squirrel buffer. In each buffer, we measured a burned area only once even if the same area burned on multiple occasions.
We included Canopy as an explanatory variable to describe the openness of the habitat, which influences dorsal pelage in rodents (Ancillotto and Mori, 2017). We calculated mean percent canopy in each squirrel buffer from the 2011 NLCD Tree Canopy Cartographic dataset (Homer et al., 2015). Each 30 m 2 pixel in this raster layer contains a single value for percent canopy cover, which we averaged for each squirrel buffer. We included Precip and Elevation as explanatory variables in our model because both variables have been correlated with fox squirrel pelage variation in previous studies (Kiltie, 1989). We used a digital elevation model (DEM; Gesch et al., 2018) with 30 m 2 resolution to calculate Mean_Elevation within a squirrel buffer. To calculate Precip, we averaged the annual precipitation from 1981-2010 (PRISM Climate Group, 2020) within a squirrel buffer. We retrieved soil data from the Web Soil Survey (USDA Natural Resources Conservation Service, 2020), and calculated the proportion of all soil orders in a fox squirrel dispersal buffer. To reduce the number of parameters in our model, we included only the most common (greatest proportion in a buffer) soil order for each squirrel in our analysis.

Statistical Analysis
To better understand the environmental variables that influence fox squirrel polymorphism, we subset the pelage color data into four non-exclusive groups. The first subset (Specimen, Full Color [SFC]) contained the L * , a * , and b * color data from all the locations we measured on each specimen. The second subset (Specimen, Brightness [SB]) included only the L * axis from each of the measured body components. This allowed us to explicitly test our prediction that darker pelage would positively correlate with increasing canopy cover in a fox squirrel's dispersal buffer. The third (Dorsum, Full Color [DFC]) and fourth (Dorsum, Brightness [DB]) data subsets also included the L * , a * , and b * , and the L * measurements only, respectively, but only used the three measurements from dorsal locations (shoulder, dorsum, rump). These subsets allowed us to compare our results to past research that focused only on fox squirrel dorsal pelage variation (Kiltie, 1989(Kiltie, , 1992a. To identify the environmental drivers of pelage variation, we conducted a redundancy analysis (RDA) for each subset of squirrel pelage color values. An RDA extends multiple regression to test relationships between a multivariate response matrix and multiple predictor variables. For analysis the response matrix was a squirrel-by-color measurement matrix for the particular group (i.e., SFC SB, DFC, or DB) and the response variables were the seven environmental variables we detailed above (van den Wollenberg, 1977). We standardized all continuous response and predictor variables to have a mean of 0 and standard deviation (SD) of 1. For each RDA, we used a forward selection approach to identify the model terms that create the most parsimonious model of predictor variables (Blanchet et al., 2008). Forward selection iteratively models individual predictor variables to the response variables, and additively combines significant predictor variables until additional predictors do not significantly improve model fit (R 2 ). We used α ≤ 0.05 for determining significant variables. Using both α and R 2 as stopping criteria in forward selection reduces type I error, a common problem with forward selection approaches (Blanchet et al., 2008). If forward selection identified more than one significant predictor variable, we used variance partitioning to quantify the unique amount of variation explained by each term (Peres-Neto et al., 2006).
For all RDAs we used a permutation test with 999 permutations to test for significance of the model terms and canonical axes (Legendre et al., 2011). We calculated an adjusted R 2 (R 2 adj) using Ezekiel's formula (Legendre et al., 2011) to determine the overall amount of variation explained by each model. To ensure model results were not influenced by spatial autocorrelation (Tobler, 1970) we plotted residuals from each model as a function of spatial distance between specimens (Wagner, 2004). We found that model residuals showed no dependency on distance and model results were not affected by spatial autocorrelation (Wagner, 2004). We therefore did not include any spatial covariates in our model. We performed all multivariate statistical analysis in R Version 3.5.0 (R Core Team, 2018) using the package Vegan (Oksanen et al., 2017).

Pelage Values
We measured pelage coloration from 135 fox squirrel specimens. Fox squirrels showed the greatest variation along the L * axis (brightness), followed by the b * axis (blue to yellow coloration), and the least amount of variation along the a * axis (green to red coloration). Fox squirrels exhibited the least amount of variation overall in head coloration as indicated by the narrowest range and lowest SD along the L * (mean = 18.8 ± 3.10 SD), a * (mean = 0.50 ± 0.47 SD), and b * (mean = 1.00 ± 1.00 SD) axes (Figure 4). Fox squirrel venters had the widest range and largest SD for all color axes (L * : mean = 62.90 ± 11.40 SD; a * : mean = 5.40 ± 2.69 SD; b * : mean = 16.5 ± 5.57 SD). The three dorsal metrics, shoulder, dorsum, and rump were similar along the a * and b * axis (a * Shoulder : mean = 2.17 ± 1.09 SD; b * Shoulder : mean = 7.72 ± 3.20 SD; a * Dorsum : mean = 1.94 ± 1.06 SD; b * Dorsum : mean = 7.55 ± 3.41 SD; a * Rump : mean = 2.13 ± 0.98 SD; b * Rump : mean = 8.66 ± 2.84 SD). The L * value was similar for shoulder and dorsum (L * Shoulder : mean = 35.30 ± 8.34 SD; L * Dorsum : mean = 36.20 ± 9.58 SD), but was greater and had a lower SD for rump (L * Rump : mean = 40.40 ± 7.59 SD). Greater L * values for the rump compared to the dorsum and shoulder indicate that fox squirrel dorsal coloration lightens (decreased melanism) between the central portion of the dorsum and the base of the tail. Qualitatively, increased a * and b * dorsal values represent agouti-colored pelage, while lower L * values represent melanistic pelage (Figure 4).

RDAs
Using the SFC data subset, the forward selection identified the proportion of burned area (Prop_Burned) as the only significant model term (p < 0.001) in the RDA. Since there was only one predictor variable, there was only one RDA axis that explained pelage variation (p = 0.003). The model adjusted R 2 = 0.03. In general, a greater proportioned of burned area positively correlated with greater values of a * (increased red) and b * (increased yellow), and negatively correlated with increasing L * (lighter colors) for the venter (Figure 5A).
Using the SB data subset, the forward selection identified an additive model with the predictor variables Crops (p = 0.01), Precip (p = 0.03), and Canopy (p = 0.014) as the most parsimonious model. Using variance partitioning we found that Crops, Precip, and Canopy accounted for 3.8, 1.9, and 1.9% of the variation in pelage brightness, respectively. Only the first RDA axis explained a significant amount of variation in the pelage data (RDA1, p = 0.002). Crops loaded most on RDA1 (0.67) followed by Precip (−0.46) and Canopy (0.15). Precip loaded most heavily (0.81) on RDA2, but RDA2 did not explain a significant amount of variation in the data (p = 0.31). The model adjusted R 2 = 0.05. Along RDA1, increasing proportion of row crops was positively associated with increasing L * values for all body components except head (Figure 5B).
The forward selection for the DFC subset identified only Prop_Burned as a significant model term (p < 0.001). Since there was only one predictor variable, there was only one RDA axis that explained pelage variation (p < 0.001) with model adjusted R 2 = 0.05. Increasing Prop_Burned was most associated with increased rump a * coloration, as well as increased rump b * and dorsum a * (Figure 5C).
The forward selection for the DB subset identified only Crops as a significant model term (p = 0.016). Since there was only one predictor variable, there was only one RDA axis that explained pelage variation (p = 0.009) with model adjusted R 2 = 0.03. Increasing proportion of row crops was associated with increasing L * in shoulder, dorsum, and rump coloration ( Figure 5D).

DISCUSSION
We found that variation in fox squirrel pelage coloration was associated with differences in spatial and temporal environmental heterogeneity. Fires can increase heterogeneity in the environment at a variety of scales (McKenzie et al., 1995), while cropland remains comparatively stable. Fire may therefore promote pelage variation while crops promote a more uniform light-colored pelage morph. Our results support previous findings that environmental heterogeneity shapes pelage polymorphism by manifesting a variety of different colored backgrounds for prey to match (Kettlewell, 1955;Merilaita et al., 1999;Bond and Kamil, 2006).
We found support for our prediction that agouti pelage was positively associated with fire (Figures 3, 5). Previous research found that dorsally pale-agouti and intermediateagouti (partial melanism) colored fox squirrels were superiorly camouflaged against a variety of backgrounds in fire-prone ecosystems (Kiltie, 1992a,b). Compared to other fox squirrel pelage morphs, agouti pelage is the best overall at background matching against the ground and common trees in the southeast (Kiltie, 1992b). The large amount of variation within the agouti morph may make it well-adapted for concealment in the spatially and temporally heterogeneous habitats found in fire-driven ecosystems (McKenzie et al., 1995). For example, in the open pine forests of the southeastern United States, frequent fire can create a bare understory (Darracq et al., 2016), in which pale agouti rodents are most cryptic (Ancillotto and Mori, 2017). Fox squirrels in open pine forests rely on interspersed hardwoods for food and nest sites (Perkins et al., 2008), and hardwoods can create fire shadows which prevent complete burning (Jack and Pecot, 2017), increasing heterogeneity at a fine spatial scale. In contrast with burned areas, croplands in the Southeast are temporally and spatially stable, likely favoring a more consistent, uniformly colored population. We found that cropland was associated with decreased melanism in fox squirrels, possibly indicating that fox squirrels in human-engineered landscapes converge toward a uniform light-colored morph.
Dark-colored rodents are more cryptic in darker environments (Wauters et al., 2004;Ancillotto and Mori, 2017). However, we found no support for our prediction, that dark-colored pelage would be positively associated with closed canopy ecosystems. Fox squirrels are primarily found in open canopy forests (Moore, 1957;Weigl et al., 1989). Although fox squirrels use closed canopy hammocks occasionally (Weigl et al., 1989;Conner et al., 1999), periodic usage of closed canopy areas may not be enough to influence pelage coloration. Thus, fox squirrels may not use poorly lit closed canopy forests enough for dark pelage to have a cryptic advantage in these areas.
We did not find that countershading was associated with open canopy forest. Although forward selection identified canopy cover as a significant term in the SB model, canopy cover explained relatively little variation along the only significant RDA axis (RDA1). Furthermore, increased canopy cover was associated with lighter colored pelage in fox squirrel dorsal and ventral coloration (Figure 4). Countershading may not be an effective antipredator adaptation in terrestrial environments (Kiltie, 1988;Ruxton et al., 2004) and thus, countershading is unlikely to convey an antipredator advantage to fox squirrels.
Our analysis was limited by measuring pelage coloration in CIELab color space, which is based on human vision. Most animals have different visual acuity and color perception than humans (Endler, 1978;Cuthill et al., 2017), and therefore may not respond to the colors as we measured them; however, it was beyond the scope of this study to incorporate predator vision and behavior into our analysis. Our results explained only a small amount of pelage variation, suggesting that additional environmental factors not included in our analysis are likely driving pelage variation. Animal coloration is driven by additional factors other than predation risk (Smith et al., 1972), such as communication and thermoregulation (Stoner et al., 2003;Caro, 2005), hair strength (Bonser, 1995), or a combination of multiple factors (Cloudsley-Thompson, 1999). Fox squirrels are capable of frequent long-distance movements (Potash et al., 2018) and may therefore encounter a variety of ecological conditions that influence pelage coloration. Thus, there are likely many factors driving pelage coloration, preventing any individual variable we included in our analysis from explaining large amounts of pelage variation.
In this study we showed clear links exist between environmental heterogeneity and color polymorphism. , the x-axis shows the only significant RDA axis from forward selection, while the y-axis shows the first principal components (PCA). The PCA axis represents the residuals of the constrained ordination.
To preserve diverse animal color patterns, management and conservation efforts should consider the spatial and temporal scale of all management actions. Pelage adaptations to heterogeneity may be ineffective and lost in regularly managed ecosystems. Managing ecological processes, such as prescribed fire to promote heterogeneity rather than uniformity (Darracq et al., 2016) may be key for maintaining animal color diversity.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because the animals used in this study were previously collected opportunistically as roadkill and were part of the Florida Museum of Natural History's mammal collections.

AUTHOR CONTRIBUTIONS
AP, DG, LC, and RM conceived the original idea. AP, DG, and VM collected the data and developed data collection protocols. AP, BB, and RM developed the statistical analyses. All authors contributed to the final manuscript.