Playing by the rules? Phenotypic adaptation to temperate environments in an American marsupial

Museum specimens

We examined 352 (163 females and 189 males) study skins of Virginia opossum museum specimens to measure body and extremity dimensions, and 345 (159 females and 186 males) specimens for skin pigmentation and coat coloration analysis (Table S1). The specimens belong to five different natural history collections and were collected over 145 years (1865–2010; Table S2) along the geographic range of the species, from southern Nicaragua to the northern east coast of the United States (US), including some specimens from the western US in California (Fig. 1). We obtained the coordinates for each specimen collecting locality using two different approaches: (1) the coordinates provided by the original collectors; or (2) we located the collecting site provided by the original collector using the Google Earth program and extracted the geographic coordinates. Table S2, provides the geographic coordinates, year and country in which the specimens were collected. To avoid the confounding effects of age, all specimens measured for this study were adults. We determined the age following Gardner & Sunquist’s (2003) age estimation method based on tooth eruption. For each specimen, we took measurements to identify variation in body size, body extremities size, proportion of pigmented skin of the extremities and, body and facial coat coloration. Twelve phenotypic traits were measured in total, and because some specimens had damaged body parts, not all specimens were measured for all traits (see Table S1 ). The variation in some of the traits is shown in Fig. 2. This study was done following UCLA’s Office of Animal Research Oversight protocols #2011-121-02 and #2011-121-03.

Phenotypic measurements

To quantify variation in body size, we measured body length (head and body; Fig. 2A) and left hindfoot length (heel to tip of middle digit). Hindfoot size is often used as an alternative to represent body size because these traits frequently co-vary (Suttie & Mitchell, 1983; Martin et al., 2013). Consequently, since these measurements were positively correlated with each other (r = 0.547; P < 0.001), we used hindfoot length as an additional body size measurement in subsequent analysis, instead of as a measurement of extremity size. For appendage size, we measured the tail length (base to tip; Fig. 2B) and posterior aspect of the ear length (from base to most superior part of the pinna). We measured the posterior aspect of the ear because most specimens had their ears folded and it was impossible to measure the anterior part of the pinnae (the more commonly used measurement). For simplicity, hereafter, we will refer to the posterior aspect of ear length as ear length. All measurements were taken in centimeters to the nearest 0.1 decimal.

We measured four skin pigmentation traits by recording the proportion of the tail, ear, ventral and dorsal aspect of the hindfoot’s middle digit that was visibly pigmented (Figs. 2B, 2C). We measured the proportion of pigmentation on both the ventral and dorsal part of the middle digit because populations in the US have hindfeet with light pigmentation on the sole and ventral part of the toes, and in the most northern populations, the light pigmentation reaches the dorsal part of the toes (Gardner, 1973).

Body and facial coat coloration measurements were taken using the tristimulus colorimeter Minolta ChromaMeter CR-200 (Minolta, Osaka, Japan). This colorimeter measures the reflectance of a xenon flash light and records the color in the three-dimensional color space CIELab. The lightness axis (L*) expresses color brightness, with values ranging between 100 and 0 for white and black surfaces, respectively, whereas the a* and b* axes are the color coordinates (Fullerton & Keiding, 1997; Clarys et al., 2000). Since the phenotypes of coat coloration in the Virginia opossum range from dark to light (Fig. 2A) we only used the information of the lightness axis (L*) in our analysis. The probe with the light source was held at 90° angle to the surface parallel to the body axis, and the reflectance readings were recorded manually. All reflectance measurements in each site were measured in triplicate and averaged to calculate the total lightness value per site. We recorded reflectance values of lightness from three facial traits (rostrum, cheek and temporal region) and from the torso of the body in nine sites (i.e., three on the right flank, three on the dorsal stripe, and three on the left flank; see Fig. 3). The mean of the torso measurements was used as the average torso lightness value.

Environmental data

For each locality where specimens were collected, we obtained information at 1-km spatial resolution on 12 environmental variables, which have been reported to affect body and extremity size, and pigmentation in mammal species (Lindsay, 1987; Ashton, Tracy & De Queiroz, 2000; Chaplin, 2004; Caro, 2005; Cordero & Epps, 2012; Kamilar et al., 2012; Terada, Tatsuzawa & Takashi, 2012). Eight bioclimatic variables related to temperature (annual mean temperature, mean temperature of warmest quarter, mean temperature of coldest quarter), precipitation (annual precipitation, precipitation of coldest quarter), temperature seasonality (mean diurnal range, temperature seasonality) and precipitation seasonality were obtained from the WorldClim database (Hijmans et al., 2005; Table S3). These variables maximize variation in North America while minimizing correlation (Harrigan et al., 2014). In addition, four variables were extracted from remote-sensing databases: the maximum value of the normalized difference vegetation index (NDVIMAX) that is related to vegetation density and productivity (Tucker & Sellers, 1986; Buermann et al., 2008); the vegetation continuous field product as an estimate of the percentage of tree coverage (TREECOV; Hansen et al., 2002); the monthly backscatter measurements that capture attributes related to surface moisture (ROUGH) (Long et al., 2001); and elevation (ELEV) (Table S3 ). Elevation was included in the analyses due to its relationship with temperature variation, which can influence the variance of several traits we examined (Blackburn & Ruggiero, 2001; Symonds & Tattersall, 2010; Kamilar & Bradley, 2011).

Data analysis

We carried out Shapiro–Wilk normality and Bartlett homoscedasticity tests for all the traits and analyzed them accordingly. There is important sexual dimorphism in body and appendages size in opossums (Gardner & Sunquist, 2003); however, these measurements were normally distributed and a preliminary analysis of our data showed that they have similar patterns of latitudinal variation for males and females (Figs. S1A–S1D). To insure adequate sample size for body and extremity size measurements, we controlled for the effect of sex on these traits, and used the residuals to analyze both sexes together. To test if there were differences in skin pigmentation and coat coloration traits between males and females we used analysis of variance (ANOVA) or Wilcoxon rank-sum test depending on the normality of the data. All traits had homogeneity of variance, but only cheek and torso lightness were normally distributed and were the only coloration traits statistically different between males and females (cheek lightness: F = 29.85, P < 0.001; coat lightness: F = 8.85, P = 0.003), therefore we analyzed them separately for both sexes. All tests were done using their standard functions in the R statistical framework version 3.1.3 (R Core Team, 2015).

Tests of correlation with latitude

To determine the relationship between the phenotypic traits and latitude, we carried out Pearson (r) or Spearman (r_s) correlation analyses depending on the normality of the data in the R statistical framework (R Core Team, 2015). Additionally, using the ggplot2 package in R (Wickham, 2009), we plotted the trait values against latitude, implementing the non-linear regression loess function to graphically show the trend of correlations and the latitudes for which there is a change in this trend.

Association with environmental variables

To detect if spatial autocorrelation was present in our phenotypic data, we calculated Moran’s I weighted by the Euclidean distance between two points for each trait using the ape package v. 4.1 (Paradis, Claude & Strimmer, 2004) in R. To identify the best environmental models for predicting the phenotypic variation in the Virginia opossum, we ran random forest analyses using the randomForest package v. 4.6-12 in R (Liaw & Wiener, 2002; Prasad, Iverson & Liaw, 2006). Considering measurements from museum specimens that were collected in the same geographic location could contribute to increase spatial autocorrelation in our data; therefore, we averaged the phenotypic values of those specimens (Dataset S1) for the random forest analyses. To run the random forest analyses, we used the values of each phenotypic trait as response variables, and the 12 environmental variables and two geographic variables (i.e., latitude and longitude; Table S3) as predictors. The two geographic variables were included because incorporating geographic information in random forest models allows an evaluation of how much variation in response is explained by environmental variation as compared to simply geographic proximity (Evans et al., 2011; Mascaro et al., 2014). Decision trees (regression or classification) and random forest methods have no a priori assumptions about the relationship between predictor and response variables, allowing for the possibility to analyze non-linear relationships with complex interactions (Breiman, 2001; Cutler et al., 2007; Strobl, Malley & Tutz, 2009; Evans et al., 2011). Random forests analyses are iterations of large number of decision trees, which recursively partition the data into binary homogeneous groups splitting the response variable by the predictor variable explaining most of the remaining variance. The amount of variation in the response variable explained by each predictor is incorporated in the model. Applying a randomized bootstrapping (bagging) method, random forest analysis uses a subset of both response and predictor variables randomly permuted to construct each regression tree and assess the robustness of the model based on the remaining data not included in the tree. If the accuracy of the model decreases appreciably when a variable is left out of the model, that variable is considered an important predictor of the data (Breiman, 2001; Prasad, Iverson & Liaw, 2006; Strobl, Malley & Tutz, 2009).

We ran 10,000 regression trees for each random forest run, and after each run we eliminated the least important variables in the model and re-ran random forest until we identified the most predictive, least complex models, which were composed of three to five predictors that explained the largest amount of variation for each trait. We compared these combined models (which included environmental and geographic variables) with models composed only by the two geographic variables and models including environmental variables only. This was done to detect the role that geography or environment alone play in explaining phenotypic variation in the opossum. Finally, we also tested for spatial autocorrelation in the residuals of the combined models using the Moran’s I statistic to be confident that spatial autocorrelation has not affected the results of those models.

Latitudinal variation

Body and extremities size

We found correlation with latitude for all body dimension measurements (see Table 1; Figs. 4A–4D). The correlations with latitude for body size and hindfoot length were positive (r = 0.314, P < 0.001 and r = 0.284, P < 0.001, respectively), increasing after 27°N (Figs. 4A, 4B). In contrast, the correlations for tail length and ear length were negative (see Table 1). Tail length, was smaller above 27°N, whereas for ear length smaller values occur only in latitudes above 37°N (Figs. 4C, 4D). In general, opossums with smaller bodies but larger tails were distributed in latitudes below 27°N, and larger individuals with shorter tails were found in higher latitudes. This pattern of variation follows both Bergmann’s and Allen’s rule for body and extremity size, respectively.

Table 1:

Pearson and Spearman correlation between Didelphis virginiana phenotypic traits and latitude.

Trait	Sample size	Pearson correlation	Spearman correlation
Body length	348	0.314**	n.a.
Hindfoot length	345	0.284**	n.a.
Tail length	348	−0.613**	n.a.
Ear length	185	−0.172*	n.a.
Tail pigmentation	340	n.a.	−0.701***
Ear pigmentation	293	n.a.	−0.559***
Toe ventral pigmentation	334	n.a.	−0.583***
Toe dorsal pigmentation	338	n.a.	−0.548***
Rostrum lightness	345	n.a.	0.679***
Temporal region lightness	345	n.a.	0.719***
Cheek lightness
F	159	0.517***	n.a.
M	186	0.512***	n.a.
Torso lightness
F	159	0.298**	n.a.
M	186	0.251**	n.a.

DOI: 10.7717/peerj.4512/table-1

Notes:

F

Females

M

Males

n.a

Not applicable

*P < 0.05.

**P < 0.001.

***P < 0.0001.

Association with environmental variation

Concordant with our results of latitudinal variation, we found that all phenotypic traits show positive spatial autocorrelation (Table S4) and the most important variable for the models of these traits was latitude, whereas longitude was not consistently among the top variables (Fig. 5). Despite this association between morphology and latitude, environment also helped predict morphological characteristics; models with environmental or geographic variables alone explained less of the Virginia opossum phenotypic variation than the combined models including environmental and geographic variables together (Table 2). Further, the residuals of these combined models showed small negative Moran’s I values, most of them were slightly significant or not significant (Table S5), which indicates the models are appropriate for our analyses since they control for spatial autocorrelation. As such, we explore and explain these more predictive, combined models below.

Skin and coat pigmentation

Combined models explained 31.17–62.56% of the variation in skin pigmentation, 13.75–78.34% in face coloration, and 14.41–19.88% in torso coloration (Table 2). The most important and consistent environmental predictors for explaining phenotypic variation in skin pigmentation and face coloration traits were temperature seasonality, temperature of coldest quarter, along with two precipitation variables: precipitation seasonality (Bio15) and precipitation of the coldest quarter (Bio19; Table 2; Figs. 5E–5L). According to the distribution of these predictors and of the opossums, individuals with more depigmented skin and lighter face (distributed above 27°N) were found in localities where temperatures during the coldest quarter (i.e., winter; Fig. 6B) and precipitation seasonality are low (Fig. 6C), while temperature seasonality (Fig. 6A) and precipitation of the coldest quarter are higher (Fig. 6D).

In contrast with other pigmentation traits, the combined models explained relatively less of the variation in torso lightness, for males (14.41%) and females (19.88%), with temperature seasonality, temperature of the coldest quarter and annual precipitation (Bio12) being the most important environmental predictors (Table 2; Figs. 5M, 5N). A pattern of lighter torso coloration was found where temperature seasonality was higher and temperature of the coldest quarter and annual precipitation were lower. These results together with our latitudinal distribution results suggest opossums have high proportion of skin pigmented and darker pelage on the face and torso in humid tropical environments below 26°N, where conditions are warmer, less seasonal and with higher precipitation variability, whereas they are less pigmented towards seasonal, drier and colder habitats above 31°N.

Adaptation to temperate and seasonal environments

Our results indicate that there is a dramatic phenotypic change between tropical and temperate Virginia opossum populations. This change may be driven by natural selection since they adjust to three ecogeographic patterns. The facts that the traits show spatial autocorrelation (Table S4) and that latitude was the most important variable for most of the models (Fig. 5) is not surprising, given the latitudinal gradient we observed in the pattern of phenotypic variation (Table 1) and that environmental variables are also correlated with latitude. Therefore, it is difficult to distinguish the effects of pure spatial autocorrelation when environmental variables that change along similar gradients, which is a limitation of studies perform at large spatial scales. However, the inclusion of latitude and longitude as predictor variables allows us to better tease apart the contributions of geography and ecology separately in our random forest models (Evans et al., 2011; Mascaro et al., 2014) and control for spatial autocorrelation (Table S5). The results of these models support the idea that phenotypic variation in the opossum follow ecogeographic patterns that vary with latitude and identify the most important environmental variables responsible for those patterns.

Koch (1986) detected a similar trend as we observed, that of increasing body size with latitude using the lower first molar area as body size surrogate, but the association to environmental variables was not tested in that study. The association we found of larger hindfoot and body size with environmental predictors of seasonality (Table 2; Figs. 5A, 5B), is consistent with the resource seasonality (also known as fasting endurance) hypothesis as an explanation for Bergmann’s rule. This hypothesis suggests natural selection favors larger individuals in regions of greater seasonality where food availability and energy demands are less predictable (Boyce, 1979; Lindstedt & Boyce, 1985). For the opossum, one of the most important causes of mortality is starvation due to harsh climate conditions in the winter, especially in the northern part of its distribution (Kanda & Fuller, 2004; Kanda, 2005). The main factors predicting whether an opossum will survive the winter are body weight and size (Brocke, 1970; Kanda, 2005). Larger individuals save more energy than smaller ones, owing to a lower metabolism at low temperatures (McNab, 1978). Moreover, the seasonality hypothesis predicts that larger individuals accumulate more fat and metabolize it at lower rates than smaller ones, having greater fasting endurance and survival probability (Lindstedt & Boyce, 1985; Millar & Hickling, 1990). In accordance, under fasting conditions, opossums change from using carbohydrates as energy source to lipid storage (Weber & O’Connor, 2000). Our data suggests that the selective pressure of weather-driven food seasonality and availability may have resulted in a phenotypic adaptation of the Virginia opossum towards larger bodies in seasonal and temperate environments above 27°N. This finding is supported by the fact that at least some Australian marsupials conform with Bergmann’s pattern in response to temperature variability and food availability (Yom-Tov & Nix, 1986; Quin, Smith & Norton, 1996).

The variation of body extremities, especially for tail length, follow Allen’s rule (Figs. 4C, 4D) and is associated with temperature variables (Table 2; Figs. 5C, 5D). Allen’s pattern suggests that natural selection favors individuals with larger body appendages that increase surface area to dissipate heat via conduction in warmer climates, whereas in colder climates shorter appendages would be favored to reduce heat loss (Millien et al., 2006; Tattersall et al., 2012). This pattern is found in two kangaroo species (Yom-Tov & Nix, 1986) while other three Australian marsupials do not show this pattern (Yom-Tov & Nix, 1986; Lindenmayer et al., 1995; Cooper, 1998). Physiological studies have proposed the Virginia opossum is poorly adapted to cold climates due to its high conductivity (i.e., its skin facilitates heat transfer) and low metabolism, and because it relies on behavioral and heat production mechanisms instead of heat conservation as a primary thermoregulatory strategy (Lustick & Lustick, 1972; Hsu, Harder & Lustick, 1988). However, opossums at high latitudes show heat conservation adaptations to seasonal decrease in temperature by developing a higher pelage density during fall and winter months (Gardner, 1973). All physiological studies in this species have used individuals from northern populations above 29°N (i.e., Florida (McNab, 1978), Ohio (Lustick & Lustick, 1972; Hsu, Harder & Lustick, 1988), Michigan (Brocke, 1970) and New York McManus, 1969), which our results revealed are phenotypically more similar among them than compared with southern populations. Consequently, further research is needed to explore the physiological characteristics of tropical populations. As we have shown, northern populations may have phenotypic adaptations to colder climates favoring heat conservation compared to populations in the south. It is possible that during the northwards expansion of the species the new selective pressures imposed by colder environments favored individuals with reduced extremities that conserved heat better.

Finally, the conformance to Gloger’s pattern of skin and coat pigmentation variation and its association with temperature and precipitation variables (Table 2; Figs. 5E–5N) may be driven by pathogens. The pathogen-resistant hypothesis suggests the higher pathogenic incidence in humid and warm tropical environments (Guernier, Hochberg & Guégan, 2004; Lafferty, 2009) is the selective pressure driving the increased pigmentation observed in tropical populations, because highly pigmented skin, hair or feathers confer better resistance to pathogenic infection (Mackintosh, 2001; Burtt & Ichida, 2004). Accumulating evidence suggests that melanocytes, melanosomes and melanin function as integral part of the innate immune system response against invading pathogens (Mackintosh, 2001; Elias, 2007). For example, darker humans are less prone to bacterial and fungal infections than individuals with light skin (Mackintosh, 2001), black feathers on birds are more resistant to bacterial degradation than light ones (Burtt & Ichida, 2004), and darker greenfinches (Carduelis chloris) have higher survival to protozoan infections than paler ones (Männiste & Hõrak, 2014). Similar to the other phenotypic traits we analyzed, the evolution of adaptations in pigmentation (in the opossum’s case is depigmentation) must have occurred during the expansion of the species into North America’s drier and temperate environments at the end of the last Ice Age.

Taken together, our results suggest that the strong geographic variation observed in body size, extremity size and pigmentation traits in the Virginia opossum represent phenotypic adaptations of a species of tropical origins to more seasonal, colder and drier environments. These adaptations may have arisen rapidly, around 15–11 ka (thousand years ago), during the initial phase of the range expansion of the species (Graham et al., 1996; Morgan, 2008; Graham & Lundelius, 2010), which may have facilitated the survival of individuals and the colonization of temperate North America. This is a feasible scenario since rapid adaptation to environmental changes can be accomplished within a few generations (Berteaux et al., 2004; Hairston et al., 2005; Bradshaw & Holzapfel, 2006).

Coloration differences in sexes

As torso and cheek coloration were the coat coloration traits with less variation explained by the environmental variables that we analyzed (Table 2), it is possible there are other factors we did not consider in our analysis, for example, predation, sexual selection or communication, that could also be important for explaining their variation. The fact that, for these two traits alone, females were darker than males, suggest sexual selection may be playing a role, and requires further investigation.

Phenotypic plasticity

There is the possibility that some of the patterns we observed might be achieved via phenotypic plasticity (Ghalambor et al., 2007; Pfennig et al., 2010). In a marsupial species (Sminthopsis crassicaudata), the difference in body size was found to be a response to temperature changes (Riek & Geiser, 2012). Our results for body size variation detected low temperatures as an important predictor, making this alternative hypothesis plausible, however, experimental studies should be done to further investigate the role of plasticity in opossum body size variation. Experimental studies in several mammal species (Ashoub, 1958; Harrison, Morton & Weiner, 1959; Lee, Chu & Chan, 1969; Weaver & Ingram, 1969; Serrat, 2013) have shown that body extremity size of genetically similar individuals (i.e., siblings) varies depending on the temperature at which they are reared, resulting in larger and shorter extremities in warm and cold conditions, respectively. However, the only similar experiment done in a marsupial species did not find differences in extremity size (Riek & Geiser, 2012). Finally, Siamese cats (Iljin & Iljin, 1930) among other mammals (Robinson, 1973) show acrosematic pigmentation, with darker pigmentation on the ears, feet, tail and face, whereas the rest of the body is lighter. This pattern is due to temperature differences in the skin of the appendages and face compared to the core of the body, in which the appendages are the coolest (Stern, 1968). Nevertheless, the pigmentation pattern in the Virginia opossum is the opposite of that predicted by acrosematic patterns. The role of plasticity in the Virginia opossum phenotype variation should be further explored, common garden studies rearing individuals with distinct phenotypes in different environmental conditions and obtaining backcrossed generations may allow to assess the heritability of the traits.