https://orcid.org/0009-0009-0538-3673
Abstract
Contrary to expectations regarding efficient predator education mediated by lack of ambiguity and enhanced prey recognition, aposematic signals often show considerable intraspecific variability. For example, some striped skunks (Mephitis mephitis) are almost entirely white, others have black-and-white stripes of equivalent thicknesses, yet others are mostly black. We tested the ecological correlates of this variation in patterning using 749 museum skins collected across North America. Skunks had longer white–black borders and more bilaterally symmetrical stripes in areas with a greater number of potential predator species, and this effect was more marked for mammalian than avian predators, the latter of which may be less deterred by noxious defenses. Skunks from locations with greater predator diversity were less variable in the extent of whiteness on their dorsa and less variable in the length of their white–black borders, suggesting strong selection from predators leads to greater conformity in stripe patterns, even at the same location, but weak selection from predators leads to relaxed selection on pattern conformity. Skunks exhibited greater areas of black pelage in areas of greater humidity conforming to Gloger’s rule. Our results indicate that relaxed predation pressure is key to warning signal variation in this iconic species, whereas stronger pressure leads to signal conformity and stronger signals.
Introduction
Aposematic animals possess chemical, morphological, or behavioral antipredator defenses and signal their unprofitability to predators using distinctive odors, sounds, or appearances (Howell et al., 2021; Poulton, 1890; Rojas et al., 2015). Warning signals operate by eliciting avoidance from predators by capitalizing on innate fear, by facilitating learning, or often both (Endler & Mappes, 2004; Stevens & Ruxton, 2012). Visual warning signals are often rendered conspicuous because they are bright, markings exhibit internal contrast within the body, contrast against its background, or a combination of these factors (Stevens & Ruxton, 2012). Frequently, they consist of striking colors and repetitive patterns making bearers easily recognized and remembered thereby deterring predatory attack (Lindström et al., 1999; Roper & Redston, 1987; Rowland et al., 2010; Ruxton et al., 2018). Indeed, selection is thought to enforce conformity in the pattern of warning signals across individuals within species reducing ambiguity and promoting speedy learning by the predator (Ruxton et al., 2018; Sherratt & Christopher, 2003). Deviations from an average signal appearance or pattern should theoretically increase incidences of mistaken attacks by predators and hinder predator learning. Thus, selection is expected to favor uniform warning signals and suppress variation (Poulton, 1890), a hypothesis supported by several examples from the lab (Rowland et al., 2010; Sherratt & Christopher, 2003; Stevens, 2007) and field (Chouteau et al., 2016; Comeault & Noonan, 2011; Noonan & Comeault, 2009; Ogilvie et al., 2021). In short, the protective value of aposematic signals is enhanced with greater pattern element size and reduced with increasing bilateral asymmetry (Forsman & Merilaita, 1999), suggesting selection acts on signal variation to favor larger elements, more contrasting borders, and greater pattern symmetry. Despite apparent advantages, however, signal conformity is far from the norm in many aposematic species or populations (Briolat et al., 2018), yet the ecological drivers of such variation are only starting to be addressed. The main drivers of such intraspecific variation in warning coloration are thought to be differences in predators’ perception, variation in predator species’ tolerance to defenses, differences in predator learning ability, and differences in predator motivation (Briolat et al., 2018), but, above all, differences in predation pressure (Endler & Mappes, 2004). For example, Mochida (2011) showed that local variation in predation pressure explained differences in aposematism of island and mainland newts (Cynops pyrrhogaster) such that newts performed more aposematic behavior in areas with greater predation by mammals. Competing selective pressures from sexual selection, thermoregulation, and disease risk may be involved in warning signal variation too (Briolat et al., 2018).
In this study, we use an iconic aposematic signal, the black-and-white striping of striped skunks Mephitis mephitis (von Schreber, 1776) to explore selective forces that maintain aposematic signal variation. Mephitis mephitis’s white-on-black dual stripe signal, which advertises its noxious anal gland directed-spray defense, varies greatly between populations and within the same populations; stripe patterns even vary within litters. Furthermore, stripe symmetry within an individual varies (Figure 1), thus providing three sorts of natural manipulation for study. While other studies have found that internal contrasting stripes speed aposematic learning in birds (Aronsson & Gamberale-Stille, 2012), the limited empirical evidence that exists on predator interactions with skunks suggests that mammalian predators are more hesitant to approach and attack skunks with black-and-white stripes relative to no stripes (Hunter, 2009; Schiefelbein & Stankowich, 2016; Vo, 2021) and skunks with more white fur than those with less (Fay, 2017; Vo, 2021). Using 749 georeferenced museum specimens covering most of the species’ range (Figure 2), we examined the effects of different sorts of predation pressure on three aspects of signal conformity both between individuals (total area of whiteness, total length of contrasting white–black stripe borders) and within individuals (stripe bilateral symmetry). We additionally investigated the effects of environmental factors known to influence coat appearance in other mammals. Specifically, we examined four hypotheses as follows:
Striped skunk specimens mapped according to capture origin. Depicted here are museum M. mephitis specimens that have been selected to demonstrate variation in stripe signal that occurs between and at times within populations (Delaware), as well as a group in Chico, CA (bottom right) demonstrating signal conformity with low variance that might be expected of aposematic species.
Predator risk maps showing geographic variation in predation risk. GIS maps of calculated (A) mammalian carnivore risk scores and (B) avian predator risk scores according to predator species ranges, mapped together with sampled M. mephitis specimen locales (red circles). Note that a single red circle may represent many specimens collected from a single location, resulting in fewer circles than the total sample size. Darker brown areas indicate higher potential mammalian risk and darker purple areas indicate higher potential avian risk.
Hypothesis 1: Greater predation risk favors stronger, more consistent signals, while decreased predation risk due to fewer predator species relaxes selection on warning signal conformity. Striped skunks inhabiting areas with elevated potential predation risk from mammals and/or birds should exhibit stronger signals, which we considered to be (a) a greater area of white to contrast with the black body and dark substrate and longer contrasting white–black stripe borders; (b) greater uniformity in extent of whiteness and length of contrasting stripe borders across individuals (i.e., reduced intraspecific variation); and (c) greater stripe bilateral symmetry within individuals so as to rapidly educate and reinforce signal efficacy to predators. Conversely, striped skunks with reduced potential predatory threat should exhibit weaker signals with greater variation in their striping pattern, reflecting relaxation in selection for a strong signal (Briolat et al., 2018; Lahti et al., 2009).
Hypothesis 2: Pelage coloration is tailored to type of predatory threat: mammalian predators and avian predators may impose different forms of selection. Given that mammalian predators (carnivorans) have a strong sense of smell, but avian predators of nocturnal skunks (typically owls: Huey, 1931; Lesmeister et al., 2010; Marks et al., 1999; Rashid, 2015; Wilkinson, 1913) do not (Payne, 1971; Roper, 1999), aposematic signaling might be more targeted toward mammals than birds. Thus, selection from mammalian predators favoring increased area of whiteness, longer contrast borders, more highly bilaterally symmetrical pelages, which we assume maximizes conspicuousness, may be stronger than selection from avian predators. In short, the positive effects of mammalian potential predation risk on signals should be stronger than the effects of avian potential predation risk.
Hypothesis 3a: Aposematic signals must be enhanced in particularly dark and particularly light habitats. Lighting and background coloration of the natural habitat strongly impact the effectiveness of aposematic signals (Caro & Koneru, 2021). To this end, we expect striped skunks to have, on average, greater extent of dorsal white pelage in shadier habitats (e.g., poorly lit woodlands) so as to stand out from dark backgrounds and, conversely, greater extent of dorsal black pelage (reduced whiteness) in bright snowy habitats to contrast against the white background.
Hypothesis 3b: Aposematic species attempt to use background matching in some environments to avoid detection, but still rely on bold patterns when detected to deter attack. Recent studies, including one on spotted skunks (Spilogale gracilis, Merriam, 1890), have identified some aposematic species as possessing signals that allow for crypsis at a distance but appear bold and conspicuous up close (Barnett et al., 2016; Caro, 2013; Sherratt & Christopher, 2003; Tullberg et al., 2005). If M. mephitis stripe signals perform similarly, we expect their pelages to be black-and-white for up-close encounters, but with a pattern that enables them to blend into their habitats when viewed from a distance. We therefore predict a greater extent of dorsal dark pelage (a reduced amount of white in stripes) in striped skunks from less-lighted habitats such as woodlands, and a greater extent of white dorsal striping in individuals from well-lit and colder habitats (e.g., white snowy). The former prediction was upheld in white-backed hog-nosed skunks (Conepatus leuconotus, H. Lichtenstein, 1832), where dorsa were darker in more forested habitats (Ferguson et al., 2022).
Hypothesis 4: Pelage coloration varies in accordance with humidity. In accordance with Gloger’s rule (Delhey, 2019; Gloger, 1833), striped skunks from more humid areas should have darker pelages (less white). This was the result found in white-backed hog-nosed skunks where dorsa were darker in more humid (and perhaps forested) habitats (Ferguson et al., 2022).
Methods
Stripe data collection
The central component of this study was a library of high-quality images of striped skunk dorsa generated from preserved stuffed study skins. A total of 749 specimen photographs were gathered from skin holdings at the Smithsonian'sNational Museum of Natural History (USNM), American Museum of Natural History (AMNH), Royal Ontario Museum (ROM), Los Angeles County Museum of Natural History (LACM), University of California Berkeley Museum of Vertebrate Zoology (MVZ), California Academy of Sciences (CAS), and Oklahoma State University Museum (OSUM) by HW and TS (see Supplementary Data Set). The use of museum specimens allowed for an extremely large sample size that was tailored to represent all 13 recognized subspecies of striped skunk and the species’ entire North American geographic range (Hall, 1981). We argue that skins are a representative sample of living skunks in an area because they are normally trapped using bait by museum personnel, rarely taken from road kills, and are not trophy species that may be overrepresented by one sex or larger sizes in museum collections.
At each museum, photographs of selected specimens were captured using a whiteness-controlled Nikon DSLR camera. Prior to photographing, tag identification number, geographic site of collection, and date of collection were recorded for each specimen. Specimens were then combed in a cranial–caudal direction to smooth hairs flat and laid out on gridded (2.5 cm × 2.5 cm) paper next to a color standard card. Photographs capturing dorsal views of each specimen (Figure 1) were taken from a camera that had been mounted on a tripod and secured 1.5 m from the floor, with its lens pointed directly down and exactly perpendicular to the study skin. Any skins that appeared significantly twisted, overly faded, missing large patches of fur covering stripes, or in otherwise poor condition were excluded from the study.
ImageJ (Schneider et al., 2012) was used to digitally trace the perimeter of each specimen’s set of dorsal stripes. Photos were calibrated using the 2.5 cm2 gridded background paper for reference. Body length for the specimen pictured was estimated using the line segment feature, starting at the nose tip and ending at the tail base. The left white stripe was then traced using the freehand selection tool, followed by the right white stripe. Stripes joined at the head were bisected at the body midline. The Measure tool was then used to calculate the area (square millimeter) and perimeter length (millimeter) of each stripe.
Following Ferguson et al (2022), we calculated a unitless Whiteness Index (WI) to capture the area of white fur on the dorsum for each specimen:
We also developed a new unitless measure of Contrast Index (CI) for each specimen to capture the total length of the black–white border:
Finally, we created a measure of stripe bilateral Asymmetry (Asym) to capture symmetrical appearance:
which results in a value of 0 for perfect bilateral symmetry and 0.5 for complete asymmetry.
These three measures were calculated for each skin and used as the dependent measures of signal strength. For WI and CI, body length was used as a standard rather than body area because preserved study skin girths are often distorted with preparation techniques and manners of storage. WI and CI were positively correlated (Pearson’s r = .63, Supplementary Figure S5). At moderate and high values of WI, CI plateaued due to stripes typically extending the full length of the dorsum.
Specimen mapping
All museum specimens were georeferenced and mapped in ArcGIS 10.3 mapping software according to their capture locales. Because a significant proportion of capture locations listed on museum specimen tags were only specific to counties, estimates for each specimen were calculated at the county level.
Predation risk
Evidence of predation rates on skunks is scarce, but focused studies of predator diets suggest that skunks are uncommon to exceedingly rare in mammalian predator diets (usually less than 5% occurrence in scats); although striped skunk fur does occasionally appear in scats from coyotes, Canis latrans Say, 1823 (Reed et al., 2023), mountain lions, Puma concolor Linnaeus, 1771 (Cassaigne et al., 2016), jaguars, Panthera onca Linnaeus, 1758 (Cassaigne et al., 2016), and bobcats, Lynx rufus von Schreber, 1777 (Tewes et al., 2002). Skunks are, however, somewhat regularly preyed upon by great horned owls (Bubo virginianus, Gmelin, 1788), and very rarely by other raptors (del Hoyo et al., 1994, 1999; Platt & Rainwater, 2012); although most evidence of predation on skunks by mammals and birds comes from predator stomach and scat contents, conflating predation and scavenging. We adapted techniques established in Stankowich et al. (2014) to geographically estimate potential predation risk throughout North America for comparison to museum M. mephitis specimen capture locales. Working under assumptions that for another species to pose a threat to a striped skunk, that predator must live in the same geographic range, (a) live in the same type of habitat, (b) have a diet that includes mammals, (c) be large enough in size to be able to capture and kill a striped skunk, and (d) be active at the same time of day, we estimated potential risk posed by each North American predator species to striped skunks by multiplying together numerical estimates for each of these factors. Briefly (full details are available in Supplementary Material), we first overlaid georeferenced M. mephitis specimens onto IUCN range maps for North American mammalian carnivore species and bird of prey (Strigiformes, Acciptriformes, and Falconiformes) species using ArcGIS software. Then for each specimen, we identified all carnivoran and avian species overlapping in range with its point of capture. For each predator species with geographic overlap, we calculated its potential risk to a striped skunk by multiplying together numerical estimates of the remaining assumptions: (a) share at least one habitat type in common (1/0), (b) activity overlap (i.e., the degree of nocturnality of the predator: 0–1), (c) proportion of the predator’s diet that includes mammals (0–1), and (d) probability that the predator can kill a prey the size of a skunk based on the size of the largest prey species that predator is known to eat (see Supplementary Material for full details on how (a)–(d) were calculated). For each striped skunk potential predator species, we calculated the potential risk score (a × b × c × d; range 0-1) posed to striped skunks (see Supplementary Material for the full risk data set). For each M. mephitis specimen capture locale, we then summed the predation risk scores for all local mammalian predators into a potential mammalian carnivore risk score (Figure 2A) and for all local avian predators into a potential avian risk score (Figure 2B). Potential mammalian and avian risk scores were moderately positively correlated (Pearson’s r = .645; Supplementary Table S2). Although this method is agnostic to predator abundance, we think that it is appropriate as predator numbers have changed considerably since the time that many of the skunk specimens were collected and, in most areas of North America, are unknown (i.e., local predator abundance data are simply not available for any species across the range).
Environmental factors
Using ArcGIS, museum specimens were georeferenced on detailed habitat type, mean winter temperature, average winter snow depth, and average humidity maps.
Habitat type was estimated using a published map of North American land cover from the North American Environmental Atlas (Commission for Environmental Cooperation, 2005) (Supplementary Figure S1). The map provided a harmonized view of the physical cover of Earth’s surface across the continent based on 2005 Moderate Resolution Imaging Spectroradiometer (MODIS) satellite imagery at a spatial resolution of 50 m2. Nineteen Level II land cover classes were defined using the Land Cover Classification System (LCCS). Of these possible classes, striped skunk range included all except for water (class 18). To address the camouflage/background matching abilities of striped skunks, we followed the methods of Caro et al. (2021) and utilized ~5 photos of each land cover class to estimate shadiness scores for each class on a scale of 1–10 with 1 having the most ambient light reaching the ground (brightest) and 10 being the least amount of ambient light reaching the ground (shadiest). Each land cover class was given a 1–10 score of shadiness by HW and TS, working independently. The two independent scores were averaged for each land cover class to create a single score of Shadiness for all classes.
North American winter temperature and snow depth were also mapped. We obtained a winter temperature data set from WorldClim (Hijmans et al., 2005). From WorldClimVersion 2.0’s Bioclimactic Variables, we selected mean temperature of coldest quarter of the year (Bio11) to represent winter temperature. All available years (1970–2000) were averaged together in ArcGIS to compile a single map (Supplementary Figure S2). In addition, we utilized snow depth data available as part of ERA-Interim, a global atmospheric reanalysis produced by the European Centre for Medium-range Weather Forecasts (Dee et al., 2011). We compiled a single map of average winter (December–February) snow depth (centimeter) for years 1901–1950 in ArcGIS at a resolution of 0.75 km2 (Supplementary Figure S3).
We compiled a single map of average noontime relative humidity for years 1979–2000 (temporal availability limited due to satellite basis of this particular set of data) at a resolution of 0.75 km2 from ECMWF’s CERA-20 reanalysis (Supplementary Figure S4).
For each of the above environmental datasets (shadiness, temperature, snow depth, and humidity), we used the georeferenced location of each specimen to assign it a local value at the county level for each measure to use in analysis.
Analysis
We examined which predation and environmental factors influenced the whiteness index and contrast index using linear models that allow the expected value and the standard deviation to vary by observation; both the expected value and the standard deviation were modeled as a function of covariates (Posthuma Partners, 2019). Note that we did not pool individuals into local populations to estimate standard deviation. The standard deviation model was included to account for strong heteroskedasticity and because signal conformity is best measured as the variation in whiteness and contrast indices. As an example of our model structure, whiteness index was modeled using the following equations:
where βμ and βσ are parameter vectors and Xμ and Xσ are covariate matrices. The standard deviation was modeled using a log link to prevent negative estimates. We modeled asymmetry index using a generalized linear model with a gamma distribution and log link. The gamma distribution was appropriate for the positive, continuous asymmetry index values, and the log link prevented negative values.
Covariates were included to test the hypotheses described above. Avian and mammalian predation risk were included as covariates in all three models (i.e., whiteness index, contrast index, and asymmetry index) for the expected value (main effect of each factor) and were additionally included as covariates for the standard deviation in WI and CI. We additionally included habitat shadiness, snow depth, mean temperature of the coldest quarter of the year, and humidity for the expected values in the whiteness index model, and we included habitat shadiness and snow depth for the expected values in the contrast index model. None of the covariates were strongly correlated (all Pearson’s r < .7; Supplementary Table S1). We z-score standardized all covariates to allow for comparisons among effect sizes . For models including a log link (the standard deviation of WI and CI, and the expected value of the asymmetry index), covariates were estimated on the log scale. To help with interpretation, we report the standardized effect size as one minus the exponential of these covariate estimates times 100, which is the percent increase in the dependent variable for each 1 SD increase in the covariate. For models using an identity link (the expected value of WI and CI), the standardized effect size is the raw increase in the dependent variable for each 1 SD increase in the covariate.
All models were analyzed in a Bayesian framework in the program JAGs (Plummer, 2003), called from the programming language R (R Core Team, 2012) using the rjags and jagsUI packages (Kellner, 2019; Plummer, 2018). We used 3 chains with a burn-in of 10,000 iterations, ran for an additional 10,000 iterations and thinned by 5. Priors were all drawn from a normal distribution with a variance of 10,000. Convergence was assessed using the Gelman and Rubin potential scale reduction statistic Ȓ (Gelman & Rubin, 1992) and visual inspection of the plotted chains and posteriors. The Ȓ values were ≤1.001 suggesting the model converged well. We report the probability of direction (pd) in the results, which is the probability that the parameter differs from the null expectation (usually 0) in the more common direction of the parameter estimates; pd values greater than .95 are considered to have moderate support or be statistically significant (similar to α = 0.05) and pd > .999 are considered to have very strong support.
Results
H1 and H2: predation risk influences
We found no mean effects of either mammalian or avian predation risk on WI (Figure 3A and B; Table 1, pd = .85 and .74, respectively), but pelages had relatively longer contrasting borders (CI) in areas with greater mammalian risk (Figure 3D; Table 2; pd > .999, effect size = 0.084 increase per SD) and greater avian risk (Figure 3E; Table 2; pd > .999, effect size = 0.168 increase per SD). Importantly, the standard deviations (i.e., variations) of both WI and CI were significantly higher in areas with reduced mammalian and avian risks (Figure 3C and F; Tables 1 and 2; pd > .999 for all four estimates), suggesting that relaxation in predation risk leads to increased variation in pattern. The standard deviation in CI decreased by 13.4% and 34.4% per 1 SD increase in avian and mammalian predation risk, respectively, and the standard deviation in WI decreased by 12.4% and 18.5% per 1 SD increase in avian and mammalian predation risk, respectively.
Parameter estimates for the whiteness index GLM. Variance estimates are on the log scale. Percentages are the 95% credible interval of parameter estimates, and pd is the probability that the parameter is in the more common direction (similar to a p-value). A pd ≥ .95 can be considered significant or moderate support, and a pd > .999 can be considered strong support. Note that in the text effect sizes for the variance model are reported as the percent changes, which was calculated as one minus the exponential of the mean estimate times 100.
| Main effects | Mean | 2.50% | 97.50% | pd |
|---|---|---|---|---|
| Intercept | 0.298 | 0.292 | 0.303 | >.999 |
| Avian predation risk | −0.002 | −0.008 | 0.004 | .742 |
| Mammalian predation risk | 0.004 | −0.003 | 0.011 | .845 |
| Shadiness | 0.007 | 0.002 | 0.013 | .993 |
| Temperature | 0.003 | −0.003 | 0.008 | .82 |
| Humidity | −0.019 | −0.025 | −0.013 | >.999 |
| Snow depth | −0.004 | −0.009 | 0.000 | .955 |
| Variance model | ||||
| Intercept | −2.643 | −2.692 | −2.594 | >.999 |
| Avian predation risk | −0.132 | −0.201 | −0.073 | >.999 |
| Mammalian predation risk | −0.204 | −0.276 | −0.138 | >.999 |
| Main effects | Mean | 2.50% | 97.50% | pd |
|---|---|---|---|---|
| Intercept | 0.298 | 0.292 | 0.303 | >.999 |
| Avian predation risk | −0.002 | −0.008 | 0.004 | .742 |
| Mammalian predation risk | 0.004 | −0.003 | 0.011 | .845 |
| Shadiness | 0.007 | 0.002 | 0.013 | .993 |
| Temperature | 0.003 | −0.003 | 0.008 | .82 |
| Humidity | −0.019 | −0.025 | −0.013 | >.999 |
| Snow depth | −0.004 | −0.009 | 0.000 | .955 |
| Variance model | ||||
| Intercept | −2.643 | −2.692 | −2.594 | >.999 |
| Avian predation risk | −0.132 | −0.201 | −0.073 | >.999 |
| Mammalian predation risk | −0.204 | −0.276 | −0.138 | >.999 |
Parameter estimates for the whiteness index GLM. Variance estimates are on the log scale. Percentages are the 95% credible interval of parameter estimates, and pd is the probability that the parameter is in the more common direction (similar to a p-value). A pd ≥ .95 can be considered significant or moderate support, and a pd > .999 can be considered strong support. Note that in the text effect sizes for the variance model are reported as the percent changes, which was calculated as one minus the exponential of the mean estimate times 100.
| Main effects | Mean | 2.50% | 97.50% | pd |
|---|---|---|---|---|
| Intercept | 0.298 | 0.292 | 0.303 | >.999 |
| Avian predation risk | −0.002 | −0.008 | 0.004 | .742 |
| Mammalian predation risk | 0.004 | −0.003 | 0.011 | .845 |
| Shadiness | 0.007 | 0.002 | 0.013 | .993 |
| Temperature | 0.003 | −0.003 | 0.008 | .82 |
| Humidity | −0.019 | −0.025 | −0.013 | >.999 |
| Snow depth | −0.004 | −0.009 | 0.000 | .955 |
| Variance model | ||||
| Intercept | −2.643 | −2.692 | −2.594 | >.999 |
| Avian predation risk | −0.132 | −0.201 | −0.073 | >.999 |
| Mammalian predation risk | −0.204 | −0.276 | −0.138 | >.999 |
| Main effects | Mean | 2.50% | 97.50% | pd |
|---|---|---|---|---|
| Intercept | 0.298 | 0.292 | 0.303 | >.999 |
| Avian predation risk | −0.002 | −0.008 | 0.004 | .742 |
| Mammalian predation risk | 0.004 | −0.003 | 0.011 | .845 |
| Shadiness | 0.007 | 0.002 | 0.013 | .993 |
| Temperature | 0.003 | −0.003 | 0.008 | .82 |
| Humidity | −0.019 | −0.025 | −0.013 | >.999 |
| Snow depth | −0.004 | −0.009 | 0.000 | .955 |
| Variance model | ||||
| Intercept | −2.643 | −2.692 | −2.594 | >.999 |
| Avian predation risk | −0.132 | −0.201 | −0.073 | >.999 |
| Mammalian predation risk | −0.204 | −0.276 | −0.138 | >.999 |
Parameter estimates for the contrast index GLM. Variance estimates are on the log scale. Percentages are the 95% credible interval of parameter estimates, and pd is the probability that the parameter is in the more common direction (similar to a p-value). A pd >= 0.95 can be considered significant or moderate support, a pd > 0.999 can be considered very strong support. Note that in the text effect sizes for the variance model are reported as the percent changes, which was calculated as one minus the exponential of the mean estimate times 100.
| Main effects | Mean | 2.50% | 97.50% | pd |
|---|---|---|---|---|
| Intercept | 3.221 | 3.174 | 3.27 | >.999 |
| Avian predation risk | 0.084 | 0.047 | 0.123 | >.999 |
| Mammalian predation risk | 0.168 | 0.124 | 0.213 | >.999 |
| Shadiness | −0.036 | −0.075 | 0.002 | .966 |
| Snow depth | −0.029 | −0.052 | −0.005 | .994 |
| Variance model | ||||
| Intercept | −0.594 | −0.645 | −0.543 | >.999 |
| Avian predation risk | −0.144 | −0.216 | −0.074 | >.999 |
| Mammalian predation risk | −0.423 | −0.489 | −0.357 | >.999 |
| Main effects | Mean | 2.50% | 97.50% | pd |
|---|---|---|---|---|
| Intercept | 3.221 | 3.174 | 3.27 | >.999 |
| Avian predation risk | 0.084 | 0.047 | 0.123 | >.999 |
| Mammalian predation risk | 0.168 | 0.124 | 0.213 | >.999 |
| Shadiness | −0.036 | −0.075 | 0.002 | .966 |
| Snow depth | −0.029 | −0.052 | −0.005 | .994 |
| Variance model | ||||
| Intercept | −0.594 | −0.645 | −0.543 | >.999 |
| Avian predation risk | −0.144 | −0.216 | −0.074 | >.999 |
| Mammalian predation risk | −0.423 | −0.489 | −0.357 | >.999 |
Parameter estimates for the contrast index GLM. Variance estimates are on the log scale. Percentages are the 95% credible interval of parameter estimates, and pd is the probability that the parameter is in the more common direction (similar to a p-value). A pd >= 0.95 can be considered significant or moderate support, a pd > 0.999 can be considered very strong support. Note that in the text effect sizes for the variance model are reported as the percent changes, which was calculated as one minus the exponential of the mean estimate times 100.
| Main effects | Mean | 2.50% | 97.50% | pd |
|---|---|---|---|---|
| Intercept | 3.221 | 3.174 | 3.27 | >.999 |
| Avian predation risk | 0.084 | 0.047 | 0.123 | >.999 |
| Mammalian predation risk | 0.168 | 0.124 | 0.213 | >.999 |
| Shadiness | −0.036 | −0.075 | 0.002 | .966 |
| Snow depth | −0.029 | −0.052 | −0.005 | .994 |
| Variance model | ||||
| Intercept | −0.594 | −0.645 | −0.543 | >.999 |
| Avian predation risk | −0.144 | −0.216 | −0.074 | >.999 |
| Mammalian predation risk | −0.423 | −0.489 | −0.357 | >.999 |
| Main effects | Mean | 2.50% | 97.50% | pd |
|---|---|---|---|---|
| Intercept | 3.221 | 3.174 | 3.27 | >.999 |
| Avian predation risk | 0.084 | 0.047 | 0.123 | >.999 |
| Mammalian predation risk | 0.168 | 0.124 | 0.213 | >.999 |
| Shadiness | −0.036 | −0.075 | 0.002 | .966 |
| Snow depth | −0.029 | −0.052 | −0.005 | .994 |
| Variance model | ||||
| Intercept | −0.594 | −0.645 | −0.543 | >.999 |
| Avian predation risk | −0.144 | −0.216 | −0.074 | >.999 |
| Mammalian predation risk | −0.423 | −0.489 | −0.357 | >.999 |
Scatterplots indicate the effects of (A) mammalian (brown) and (B) avian (blue) risk on Whiteness Index (WI); how changes in these risk measures influence variation in WI (C); effects of (D) mammalian and (E) avian risk on Contrast Index (CI); how changes in these risk measures influence variation in CI (F); and effects of (G) mammalian and (H) avian risk on Asym. Predation Risk scores in (C) and (F) are standardized to allow for easier comparison between avian and mammalian predator risks. The confidence bandsare the 95% credible intervals.
Asymmetry declined significantly and strongly with increased mammalian predation risk (Figure 3G; Table 3; pd > .999) and avian predation risk (Figure 3H; Table 3; pd > .999). Specifically, asymmetry declined by 14.5% and 31.1% per 1 SD increase avian and mammalian predation risk, respectively.
Parameter estimates for the asymmetry GLM. Estimates are on the log scale. Percentages are the 95% credible interval of parameter estimates, and pd is the probability that the parameter is in the more common direction (similar to a p-value). A pd >= .95 can be considered significant or moderate support, a pd > 0.999 can be considered strong support. Note that in the text effect sizes for the variance model are reported as the percent changes, which was calculated as one minus the exponential of the mean estimate times 100.
| Mean | 2.50% | 97.50% | pd | |
|---|---|---|---|---|
| Intercept | −4.312 | −4.389 | −4.236 | >.999 |
| Avian predation risk | −0.157 | −0.242 | −0.058 | >.999 |
| Mammalian predation risk | −0.372 | −0.465 | −0.286 | >.999 |
| Mean | 2.50% | 97.50% | pd | |
|---|---|---|---|---|
| Intercept | −4.312 | −4.389 | −4.236 | >.999 |
| Avian predation risk | −0.157 | −0.242 | −0.058 | >.999 |
| Mammalian predation risk | −0.372 | −0.465 | −0.286 | >.999 |
Parameter estimates for the asymmetry GLM. Estimates are on the log scale. Percentages are the 95% credible interval of parameter estimates, and pd is the probability that the parameter is in the more common direction (similar to a p-value). A pd >= .95 can be considered significant or moderate support, a pd > 0.999 can be considered strong support. Note that in the text effect sizes for the variance model are reported as the percent changes, which was calculated as one minus the exponential of the mean estimate times 100.
| Mean | 2.50% | 97.50% | pd | |
|---|---|---|---|---|
| Intercept | −4.312 | −4.389 | −4.236 | >.999 |
| Avian predation risk | −0.157 | −0.242 | −0.058 | >.999 |
| Mammalian predation risk | −0.372 | −0.465 | −0.286 | >.999 |
| Mean | 2.50% | 97.50% | pd | |
|---|---|---|---|---|
| Intercept | −4.312 | −4.389 | −4.236 | >.999 |
| Avian predation risk | −0.157 | −0.242 | −0.058 | >.999 |
| Mammalian predation risk | −0.372 | −0.465 | −0.286 | >.999 |
In general, the effects of potential mammalian predation risk on WI, CI, and Asym were stronger than the effects of avian risk. While there were no overall effects of predation risk on the expected whiteness index, the decrease in variation with stronger mammalian predation risk on WI was larger than the similar decrease in variation with stronger avian risk (Figure 3C; Table 1; pd = .897). Similarly, mammals had a stronger effect than avian predators by 100% and 194% for the expected value and standard deviation of CI, respectively (Figure 3D–F: Table 2; mean effect: pd = .960; SD effect: pd = .998). Finally, the effect of greater mammalian risk on increased stripe symmetry was 136% stronger than that of avian risk (Figure 3G and H; Table 3; pd = .993).
H3a, H3b, and H4: environmental influences
We found that skunks had a greater amount of white on their coats (Figure 4A; Table 1; pd = .993, effect size = 0.007 increase per SD) and shorter contrasting stripe borders (Figure 4C; Table 2; pd = .966, effect size = 0.036 decrease per SD) in shadier habitats. In snowy habitats, they had a smaller area of white fur (Figure 4B; Table 1; pd = .955, effect size = 0.004 decrease per SD) and shorter contrasting stripe borders (Figure 4D; Table 2; pd = .994; effect size = 0.029 decrease per SD); all of these effects were relatively weak but statistically significant. There was no effect of cold winter temperatures on WI (Table 1, pd = .82). These results provide support for enhanced signaling in dark and light habitats (H3a) and allow us to reject a background matching hypothesis (H3b).
Scatterplots indicate the effects of increased shadiness and greater snow depth on Whiteness Index (WI) (A: shade; B: snow depth) and Contrast Index (CI) (C: shade; D: snow depth). The confidence bands are the 95% credible intervals.
Humidity had a strong effect on WI, where skunks in more humid areas had less white on their dorsa (Figure 5; Table 1; pd > .999), supporting Gloger’s rule (H4). Specifically, WI decreased by 0.019 per 1 SD increase in humidity, which was over two times stronger than the effect of any other covariate on the expected WI.
Scatterplot indicating the effect of increased environmental humidity on the whiteness index. The confidence band is the 95% credible interval.
Discussion
General
Results of this study indicate that variation in striped skunk pelage coloration and patterning is principally driven by the magnitude of potential predation risk. Specifically, enhanced potential risk from mammalian carnivore predators appears to be a strong selective driver of pattern boldness (a greater extent of black–white contrast and more symmetrical patterns when risk is high) and reduced variation in these aspects of stripe signal. Conversely, specimens originating from areas with low predation, particularly from mammals, tended to be more variable in extent of whiteness and lengths of pelage contrast borders, with greater stripe asymmetry. Of the environmental factors tested, relative humidity appeared to be the strongest selective force, with M. mephitis following other mammalian species, including other skunk species (Ferguson et al., 2022), in adherence to Gloger’s rule (Delhey, 2019): specimens clearly displayed darker pelages in areas of high humidity.
Stripe signals and predation
We hypothesized that M. mephitis stripe signal variation reflects decreases in predatory risk and predicted that areas of high predation pressure would have skunks with low signal variation (thus signal conformity). Analyses confirmed that the perimeter length of white–black fur (together comprising a bold signal) increased slightly with increasing mammalian and avian predation risk and that variation in area of white and length of contrast border (signal conformity) decreased markedly with increasing predation risk, with the latter effects being much stronger than the former. Thus, higher potential predation risk favors bold signals and intraspecific uniformity. On top of this, we found that within-individual signal conformity (symmetricity of dorsal stripes to one another) increased with increasing predation risk—further suggesting the importance of aposematic signal conformity for predator education in high-risk areas and relaxation in selection in between-individual and within-individual signal repetition or conformity in areas with fewer predators. Briolat et al. (2018) reviewed mechanisms that might maintain aposematic signal variation, and our evidence supports their view that when predation risk is reduced locally, purifying selection is relaxed, allowing more extreme phenotypes (i.e., all black bodies with a white head marking through to complete white capes where the white stripes merge; Figure 1) to survive and reproduce, gradually leading to more variable signals. Conversely, in areas with high densities of predators, purifying selection should be strong, leading to more uniform signals (i.e., black bodies with moderately-thick white stripes running from head to tail to maximize contrast borders).
Although aposematic patterning is known to vary according to differences in predation and predation risk, it has been primarily studied in insects and amphibians where variant morphs are more apparent. For example, poison dart frogs (Oophaga granulifera, Taylor, 1958 and O. pumilio, Schmidt, 1857), include aposematic morphs, green cryptic morphs, and morphs that are intermediate; in this case, frogs from populations that experience high attack rates have more conspicuous coloration and behave more boldly (Briolat et al., 2018; Maan & Cummings, 2012; Willink et al., 2013). We expected to see that reduced variability in stripe pattern would be associated with greater mammal risk than avian risk, as mammals are likely more deterred by spray defenses and so payoffs for educating them are greater (Fisher & Stankowich, 2018; Wade-Smith & Verts, 1982). Analyses did indeed show that the amount of signal contrast increased significantly with mammal risk more than avian risk, although both showed positive relationships. Again, increased variation in amount of white and length of contrast borders was seen with both decreased mammalian and avian risk, but more so with mammalian risk. These results support previous empirical work (Mann, 2018) showing that a striped skunk population with more coyotes and mountain lions had longer and less variable stripes and less variable concentrations of noxious thiols than a skunk population lacking these predators. We speculate that although stronger signals may make skunks more conspicuous to avian predators (which are less deterred by their defensive sprays), longitudinal stripes might cause avian predators to misjudge their strikes (Murali & Kodandaramaiah, 2016; Murali et al., 2018; Scott-Samuel et al., 2023). Arguably, this is supported by our finding that greater potential avian risk was associated with greater CI (Figure 3E). Thus, mammalian risk and avian risk appear to both be drivers of patterning in striped skunks, but results indicate that the impacts of mammalian risk are stronger. Other studies have shown that different antipredator coloration strategies may be tailored to different predator types. For example, a study of aposematic Australian brood frogs (Pseudophryne spp.) suggests that dorsal signaling is tailored to deter avian predators that view them from above, while dorsal and ventral signals are equally moderately effective at deterring mammalian predators, which hunt through chemosensory cues and rooting behavior when approached from the side (Lawrence et al., 2018).
Our results showed that skunks exhibited relatively more white pelage but less contrasting patterns in shady habitats and relatively more dark pelage and more contrasting patterns in bright snowy areas, which is the opposite of what was found in white-backed hog-nosed skunks (Ferguson et al., 2022). That study suggested that whiter pelage in bright, exposed environments was meant to advertise more strongly to predators, and darker pelage in forests was a form of background matching, resulting in a selective balance between crypsis and aposematism. Our results suggest that patterns are favored to be more conspicuous in both habitat types. In snow, darker (less white) patterns likely result in narrower, longer stripes (more contrast border) to boost black contrast against the white snow for both predator types. The possibility that these differences might be related to different anal spray components (Wood, 1990; Wood et al., 1993) and severity of their effects on predators should not be overlooked, although Mann (2018) did not find a correlation between stripe length and concentration of noxious thiols in the anal spray.
Stripe signals and environmental factors
We also explored the possibility that striped skunk aposematic coloration may conform to Gloger’s rule, which predicts darker coloration in areas of higher humidity, and expected M. mephitis specimens from humid locations to exhibit darker pelages (Gloger, 1833). We found an increase in the extent of dark pelage with higher humidity, which is consistent with Gloger’s rule and the findings by Ferguson et al. (2022) in white-backed hog-nosed skunks. However, strong selection for the extent of dark pelage should reduce the variance in coloration, but contrary to this expectation, variation increased.
Gloger’s rule has been documented in a variety of birds and mammals (Burtt & Ichida, 2004; Ferguson et al., 2022; Kamilar & Bradley, 2011; Nigenda-Morales et al., 2018; O’Neill et al., 2017), but the evolutionary drivers remain obscure. The melanization of feathers in birds and pelage in mammals responsible for making homeotherms darker may result in a reduction in UV radiation reaching the skin, reduced abrasion, protection against bacterial infections in the tropics, or blending in with shady environments caused by luxuriant vegetation. Our findings cannot distinguish between these alternatives, but the fact that striped skunks retain their white stripes even at low latitudes suggests that aposematism is a more important evolutionary driver of pelage coloration in this species than Gloger’s rule.
Conclusions
This is the first intensive study of aposematism in the iconic striped skunk and the first to show the importance of relaxed selection on aposematic variation in mammals. Greater predation risk, mainly by mammals, favors stronger, more symmetrical, and less variable stripe patterns, with signals becoming whiter in darker habitats and darker in whiter habitats, tuned to enhance the signal in different lighting environments. We also found further support for Gloger’s rule in this species. Future studies could focus on how predators learn about skunk appearances, specifically on whether contrast intensity and pattern type are shaped by different selective forces (Vo, 2021) and whether skunks living in areas with relaxed selection from predators behave differently than those in areas with greater potential predation risk.
Data availability
The data underlying this article are available in Supplementary Material.
Author contributions
H.W. collected stripe data and environmental data, helped conduct statistical analysis, and wrote an early full draft of the manuscript. T.C. helped conceptualize the estimates of predation risk, advised on the scope of the analysis, and helped write and edit the manuscript. D.B. helped run the statistical analysis and helped edit the manuscript. A.F. advised on the creation of stripe measures and helped edit the manuscript. T.S. conceptualized the original project, collected stripe data, created estimates of predation risk, advised on the scope of the analysis, and helped write and edit the manuscript.
Funding
Research was supported in part by California State University Long Beach.
Conflict of interest
The authors declare no conflict of interest.
Acknowledgments
We wish to thank the museum staff at the Smithsonian Museum of Natural History (USNM), American Museum of Natural History (AMNH), Royal Ontario Museum (ROM), Los Angeles County Museum of Natural History (LACM), University of California Berkeley Museum of Vertebrate Zoology (MVZ), California Academy of Sciences (CAS), and Oklahoma State University Museum (OSUM) for facilitating access to their collections. For assistance with the geospatial analyses used in calculating potential predation risk, we thank J. McNamara, C. Fay, M. Kresky, P. Haverkamp and the Center for Spatial Technologies and Remote Sensing (CSTARS) at UC Davis, and J. Blossom and the Center for Geographic Analysis at Harvard University. We also thank D. Johnson for statistical advice and A. Carter, M. Cummings, and two anonymous reviewers for comments on a previous version of this manuscript.
