Next Article in Journal
Assessment of Plant Biodiversity and the Floristic Composition in the Black Irtysh River Valley (Kazakhstan)
Previous Article in Journal
Prevalence and Diversity of Plant Parasitic Nematodes in Irish Peatlands
Previous Article in Special Issue
Sex-Based Differences in Multilocus Heterozygosity in Wild Boar from Spain
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Multi-Scale Species Distribution Model for Migrating and Overwintering Western Monarch Butterflies: Climate Is the Best Predictor

by
Ashley R. Fisher
1,2,
William T. Bean
1 and
Francis X. Villablanca
1,*
1
Department of Biological Sciences, California Polytechnic State University, San Luis Obispo, CA 93407, USA
2
The Xerces Society for Invertebrate Conservation, Portland, OR 97232, USA
*
Author to whom correspondence should be addressed.
Diversity 2024, 16(10), 640; https://doi.org/10.3390/d16100640
Submission received: 17 September 2024 / Revised: 8 October 2024 / Accepted: 10 October 2024 / Published: 15 October 2024
(This article belongs to the Collection Feature Papers in Animal Diversity)

Abstract

:
Western Monarch butterflies (Danaus plexippus) migrate from inland breeding ranges to coastal overwintering grounds in California. Given that migratory individuals may make multi-scale habitat selection decisions, we considered a multi-scale species distribution model (SDM) using range-wide climatic and local landscape-level predictors of migratory and overwintering habitat and community-science presence data. The range-wide model output was included as a predictor in the local-scale model, generating multi-scale habitat suitability. The top range-wide predictor was the minimum temperature in December, contributing 83.7% to the model, and was positively associated with presence. At the local scale, the strongest predictors of presence were the range-wide output and percent coverage of low and medium levels of development, contributing > 95%, with 61–63% from the range-wide output, with local-scale suitability coinciding with the California coastal zones. Development’s positive association with overwintering monarch presence was counterintuitive. It is likely that our local-scale model is overfit to these development zones, but it is unclear whether this overfitting resulted from modeler choices, monarchs overwintering close to human development, biased detection near human development, or a combination of these factors. Therefore, alternative approaches to collecting local-scale attribute data are suggested while recognizing the primacy of climate in restricting overwinter sites.

1. Introduction

Appropriate management and conservation require knowledge of species’ habitat needs. However, discerning what characteristics define quality habitats is not as straightforward as describing the attributes at locations where a population or species is found [1,2]. Selection often occurs in a hierarchy of temporal or spatial scales [2,3,4]. As a result, the patchwork distribution of a population in one region of its range is likely driven by a different combination of habitat characteristics than those that determine a species’ entire range. Thus, explicitly discerning the spatial and temporal scales at which environmental characteristics (e.g., climate, food resources, substrate, etc.) are biologically relevant is critical to the study of habitat selection [2,3,4].
Species distribution models (SDMs) are a popular approach for exploring questions around habitat use and selection. However, low-performance or high-uncertainty SDMs can be misleading and, therefore, counterproductive to conservation efforts [5]. SDM quality and performance could be enhanced by integrating multi-scale levels of selection directly into the modeling framework [6]. Bellamy et al. [6] developed a methodology that considers predictors at their biologically relevant scales and incorporates higher-level results into more localized models. For migratory species, where, for example, the geographic range might be subdivided into breeding and overwintering portions, Bellamy et al.’s [6] framework would be consistent with the concept of multi-scale spatial and temporal habitat selection [2,3,4]. This SDM approach allows the explicit consideration of the different spatial and temporal scales at which a species interacts with its environment by discerning which attribute and at which scale they correlate with a high probability of presence while incorporating the effects of higher-order selection into lower-order selection.
We hypothesize that migratory monarch butterflies are an excellent candidate for multi-scale species distribution modeling. In the fall, millions of migrants fly hundreds or thousands of miles south from the northern United States and Southern Canada to reach overwintering groves, where they roost in clusters until remigrating to breed in spring. Migrants east of the Rocky Mountains mostly coalesce in large roosts in high-elevation oyamel fir forests in central Mexico, where millions of migrants are contained within 9 to 14 relatively proximate sites [7]. In contrast, the western population, west of the Rocky Mountains, overwinters dispersed in smaller groves that line the coast from Sonoma County, California, to Baja California Norte, Mexico [8]. Additionally, there are several smaller sites away from the coast, for example, in Arizona and near Death Valley, east of the Sierra Nevada Mountains [8]. A multi-scale SDM could be used to predict the location of unknown sites and to identify appropriate management actions at occupied sites.
Given that overwintering habitat is essential to the Western Monarch overwintering phenomenon, information to preserve this habitat is increasingly necessary and of high interest to conservationists and land managers [8,9,10]. Some relevant grove-level or landscape-level factors might be managed to improve habitat quality; however, the data and analysis would have to be at the appropriate level in order to inform management efforts correctly. But, due to the variable and dispersed nature of these smaller groves, it has been difficult to determine which habitat characteristics are selected or required by migratory monarchs for overwintering in the west and at which scale they matter [11]. Therefore, much of the work describing favorable habitat conditions details within-grove characteristics [11,12,13]. Recently, it has been shown that attributing overwintering habitat characteristics only to the within-grove level is not tenable [11]. For example, latitude has a significant effect on the microclimate at the overwintering cluster sites within individual groves, suggesting that monarchs may be selecting their climatic niche at a landscape rather than within a grove scale. It would follow that broader climatic conditions likely play a key role in limiting this ectotherm’s winter range [11,14,15,16,17]. For instance, northern California and southern Oregon may host amenable tree groves, but prolonged sub-zero temperatures likely make this region unsuitable for migrants by November [18,19,20]. These spatial and temporal considerations associated with migration can be evaluated with a multi-scale modeling approach, which would allow us to assess selection first on habitat suitability in the western region and then integrate the local needs of monarchs. Therefore, we followed [6] and generated a multi-scale nested SDM (using MaxEnt V.3.4) with two levels.
The first and broadest level model examined the full range of the western population through the migration and overwintering season and used climate-based predictors to generate a suitability map of the western United States. This broadest level is referred to as the Full Range (FR) model hereafter. The extent was defined from the known physiological limitations of monarchs [18,21,22,23] as the physiological limit of ectotherms may determine the geographic extents of their distributions [6,23,24,25], and climatic variables have been used as predictors of the geographic distribution of western overwintering groves previously [16]. We also modeled the constricted western population’s range at the peak of overwintering. This scale is based on physical local-level characteristics suggested in the literature. Hereafter, we refer to this level as the Constricted Overwintering Range (COR). The FR output of suitability was integrated into the COR model as a predictor to generate multi-scale predictions for peak overwintering habitat. Aside from the output from FR, COR predictors were landscape-level physical characteristics that could be selected within the home range of overwintering monarchs (Table 1).
The primary goal of the multi-scale model was to explore which habitat characteristics are highly associated with migratory and overwintering monarch presence (inference) and secondarily map suitable habitats [26]. This model may serve as a tool for scientists to understand better the Western Monarch’s winter distribution and for land managers to conserve, manage, and restore monarch habitat. We found that our hypothesis that migratory monarch butterflies are an excellent candidate for a multi-scale SDM is true for the FR scale but not for the COR scale. The presence of human development along the coastal zone compromises the COR scale analysis in multiple possible ways.
Table 1. Predictors considered for FR (population-wide migratory and overwintering) and COR (peak overwintering habitat) model levels. See Supplementary Materials for spatial resolutions and calculation details (Supplementary Materials Table S3).
Table 1. Predictors considered for FR (population-wide migratory and overwintering) and COR (peak overwintering habitat) model levels. See Supplementary Materials for spatial resolutions and calculation details (Supplementary Materials Table S3).
Predictors ConsideredSourceUnits
FRMinimum Temperature of the Coldest month (December) PRISM°C
Maximum Temperature of the Warmest Month (October)
Average Temperature (October to February)
Average Dew Point (humidity) (October to February)
Month with Lowest Minimum VPD (highest humidity) (December)PRISM0–100 hPa
Month with Highest Maximum VPD (lowest humidity) (Octobeer)
Total Seasonal Precipitation (October to February)PRISMmm
Annual Precipitation
Month with Lowest Average Percent Cloud Cover (October)EarthENV% Cover × 0.01
Month with Highest Average Percent Cloud Cover (December)
Average Percent Cloud Cover (October to February)
Month with Lowest Average Wind Speed (October)WORLDCLIMm/s
Month with Highest Average Wind Speed (February)
COR% Forest Land CoverNLCD% Coverage
% Barren Land; Natural Land Cover
% Urban- Low Dev. (Open Space) Land Cover
% Urban- Low/Med Development Land Cover
% Urban- High Development Land Cover
% Cultivated Land Cover
% Shrubs Land Cover
% Wetlands Land Cover
% Unclassified Land Cover
% Open Water Land Cover
% Tree Canopy Cover
Stream DensityNHD% Coverage
Stream Density, Excluding Ephemeral Water Sources
Topographic Position Index (TPI) (aggregation = median)EarthENVMeters
Topographic Ruggedness Index (TRI) (aggregation = mean)
Vector Ruggedness Measure (VRM) (aggregation = mean)EarthENVIndex (0 to 1)
Human footprint[27]Rating (1 to 10)
Logistic Predictive Output from FR ModelFR model0 to 1

2. Materials and Methods

Based on current knowledge, we formed predictions for the FR model from various climatic factors (Table 1). We predicted that minimum temperature [16,18] and high temperatures [23] could be negatively associated with overwintering monarch presence. Low temperatures can result in freezing-related mortality, and high temperature-driven lipid depletion could be physiologically expensive. Cloud cover was predicted to be positively associated with monarch presence, as it may stabilize temperature swings. Moderate precipitation was predicted to be positively associated with roosting and nectar sources, while increased precipitation could be negatively associated with monarch presence via deadly winter storms [20]. Lastly, we included humidity and wind speed as they have been hypothesized to be important within groves [13,28,29].
Land use may reflect habitat loss, a potential driver in the Western Monarch decline [9,10,30]. Therefore, for the COR model, we predicted certain land cover classifications would have positive associations with overwintering monarchs, such as forest and shrublands [30], that provide winter roost sites and winter nectar [14,23,31,32]. Conversely, we predicted that some categories, like moderate/high levels of development, could be negatively associated with monarch presence, implying low overwintering habitat availability. Proximity to freshwater drinking sources has been described as important for overwintering monarchs [14,23,26], as have some aspects of topography, such as selection for shallow canyons or south-to-west facing aspects [12,13]. Thus, we take these descriptive attributes (Table 1) and test their predicted contribution to the model.
We intended the FR model to capture the full climatic niche migrants could occupy while completing migration and overwintering (October to February). As such, the extent of the FR model was the populations’ range, including everything west of the North American Rockies but excluding habitat beyond the geopolitical boundaries of the United States (due to data limitations for Baja California Norte, Mexico). This administrative filter excluded a few presences (<1%), so the model was unlikely to be affected. Because large geographic extents can engender artificially inflated model results [33] and yield less precise relationships between predictors and presences, we used known monarch butterfly physiological limitations to constrict the background extent further [34]. For the FR extent (Figure 1), we conservatively removed regions in the west where monthly mean minimum temperatures fall below a lethal −6 °C in October [18], the warmest month during the migration season, using 4 km resolution PRISM climate data for the month [35] (Table 1).
We constructed the COR model to estimate the habitat attributes of the regions occupied at peak overwintering season. Temporally, we defined peak overwintering season as between migration in fall and remigration in spring, thus spanning from November to February [18,19,20,21,22,23]. Also, we limited the geographic extent of COR (Figure 1) to where the minimum temperature of the coldest month (December) remained above −6 °C [18].
FR included all presence records while monarchs were completing migration and overwintering (October to February), while COR only included records for peak overwintering season (November to February). Presence for both overwintering clusters and lone adult sightings was used. This captured not only places where winter clustering occurred but also pockets of habitat utilized while migrating during fall range contraction.
Presence data from community scientists were obtained from two sources. Overwintering grove presences came primarily from data compiled by Xerces from early statewide records or Thanksgiving Count surveys. Single adult sightings were obtained from Journey North, which collects the sightings of monarchs in various life stages.
Presence data were processed and cleaned using R and ArcMap 10.6 (see Supplementary Materials Table S1) [36,37]. The Maxent.jar algorithm was used to thin points to one per raster cell to reduce spatial autocorrelation [38] such that the final number of occurrence points in FR was 490 and 637 for COR.
The appropriate resolution of our predictors should be at, or smaller than, the scale at which animals selected for the predictors [6]. Based on the estimated daily average rates of migratory movement, a resolution of 4 km was selected for the FR predictors [39]. For COR, we aimed for predictor resolutions to estimate a “home range” after arriving at overwintering grounds (Mayor et al., 2009 [3]). Based on mark-recapture studies of overwintering monarch’s local movements between sites, a predictor resolution of ~5 km or smaller was justified [40,41]. Given that percent coverage values for each predictor vary depending on the modeler’s chosen resolution, two separate scales were considered: 1 km and 4 km. Both distances were feasibly traveled by peak overwintering season monarchs, so at least one would be expected to capture the scale at which selection may occur.
Based on our hypothesis and predictions, the FR climatic variables considered were monthly average wind speed, cloud cover and precipitation, monthly minimum, maximum, and average temperature, and vapor pressure deficit (VPD) from October to February (Table 1) [35,42,43,44]. However, to reduce the number of variables while still capturing variation across the winter and the landscape, we selected months with the highest and lowest average wind speed and average cloud cover, the lowest minimum temperature and VPD (December), and the highest maximum temperature and VPD (October), (see Table 1). We also selected the average (October to February) for wind speed, dew point, temperature, and precipitation, and we selected total annual precipitation.
Monthly climatic data for temperature, VPD, dew point, and precipitation were sourced from PRISM [35]. Monthly mean wind speed (m s−1) data were sourced from WORLDCLIM 2.1 [43] and the mean monthly cloud cover was sourced from EarthENV [44]. See Supplementary Materials (Table S3) for more information on the spatial and temporal resolution of the original source data. Formatting of all spatial data were performed using R and ArcMap 10.8 [36,37]. Raw predictor data (Supplementary Materials Table S2) were resampled in R using bilinear interpolation with the aggregate and disaggregate functions from the raster package to match the 4 km resolution of interest [45].
In addition to the output from FR, COR predictors were landscape-level physical characteristics monarchs are believed to select. Following our hypothesis and predictions, the metrics of topography, human footprint (i.e., development), stream density, and percent coverage of various land cover classes and canopy covers were included as predictors for the model at 1 and 4 km resolutions unless noted otherwise (Table 1). A list of the physical characteristics is offered here with details in Supplementary Materials (Table S4): topographical complexity (Vector Ruggedness Measure (VRM), Topographic Ruggedness Index (TRI), and Topographic Position Index (TPI) (Table 1) [46], human footprint (4 km resolution of [27], stream density (all surface water sources, and only non-ephemeral water sources) [47,48]. In addition, the percent land coverage was generated for ten classes (forest cover, tree canopy cover, shrub cover, low human development, etc.) sourced from the National Land Cover Database (NLCD) (Table 1) [49].
Of the potential biologically meaningful variables listed above, we eliminated several to avoid having multiple collinear predictors in the model [50,51]. For FR, we ran a univariate Maxent model with 10,000 background points and a default regularization multiplier of one to generate jackknife variable importance values. This measures how well each single predictor (i.e., windspeed) predicts monarch presences independent of the other predictor variables [52], thus allowing a ranking of individual variable predictive strength. We generated a correlation matrix between all predictive variables. Using the jackknife values and the correlation matrix, the top variable was selected, and all variables correlated (|r| > 0.7) with that top variable were removed. This process was repeated until all variables were categorized as a “top variable” to be included or discarded because of a high correlation with a top variable. The same process of retention and elimination was employed for COR.
We used Dismo (V.1.3–9) and ENMeval (V. 2.0.3) in R to partition our data and select optimal settings for the final models [53,54,55]. We partitioned background and presence points for both FR and COR via spatial blocking. This is like the random k-fold approach but uses four folds that are partitioned spatially, which is recommended for addressing spatial autocorrelation and biased sampling [56]. All models were generated using the maxent.jar algorithm [38,57,58].
We determined the optimal combination of values for the regularization multiplier (RM) and feature classes, which affect the model complexity and fit [33,52]. Proper optimization prevents model overfitting (poor generalizability and overrepresentation of bias) or underfitting (overestimation of the range of conditions and area with suitable habitat) [50,51]. The RM settings control model penalties for overfitting, with higher RM values leading to higher penalties and a less complex model. Feature classes are the potential shapes the response curves take [53]. Certain feature classes, like threshold (T) or hinge (H), allow for more complexity relative to the linear (L) or quadratic (Q) classes. We used the ENMEval (V. 2.0.3) function to determine which combination of RM value and feature classes performed best. We generated and compared models that combined RM values from 1 to 5 in 0.5 increments and all possible combinations of feature classes available (TLHQP). The combination resulting in the model with the lowest AICc score was selected for use in the final model. If a model with a differing combination fell within 2 AICc points of the top model, it was considered statistically equivalent, and the results were also presented.
We report metrics for the assessment of both model levels, FR and COR. The Area Under the Curve (AUC) measures how well (from 0 to 1) the model can differentiate between conditions at withheld (test) occurrences and of background localities [54]. FR and COR models, run with optimal settings but fit with the full dataset (no partitioning), generate default “cloglog” predictive outputs across geographic extents that identify regions most likely to harbor overwintering monarchs. The cloglog is interpreted as the relative probability of presence (from 0 to 1), where higher suitability values indicate a higher relative probability of monarch presence [52,57]. Additionally, a difference map was generated by subtracting the cloglog output of COR from FR to visualize and compare the change in the suitability of regions from one scale to the other. Univariate response curves are also presented. They are created by fitting a model with a singular predictor to illustrate the relationship of that predictor to the estimated suitability output [38,57].

3. Results

3.1. Final Model Settings and AUC Scores

FR settings with the lowest AICc used an RM of 1, with L and T features, so the model was run with these settings. The mean validation AUC score was 0.93, indicating high accuracy in discerning test observation localities from the background (Table S4). COR settings with the lowest AICc used an RM of 1.5, with L and T features. The next best model had a ΔAICc of 1.4 and used an RM of 1.5, with L T and Q features (Table S4). Both combinations of settings were considered competitive, and results for both are reported below (as LT and LQT models). The LT and LQT models both had validation AUC scores of 0.96.

3.2. Final Models Relative Habitat Suitability: Geographic Heat Maps and Difference Map

The FR cloglog prediction output showed locations with the highest probability of presence situated along the coast of California, within the bays and inlets of northern California, and some inland southern California regions, and less so in Arizona (Figure 2A and Figure 3 for detail). The northernmost hotspot along the coast was near Arcata, California. But consistently high probabilities of presence along the coast are not contiguous until roughly 39° latitude, just north of San Francisco Bay, and then extend to Southern California. The contiguous high probability of presence extended from the coastal regions further inland in southern California than in northern California. The Channel Islands, off the coast of southern California, had a high probability of presence. There were elevated probabilities of presence east of San Diego in the foothills east and west of the Salton Sea and in and around Phoenix, Arizona. Other inland areas showing a slight elevation in the probability of presence included the Central Valley in northern California (near Sacramento), the Colorado River south of Las Vegas, Nevada, and the western foothills of the Sierra Nevada.
The overall area with a high probability of presence was reduced in COR relative to FR (Figure 2). The COR predictive output retained the high probability of presence along the coast from FR but was restricted to nearer the coastline (Figure 2B). The northernmost high probability area near Arcata, California, was retained but reduced. Inland areas with elevated suitability near Phoenix and the foothills to the west of the Salton Sea were also retained but reduced. The moderate probability areas in FR were reduced or eliminated in COR, including the Colorado River section, much of the Central Valley of California, and the western foothills of the Sierra Nevada (Figure 2B). Less developed areas near the coast, such as Big Sur (Figure 4B, between San Luis Obispo and Salinas), the Channel Islands, and the Santa Barbara coast, which were all highly suitable in the FR model, were considered much less suitable in the COR model (Figure 2B and Figure 4). Conversely, inland urban centers such as Las Vegas and Phoenix and all the urban centers in the Central Valley of California (Sacramento, Stockton, Modesto, Fresno, Bakersfield) were ranked as much more suitable in COR than in FR.

3.3. Top Predictors–Percent Contribution and Jackknife Training Gain

The top predictor for FR was the minimum temperature in December, with an 83.7% contribution (Table 2). Cloud cover in December and total precipitation for the overwintering season were the next highest contributing variables, followed by the VPD in December and the average wind speed in February (Table 2). For the Jackknife training gain values, the minimum December temperature was again by far the best predictor of species presence. The next highest contributing variables were VPD and cloud cover in December. Total precipitation and mean wind speed in February had the lowest training gains (Table 2, Figure S3). All predictors not listed in Table 2 present in Table 1 were eliminated due to high correlation.
For COR, the predictor with the top percent contribution of 61.3% and 63.3% across both top models (LT and LQT, respectively) was the nested predictor cloglog output from FR (Table 3). Contrary to our prediction, the second highest contributor at 30.1% and 31.8% for both models was the percent coverage of medium development (rather than the percentage of tree cover). The remaining COR LT predictors all had importance values of 2% or less, while the LQT predictors had importance values of 2.1% or less, but not in the relative order of importance values as for LT. Percent cover low development had a contribution of 1–2% in both COR LT and LQT. All remaining variables had <0.9% contribution (Table 3).
The jackknife training gain values identified the FR output as the strongest predictor of species presence at the COR scale (Table 3). Variables with the next highest training gain (>0.5 and in descending order by training gain in the LT model) were the percent coverage of medium, high, and low-density development, human footprint, percent coverage of cultivated land, tree cover, and open water. Variables with lower training gain values (0.1 to 0.4 and in descending order by training gain in the LT model) include percent coverage of shrubs, TRI, percent coverage of barren land, then wetlands, and stream densities. All other variables had training gain values < 0.1 (Table 3, Figure S3). The LQT model was effectively the same, with only minor variation in the training gain values relative to the LT model.

3.4. Top Predictors–Univariate Response Curves

For the FR model, the minimum temperature in December showed a sharp positive sigmoidal relationship with the probability of presence, with suitability increasing most steeply > 3.0 °C (Figure 5). Both with cloud cover (December) and precipitation (October–February), suitability peaked at very low values and then declined as values increased. Minimum vapor pressure had a positive relationship with the probability of presence at low values and a neutral relationship at high values. Lastly, average wind speed had the most erratic relationship with the relative probability of presence, possibly due to weak predictive power (Figure 5, Table 2).
For the COR model, most of the response curves did not merit interpretation in both analyses due to the exceedingly low percent contribution, except for one attribute-medium development (Table 3). Unexpectedly, increasing values of medium development showed a positive relationship with the probability of monarch presence. The next highest percent contribution values were substantially lower, and given they were for human footprint and percent coverage low development, they were likewise unexpected. Both showed a positive relationship with the probability of monarch presence (Figure 6). Surprisingly, given biological expectations, the percent contribution for percent tree cover was very low (0.2), and the response curves showed a negative association with the probability of presence. (Figure 6 and Figure S4).

4. Discussion

The climatic range-wide model (FR) found that the higher minimum December temperature and lower levels of cloud cover were the most important predictors of monarch presence for the western overwintering populations. At the finer scale, in the peak-overwintering model (COR), urban development was unexpectedly positively correlated with monarch presence despite its negative cost to forested habitat.
A combined view of the FR percent importance, training gain, and univariate response curves led to a comprehensive climatic characterization of the overwintering habitat. Unsurprisingly, minimum temperature was the major predictor, with suitability increasing most quickly above 3.0 °C. This is unsurprising because it has been previously identified as important for monarch butterflies [16,18,59], for overwintering monarch butterfly ecophysiology [23] in particular, and for ectotherms in general [25,60,61,62].
Indeed, monarchs are hypothesized to come to the coast to preserve lipids [23], meaning that they seek the cold while simultaneously avoiding hard freezes. As such, a prediction is that they might seek the coldest non-freezing temperatures in order to preserve the maximum proportion of their lipids. But an important insight from our analysis is that the probability of occurrence was positively correlated with increasing minimum temperature (in the coldest month–December), which reached up to 11.4 °C in the analyzed region, just below a monarch’s flight threshold (Figure 6), meaning there was no evidence that the pressure to preserve lipids was so extreme that risking freezing was a reasonable trade-off in this region.
But being warm runs counter to the hypothesis that staying cool benefits winter energetics as suggested by Chaplin and Wells (1982) [23]. Indeed, we had predicted that high maximum mean daily temperatures would be negatively associated with migrants [16,18,23,59]. But that variable was removed before the final model due to collinearity with the minimum temperature of the coldest month, average temperature, and dew point across all months (October to February). However, in support of our original hypothesis, and relevant here, the univariate response curve for maximum average daily temperatures in the preliminary model runs showed a strong positive trend with the highest probability of presence at ~22 °C and a strong negative relationship with maximum temperatures above ~32 °C (Figure S5). Therefore, it appears that overwintering monarchs have the highest probability of occurrence with higher minimum temperatures and lower maximum temperatures. The lower maximum supports the hypothesis of lipid conservation, and the higher minimum suggests a more reasonable tradeoff—avoiding freezing by staying warm and retaining lipid stores by staying cool.
Both of these optima might be simultaneously achieved with proximity to the coastline. This is because the heat capacity and thermal inertia of the ocean moderate both the minimum and maximum air temperatures at locations proximate to the coast [63,64]. Our analysis provides evidence that this tradeoff could be a mechanistic explanation for why, on a larger regional scale, the Western Monarch overwintering sites are primarily located along the coast of California and for the hypothesis that monarchs migrate along riparian corridors [34,65]. Recognizing that local-scale environmental attributes, including topography, vegetation cover, and cloud cover, affect the thermal inertia and energy flux of coastal waters [64] may provide further insights into the energy dynamics of monarch overwintering on the coast, even potentially defining thermal differences across individual overwintering sites. Collectively, this analysis showed that habitat attributes that exist well beyond the edge of an overwintering grove, such as minimum temperature and cloud cover, show a strong correlation with overwintering presence.
A second characterization of the range-wide (FR) overwintering habitat was that the low levels of cloud cover were positively associated with the overwintering presence. This could be due to higher solar radiation’s contribution to winter energetics, as suggested by Chaplin and Wells (1982) [23] and supported by Saniee, K.; Villablanca (2022) [11]. The specific potential thermoregulatory importance of basking opportunities is discussed in Saniee, K.; Villablanca (2022) [11], where they argue empirically that monarch clustering sites within groves are exposed to short bursts of high-intensity light from the direction of the sunrise. This energy input would reduce the need to use stored energy to bring a monarch’s body up to flight temperature. In addition, low levels of cloud cover could be positively associated with overwintering presence if there is a negative association with storms. Cloud cover also impacts energy flux from seawater to air. Thus, there are multiple mechanistic hypotheses related to low cloud cover.
It was surprising that other climatic predictors, such as precipitation, vapor pressure deficit, and wind speed, had moderate to low additional influence on the population-wide model (Table 2, Figure 5). At this large scale, these variables may be less important than currently thought. Instead, it is still possible they may be important at a finer scale, where locally intense storms are deadly, and local high winds potentially discourage the formation of aggregations [20,28,29].
In contrast to the FR analysis, the nested peak-overwintering habitat model (COR) contained unexpected results that likely do not reflect monarchs’ true habitat preferences. We predicted that a high percent coverage of medium development would negatively influence the probability of species occurrence but found the opposite relationship. Habitat suitability increased even at a 94% cover of medium development (Figure 6). One interpretation leads to the hypothesis that urban development directly enhances the region’s quality as an overwintering habitat. However, this is not likely, given that urban cover increases at a cost to grove habitat.
The alternative and more likely explanation is that the model was overfit, a common problem where training data fit so well that the model does not generalize to novel environmental spaces, thus over-representing bias [66]. In this case, the primary bias overrepresented is likely the association of presence points with urban development. Evidence of this is provided by the difference map (Figure 4). It shows that areas with the largest drop in suitability from FR to COR were Big Sur, rural portions of the Santa Barbara Coast, and the California Channel Islands. These less developed areas are known to support large aggregations of monarchs [13], have had long-term overwinter occupancy [34], and contain 13 of the top 50 sites as ranked by overwintering monarchs’ abundance [9]. The difference map also showed that the suitability increased for COR in the urban areas of Las Vegas, Tucson, Phoenix, the California Central Valley, Sacramento, Stockton, Modesto, Fresno, and Bakersfield. In stark contrast to Big Sur, the Santa Barbara Coast, and the California Channel Islands, and despite there being an abundance of possible observers, only Phoenix [67] and Bakersfield [68] are known to have even a small number of, and certainly not an abundance of, overwintering monarchs.
There are several potential reasons why the COR model was biased towards urban development. One explanation is that the model simply captured the correlation between where monarchs congregate and where humans congregate in the coastal zone. Overwintering monarchs may have a true preference for the coastal zone as most overwintering sites are within five miles of the ocean [13], and as outlined above, there may be a mechanistic reason for this. But given that coastal California is heavily developed (Figure S2), it is possible that the model confounded this coastal association with urban development and extrapolated it to the rest of the extent.
The overfit of the COR level model could also be partly due to observer bias in the presence-only datasets, as community scientist occurrence data are prone to be biased toward heavily populated zones and infrastructure [69]. Bias may explain why the COR heat map seemingly ranks all areas with a high percentage of human development as highly suitable (Figure S2). It is also known that there are non-migratory resident populations in southern, central, and potentially northern California, where breeding occurs year-round primarily on non-native tropical milkweed [70,71,72], meaning that non-migrant lone adults would most likely be spotted in populous urban areas. An observer bias potentially also explains the weak negative association between the probability of presence and high levels of tree cover (Figure 6). As stated above, the rural and remote Big Sur coastline is ranked as a place with high suitability in FR and has several known sites, yet it was predicted unsuitable in COR (Figure 3C and Figure 4), and likewise for the relatively inaccessible Channel Islands [73].
Modeler choices may also have ratcheted the bias overrepresenting a positive association between medium development and overwintering habitat. Specifically, the decision to convert categorical land cover variables to percent cover at chosen scales of 1 km and 4 km; this approach prevented certain issues while potentially exacerbating the poor extrapolation of the COR model. We suspect this because the percent coverage of medium development (e.g., housing) within 1 km and 4 km of overwintering sites, which commonly occur in patches of open space near development, is likely greater than that of the other land cover types such as open grassland or trees. Indeed, where the opposite is true (e.g., Channel Islands and Big Sur), the COR model greatly reduces the probability of occurrence. This was in spite of the fact that we purposefully tried to penalize overfitting through RM while reducing the potential for overfitting by only using biologically meaningful predictors or their proxies and by reducing the number of predictors by removing highly collinear ones.
While distribution models are often positively correlated with metrics of habitat quality, such as density or survival [74,75], suitability and quality are not equivalent. It is possible that urban areas represent ecological traps for monarch sites that promote habitat selection but represent low-quality habitats [76]. For example, they may select sites close to non-native milkweed, which provides oviposition resources throughout the winter but also promotes disease [71,72]. In this case, the high suitability reflected by the COR model might be accurate in assessing a higher probability of presence but misleading in assessing habitat quality [77]. Crone and Schultz (2021) [78] have explored the possibility that a non-migratory monarch population segment that breeds on non-native milkweed could be a cause of the Western Monarch’s decline but find the evidence to date inconclusive.
Regardless, it is impossible to determine whether the cause of this association is due to (the unlikely possibility that) urban areas directly enhance a region’s quality as overwintering habitat, the presence of data bias, a true preference for the coastal zone leading to a confounding relationship with urban development, or some combination of these non-mutually exclusive hypotheses. Thus, we conclude that overall, the COR output is likely not reflective of monarchs’ true habitat preference.
To ameliorate overfitting to urban cover in the peak-overwintering habitat model (COR), we ran a post-hoc analysis where human footprint was incorporated as a bias file rather than a predictor (Supplementary Materials Figures S6–S8). The results showed that the percent contribution for medium development cover dropped from 31 to ~18%. However, the predicted suitability across the West appeared to be heavily still weighted toward developed locations, and relative contributions from other predictors remained similarly low. For example, medium, high, and low development were still the most important predictors (following the FR contribution). Unfortunately, our ad hoc approach did not mitigate overfitting.
Unfortunately, we reject our hypothesis that migratory monarch butterflies are an excellent candidate for multi-scale species distribution modeling at this time. Collectively, we would urge caution in interpreting the COR level analysis given the potential for overfitting and lack of generalization, correlations among attributes in the monarch/human coastal congregation zone, scales of available data, and/or of formatting to categorical variables to percent coverage at 1 km and 4 km or otherwise, and presence-only community-science sampling bias.
As an alternative approach, we suggest that simply mapping the spatial intersection of groves of trees with the FR model would provide useful multi-scale modeling. An intersection map would provide the predicted locations of possible overwintering sites (known and not yet known). It would include the FR output in areas with values of 0.7 or higher, which collectively contributed the most to the COR model, and a relevant habitat attribute (tree grove) at whatever scale each grove presents itself. Furthermore, there could be significant benefits from systematic and/or randomized sampling from this intersection map for overwintering habitats outside of heavily populated zones. This approach has the potential to improve our understanding of species distribution in a non-biased manner and as a complement to the current community-science-based datasets.
According to the FR output, the extent to which inland highly suitable habitats extend from the coastline is geographically variable. In the north, it appears that highly suitable habitats follow the coast narrowly, including bays and their coastal valleys (Figure 2, Figure 3 and Figure S1). As an example, San Jose, which lies within San Francisco Bay, appears to have suitability values comparable to those found close to the Pacific coast. Further south in Santa Cruz, Monterey, and San Luis Obispo counties, highly suitable regions appear close to the coast where they are possibly constrained by coastal mountains (i.e., the Santa Lucia range and the Santa Ynez mountains). Where there are no coastal mountains, suitability comparable to that found close to the coast extends into coastal valleys (examples include the Salinas, Santa Maria, and Lompoc valleys) (Figure 2 and Figure 3). In southern coastal California, suitable habitat reaches further inland, especially in low-lying areas and valleys like the Los Angeles basin, where it extends farther inland than in any other region adjacent to coastal habitat.
Other areas of notably high probability of presence, and therefore high interest, include the Channel Islands and mountains west of the Salton Sea (Figure 2 and Figure 3). Both areas were lower in population density but were still selected. On the northern/opposite end of the population’s range, Arcata, California, also ranked as a location that had high suitability (Figure 2 and Figure 3). This could be plausible given its location within a bay, but though our model did not directly capture the intensity or frequency of winter storms, we anticipate they could be a limiting factor at this latitude. It is possible that as a mitigating factor, some of these populations (if present) might be located right along the coast for thermal benefits, in dense forests, for wind abatement benefits, or both.
We suggest that all the inland, island, and high latitude regions that previously have not been given high priority for monarch surveys be given more attention in a year with high population numbers. The discovery of even transitional sites in these “unorthodox” locations could shed insights into the breadth of the ecological niche overwintering monarchs select. Considering that inland sites in Death Valley, Bakersfield, and Arizona exist, our model may provide a framework for predicting where systematic sampling for other overwintering sites or overwintering site attributes might be fruitful. Alternatively, this framework could come not only from the FR model but from a combination of mapping the FR model and its spatial intersection with the groves of trees.
As with all models, ours should be interpreted with its assumptions and limitations in mind. Not all predictors potentially relevant to overwintering monarch habitat could be included. For example, storm frequency, intensity, and pesticide application were all considered because of their potential relationship with monarch mortality [29,79]. However, no sources were found that had the appropriate extent or resolution for this analysis.
There is no way for us to truly discern which lone adult occurrences are true migrants and which are potentially non-migratory year-round breeders that do not exhibit overwintering behavior or use overwintering habitat [41,71,72,73]. This issue may be of most concern in and near urban areas, where tropical milkweed is commonly planted and available year-round [71,72].

5. Conclusions

Western Monarchs are under threat of a long-term decline [9,80]. Protecting overwintering habitat is a key step in conservation planning for the long-term success and recovery of the western migratory population [10,81]. The multi-scale species distribution model presented here took a large-scale approach to identifying the climatic and local landscape-level predictors of quality overwintering habitat. We found that minimum temperature is the strongest climatic predictor of presence across the species’ western migration and winter range. We continue to find that the habitat with the highest probability of occurrence lies primarily along the California Coast. We argue that this makes sense given monarch butterfly energetic constraints. In addition, we continue to develop evidence that the salient characteristics of overwintering groves may be defined by the attributes found outside of the grove. On a more localized scale, the strongest predictors were the percent coverage for low and medium levels of development, which had a strong positive association with the probability of overwintering monarch presence. This result is likely an artifact of overfitting, and we infer that the local-scale model is not a good representation of quality overwintering habitats. The multiple possible factors leading to overfitting could be mitigated by more systematic sampling of less developed regions for overwintering monarch presence and habitat attributes. This might effectively be conducted by overlaying our climatic and landscape (FR) model and the location of groves to identify possible sampling sites. Subsequent randomized or systematic sampling for overwintering monarch presence would help us gain an unbiased understanding of overwintering habitat preference.

Supplementary Materials

The following supporting information can be downloaded at https://www.mdpi.com/article/10.3390/d16100640/s1: Figure S1. FR (A) and COR (B) predicted suitability estimates in Central and Southern California, and parts of Arizona, with major cities and roads overlayed for reference. Figure S2. Percent coverage Medium Development (A) versus COR predictive output (B). Figure S3. Jackknife Variable Importance under optimal settings for FR and COR. Figure S4. Marginal and Univariate Response Curves from the peak overwintering range (COR) model, run with L Q and T feature classes. Figure S5. Univariate Response Curve for maximum temperature in the hottest month, October (tmax_10), from the preliminary model run with all candidate variables. Figure S6. Percent contribution of the COR predictors in the run with human footprint used as a bias file instead of a predictor. Figure S7. Jackknife Variable importance results of COR predictors run with human footprint used as a bias file instead of a predictor. Figure S8. Cloglog Predictive output of COR predictors run with human footprint used as a bias file instead of a predictor. Figure S9. Comparison of the output from the two top-performing models from the COR Model. Table S1. Presence data provided and subsequently removed through cleaning processes described in the methods. Table S2. Source and resolution of raw predictor data. Table S3. Reclassification for National Landcover Database (NLCD) categories into our 10 Percent coverage Categories. Table S4. Optimal Model Settings and performance of top five model setting combinations for FR and COR candidate models. Model settings altered between runs were the Regularization Multiplier (from 1 to 5 in 0.5 increments) and Feature Classes, including Linear (L), Quadratic (Q), Threshold (T), and Hinge (H). The average validation AUC for each Model is reported. Optimal model settings, denoted with a (*), were those with the lowest AICc score or within 2 points of the lowest AICc score.

Author Contributions

Conceptualization, A.R.F., W.T.B. and F.X.V.; methodology, A.R.F., W.T.B. and F.X.V.; software, A.R.F., W.T.B. and F.X.V.; validation, A.R.F., W.T.B. and F.X.V.; formal analysis, A.R.F.; investigation, F.X.V.; resources, F.X.V.; data curation, A.R.F.; writing—original draft preparation, A.R.F.; writing—review and editing, W.T.B. and F.X.V.; visualization, A.R.F.; supervision, F.X.V.; project administration, F.X.V.; funding acquisition, F.X.V. All authors have read and agreed to the published version of the manuscript.

Funding

Funding was provided by the U.S.F.W.S Coastal Program (Grant #F20AC11996-00) and by The College of Science and Mathematics at Cal Poly State University, San Luis Obispo.

Institutional Review Board Statement

Not applicable.

Data Availability Statement

Restrictions apply to the availability of the distribution data. Data were obtained from Xerces and Journey North and are available from the authors with permission from Xerces and Journey North. The predicted model outputs are hosted at Figshare. Data for the raw predictor layers can be found in the cited sources.

Acknowledgments

The Xerces Society for Invertebrate Conservation and Journey North provided the presence data, which we acknowledge were collected by a small army of volunteer community scientists to whom we are grateful. Emma Pelton and Isis Howard of Xerces, Jessica Griffiths of Colorado State University’s Center for Environmental Management of Military Lands, Dan Meade of Althouse and Meade, and Samantha Marcum of the U.S. Fish and Wildlife Service contributed to the conceptualization of this research.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Morrison, M.L.; Marcot, B.G.; Mannan, R.W. Wildlife-Habitat Relationships: Concepts and Applications, 3rd ed.; Island Press: Washington, DC, USA, 2006. [Google Scholar]
  2. Johnson, D.H. The Comparison of Usage and Availability Measurements for Evaluating Resource Preference. Ecology 1980, 61, 65–71. [Google Scholar] [CrossRef]
  3. Mayor, S.J.; Schneider, D.C.; Schaefer, J.A.; Mahoney, S.P. Habitat selection at multiple scales. Ecoscience 2009, 16, 238–247. [Google Scholar] [CrossRef]
  4. Wiens, J.A. Spatial Scaling in Ecology. Funct. Ecol. 1989, 3, 385. [Google Scholar] [CrossRef]
  5. Loiselle, B.A.; Howell, C.A.; Graham, C.H.; Goerck, J.M.; Brooks, T.; Smith, K.G.; Williams, P.H. Avoiding Pitfalls of Using Species Distribution Models in Conservation Planning. Conserv. Biol. 2003, 17, 1591–1600. [Google Scholar] [CrossRef]
  6. Bellamy, C.; Boughey, K.; Hawkins, C.; Reveley, S.; Spake, R.; Williams, C.; Altringham, J. A sequential multi-level framework to improve habitat suitability modelling. Landsc. Ecol. 2020, 35, 1001–1020. [Google Scholar] [CrossRef]
  7. Rendón-Salinas, E.; Martínez-Meza, F.; Mendoza-Pérez, M.; Cruz-Piña, M.; Mondragon-Contreras, G. Superficie Forestal Ocupada Por Las Colonias De Mariposas Monarca En México Durante La Hibernación De 2018–2019; World Wildlife Fund México: Mexico City, Mexico, 2019. [Google Scholar]
  8. Xerces Society Western Monarch Count. Western Monarch Thanksgiving Count and New Year’s Count Data, 1997–2021. 2022. Available online: www.westernmonarchcount.org (accessed on 10 October 2024).
  9. Pelton, E.; Jepsen, S.J.; Schultz, C.B.; Fallon, C.; Black, S.H. State of the Monarch Butterfly Overwintering Sites in California|Xerces Society; The Xerces Society for Invertebrate Conservation: Portland, OR, USA, 2016; Available online: https://www.xerces.org/publications/scientific-reports/state-of-monarch-butterfly-overwintering-sites-in-california (accessed on 21 November 2023).
  10. Pelton, E.M.; Schultz, C.B.; Jepsen, S.J.; Black, S.H.; Crone, E.E. Western Monarch Population Plummets: Status, Probable Causes, and Recommended Conservation Actions. Front. Ecol. Evol. 2019, 7, 258. [Google Scholar] [CrossRef]
  11. Saniee, K.; Villablanca, F. Hierarchy and Scale Influence the Western Monarch Butterfly Overwintering Microclimate. Front. Conserv. Sci. 2022, 3, 844299. [Google Scholar] [CrossRef]
  12. Lane, J. Overwintering Monarch Butterflies in California: Past and Present. In The Monarch Butterfly: Biology and Conservation; Malcolm, S.B., Zalucki, M.P., Eds.; Cornell University Press: Ithaca, NY, USA, 1993. [Google Scholar]
  13. Leong, K.L.H.; Sakai, W.; Bremer, W.; Feuerstein, D.; Yoshimura, G. Analysis of the Pattern of Distribution and Abundance of Monarch Overwintering Sites Along the California Coastline. In The Monarch Butterfly: Biology and Conservation; Oberhauser, K., Ed.; Cornell University Press: Ithaca, NY, USA, 2004. [Google Scholar]
  14. Tuskes, P.M.; Brower, L.P. Overwintering ecology of the monarch butterfly, Danaus plexippus L., in California. Ecol. Entomol. 1978, 3, 141–153. [Google Scholar] [CrossRef]
  15. Batalden, R.V.; Oberhauser, K.; Peterson, A.T. Ecological niches in sequential generations of Eastern North American monarch butterflies (Lepidoptera: Danaidae): The ecology of migration and likely climate change implications. Environ. Entomol. 2007, 36, 1365–1373. [Google Scholar] [CrossRef]
  16. Fisher, A.; Saniee, K.; Van der Heide, C.; Griffiths, J.; Meade, D.; Villablanca, F. Climatic Niche Model for Overwintering Monarch Butterflies in a Topographically Complex Region of California. Insects 2018, 9, 167. [Google Scholar] [CrossRef]
  17. Oberhauser, K.; Peterson, A.T. Modeling current and future potential wintering distributions of eastern North American monarch butterflies. Proc. Natl. Acad. Sci. USA 2003, 100, 14063–14068. [Google Scholar] [CrossRef] [PubMed]
  18. Anderson, J.B.; Brower, L.P. Freeze-protection of overwintering monarch butterflies in Mexico: Critical role of the forest as a blanket and an umbrella. Ecol. Entomol. 1996, 21, 107–116. [Google Scholar] [CrossRef]
  19. Brower, L.P.; Williams, E.H.; Fink, L.S.; Slayback, D.A.; Ramírez, M.I.; García, M.V.L.; Zubieta, R.R.; Weiss, S.B.; Calvert, W.H.; Zuchowski, W. Overwintering Clusters of the Monarch Butterfly Coincide with the Least Hazardous Vertical Temperatures in the Oyamel Forest. J. Lepidopterists' Soc. 2011, 65, 27–46. [Google Scholar] [CrossRef]
  20. Calvert, W.H.; Zuchowski, W.; Brower, L.P. The Effect of Rain, Snow and Freezing Temperatures on Overwintering Monarch Butterflies in Mexico. Biotropica 1983, 15, 42–47. [Google Scholar] [CrossRef]
  21. Urquhart, F.A.; Urquhart, N. The overwintering site of the eastern population of the monarch butterfly (Danaus p. plexipus; Danaidae) in southern Mexico. J. Lepid. Soc. 1976, 30, 153–158. [Google Scholar]
  22. Urquhart, F.A.; Urquhart, N. Overwintering Areas And Migratory Routes Of The Monarch Butterfly (Danaus p. plexippus, Lepidoptera: Danaidae) In North America, With Special Reference To The Western Population. Can. Entomol. 1977, 109, 1583–1589. [Google Scholar] [CrossRef]
  23. Chaplin, S.B.; Wells, P.H. Energy reserves and metabolic expenditures of monarch butterflies overwintering in southern California. Ecol. Entomol. 1982, 7, 249–256. [Google Scholar] [CrossRef]
  24. Baker, R.H.A.; Sansford, C.E.; Jarvis, C.H.; Cannon, R.J.C.; MacLeod, A.; Walters, K.F.A. The role of climatic mapping in predicting the potential geographical distribution of non-indigenous pests under current and future climates. Agric. Ecosyst. Environ. 2000, 82, 57–71. [Google Scholar] [CrossRef]
  25. Forister, M.L.; Halsch, C.A.; Nice, C.C.; Fordyce, J.A.; Dilts, T.E.; Oliver, J.C.; Prudic, K.L.; Shapiro, A.M.; Wilson, J.K.; Glassberg, J. Fewer butterflies seen by community scientists across the warming and drying landscapes of the American West. Science 2021, 371, 1042–1045. [Google Scholar] [CrossRef]
  26. Zurell, D.; Franklin, J.; König, C.; Bouchet, P.J.; Dormann, C.F.; Elith, J.; Fandos, G.; Feng, X.; Guillera-Arroita, G.; Guisan, A.; et al. A standard protocol for reporting species distribution models. Ecography 2020, 43, 1261–1277. [Google Scholar] [CrossRef]
  27. Leu, M.; Hanser, S.E.; Knick, S.T. The Human Footprint in the West: A Large-Scale Analysis of Anthropogenic Impacts. Ecol. Appl. 2008, 18, 1119–1139. [Google Scholar] [CrossRef] [PubMed]
  28. Leong, K.L.H. Microenvironmental Factors Associated with the Winter Habitat of the Monarch Butterfly (Lepidoptera: Danaidae) in Central California. Ann. Entomol. Soc. Am. 1990, 83, 906–910. [Google Scholar] [CrossRef]
  29. Leong, K.L.H.; Frey, D.; Brenner, G.; Baker, S.; Fox, D. Use of Multivariate Analyses to Characterize the Monarch Butterfly (Lepidoptera: Danaidae) Winter Habitat. Ann. Entomol. Soc. Am. 1991, 84, 263–267. [Google Scholar] [CrossRef]
  30. Crone, E.E.; Pelton, E.M.; Brown, L.M.; Thomas, C.C.; Schultz, C.B. Why are monarch butterflies declining in the West? Understanding the importance of multiple correlated drivers. Ecol. Appl. 2019, 29, e01975. [Google Scholar] [CrossRef] [PubMed]
  31. Griffiths, J.; Villablanca, F. Managing monarch butterfly overwintering groves: Making. Calif. Fish Game 2015, 101, 11. [Google Scholar]
  32. Masters, A.R.; Malcolm, S.B.; Brower, L.P. Monarch Butterfly (Danaus plexippus) Thermoregulatory Behavior and Adaptations for Overwintering in Mexico. Ecology 1988, 69, 458–467. [Google Scholar] [CrossRef]
  33. Merow, C.; Smith, M.J.; Silander, J.A. A practical guide to MaxEnt for modeling species’ distributions: What it does, and why inputs and settings matter. Ecography 2013, 36, 1058–1069. [Google Scholar] [CrossRef]
  34. Dingle, H.; Zalucki, M.P.; Rochester, W.A.; Armijo-Prewitt, T. Distribution of the monarch butterfly, Danaus plexippus (L.) (Lepidoptera: Nymphalidae), in western North America: Monarch Butterflies in Western North America. Biol. J. Linn. Soc. 2005, 85, 491–500. [Google Scholar] [CrossRef]
  35. PRISM Climate Group. PRISM Climate Group, Oregon State University. 2014. Available online: https://prism.oregonstate.edu/ (accessed on 15 June 2021).
  36. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2013; Available online: http://www.R-project.org/ (accessed on 12 October 2022).
  37. ESRI. ArcMap, version 10.8; Environmental Systems Research Institute: Redlands, CA, USA, 2020. [Google Scholar]
  38. Phillips, S.J.; Anderson, R.P.; Schapire, R.E. Maximum entropy modeling of species geographic distributions. Ecol. Model. 2006, 190, 231–259. [Google Scholar] [CrossRef]
  39. Billings, J. Opening a window on southwestern monarchs: Fall migrant monarch butterflies, Danaus plexippus (L.), tagged synchronously in southeastern Arizona migrate to overwintering regions in either southern California or Central Mexico. J. Lepid. Soc. 2019, 73, 257–267. [Google Scholar] [CrossRef]
  40. Griffiths, J.L. Monarch butterfly (Danaus plexippus) Tree Preference and Intersite Movement at California Overwintering Sites. Master’s Thesis, California Polytechnic State University, San Luis Obispo, CA, USA, 2014. Available online: https://digitalcommons.calpoly.edu/cgi/viewcontent.cgi?article=2308&context=theses (accessed on 14 December 2014).
  41. James, D.G.; Kappen, L. Further Insights on the Migration Biology of Monarch Butterflies, Danaus plexippus (Lepidoptera: Nymphalidae) from the Pacific Northwest. Insects 2021, 12, 161. [Google Scholar] [CrossRef] [PubMed]
  42. Daly, C.; Halbleib, M.; Smith, J.I.; Gibson, W.P.; Doggett, M.K.; Taylor, G.H.; Curtis, J.; Pasteris, P.P. Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States. Int. J. Climatol. 2008, 28, 2031–2064. [Google Scholar] [CrossRef]
  43. Fick, S.; Hijmans, R. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 2017, 37, 4302–4315. [Google Scholar] [CrossRef]
  44. Wilson, A.M.; Jetz, W. Remotely Sensed High-Resolution Global Cloud Dynamics for Predicting Ecosystem and Biodiversity Distributions. PLoS Biol. 2016, 14, e1002415. [Google Scholar] [CrossRef] [PubMed]
  45. Hijmans, R.J. Raster: Geographic Data Analysis and Modeling. (R Package). 2022. Available online: https://CRAN.R-project.org/package=raster (accessed on 5 January 2023).
  46. Amatulli, G.; Domisch, S.; Tuanmu, M.-N.; Parmentier, B.; Ranipeta, A.; Malczyk, J.; Jetz, W. A suite of global, cross-scale topographic variables for environmental and biodiversity modeling. Sci. Data 2018, 5, 180040. [Google Scholar] [CrossRef]
  47. Moore, R.B.; McKay, L.D.; Rea, A.H.; Bondelid, T.R.; Price, C.V.; Dewald, T.G.; Johnston, C.M. User’s Guide for the National Hydrography Dataset Plus (NHDPlus) High Resolution U.S. Geological Survey Open-File Report. In Open-File Report (No. 2019-1096); U.S. Geological Survey: Reston, VA, USA, 2019. [Google Scholar] [CrossRef]
  48. U.S. Geological Survey. National Hydrography Dataset for California State, 2.2.1; U.S. Geological Survey: Reston, VA, USA, 2020.
  49. Dewitz, J. National Land Cover Database (NLCD) 2016 Products: U.S. Geological Survey Data Release; U.S. Geological Survey: Reston, VA, USA, 2019. [CrossRef]
  50. Warren, D.L.; Wright, A.N.; Seifert, S.N.; Shaffer, H.B. Incorporating model complexity and spatial sampling bias into ecological niche models of climate change risks faced by 90 California vertebrate species of concern. Divers. Distrib. 2014, 20, 334–343. [Google Scholar] [CrossRef]
  51. Warren, D.L.; Seifert, S.N. Ecological niche modeling in Maxent: The importance of model complexity and the performance of model selection criteria. Ecol. Appl. 2011, 21, 335–342. [Google Scholar] [CrossRef]
  52. Elith, J.; Phillips, S.J.; Hastie, T.; Dudík, M.; Chee, Y.E.; Yates, C.J. A statistical explanation of MaxEnt for ecologists. Divers. Distrib. 2011, 17, 43–57. [Google Scholar] [CrossRef]
  53. Hijmans, R.J.; Phillips, S.; Leathwick, J.; Elith, J. dismo: Species Distribution Modeling, R package version 1.3-9; 2023; Available online: https://github.com/rspatial/dismo (accessed on 10 October 2024).
  54. Muscarella, R.; Galante, P.J.; Soley-Guardia, M.; Boria, R.A.; Kass, J.M.; Uriarte, M.; Anderson, R.P. ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods Ecol. Evol. 2014, 5, 1198–1205. [Google Scholar] [CrossRef]
  55. Warren, D.L.; Glor, R.E.; Turelli, M. ENMTools: A toolbox for comparative studies of environmental niche models. Ecography 2010, 33, 607–611. [Google Scholar] [CrossRef]
  56. Roberts, D.R.; Bahn, V.; Ciuti, S.; Boyce, M.S.; Elith, J.; Guillera-Arroita, G.; Hauenstein, S.; Lahoz-Monfort, J.J.; Schröder, B.; Thuiller, W.; et al. Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. Ecography 2017, 40, 913–929. [Google Scholar] [CrossRef]
  57. Phillips, S.J.; Dudík, M. Modeling of species distributions with Maxent: New extensions and a comprehensive evaluation. Ecography 2008, 31, 161–175. [Google Scholar] [CrossRef]
  58. Phillips, S.J.; Anderson, R.P.; Dudík, M.; Schapire, R.E.; Blair, M.E. Opening the black box: An open-source release of Maxent. Ecography 2017, 40, 887–893. [Google Scholar] [CrossRef]
  59. Brower, L. Understanding And Misunderstanding The Migratnon of the Monarch Butterfly (Nymphalidae) in North America: 1857–1995. J. Lepid. Soc. 1995, 49, 304–385. [Google Scholar]
  60. Hodek, I.; Hodková, M. Multiple role of temperature during insect diapause: A review. Entomol. Exp. Appl. 1988, 49, 153–165. [Google Scholar] [CrossRef]
  61. Huey, R.B.; Kingsolver, J.G. Climate Warming, Resource Availability, and the Metabolic Meltdown of Ectotherms. Am. Nat. 2019, 194, E140–E150. [Google Scholar] [CrossRef] [PubMed]
  62. Koštál, V. Eco-physiological phases of insect diapause. J. Insect Physiol. 2006, 52, 113–127. [Google Scholar] [CrossRef]
  63. Abraham, J.; Cheng, L.; Mann, M.E.; Trenberth, K.; von Schuckmann, K. The ocean response to climate change guides both adaptation and mitigation efforts. Atmos. Ocean. Sci. Lett. 2022, 15, 100221. Available online: https://www.sciencedirect.com/science/article/pii/S1674283422000964 (accessed on 22 July 2024). [CrossRef]
  64. Kahle, A.B.; Schieldge, J.P.; Alley, R.E. Sensitivity of thermal inertia calculations to variations in environmental factors. Remote Sens. Environ. 1984, 16, 211–232. Available online: https://www.sciencedirect.com/science/article/abs/pii/0034425784900658 (accessed on 22 July 2024). [CrossRef]
  65. Pyle, R.M. Chasing Monarchs: Migrating with the Butterflies of Passage; Yale University Press: New Haven, CT, USA, 2014. [Google Scholar]
  66. Ying, X. An Overview of Overfitting and Its Solutions. J. Phys. Conf. Ser. 2019, 1168, 022022. [Google Scholar] [CrossRef]
  67. Morris, G.M.; Kline, C.; Morris, S.M. Status of Danaus plexippus population in Arizona. J. Lepid. Soc. 2015, 69, 91–107. [Google Scholar] [CrossRef]
  68. Frey, D.; Schaffner, A. Spatial and Temporal Pattern of Monarch Overwintering Abundance in Western North America. In The Monarch Butterfly: Biology and Conservation; Cornell University Press: Ithaca, NY, USA, 2004; pp. 167–176. [Google Scholar]
  69. Geldmann, J.; Heilmann-Clausen, J.; Holm, T.E.; Levinsky, I.; Markussen, B.; Olsen, K.; Rahbek, C.; Tøttrup, A.P. What determines spatial bias in citizen science? Exploring four recording schemes with different proficiency requirements. Divers. Distrib. 2016, 22, 1139–1149. [Google Scholar] [CrossRef]
  70. James, D.G.; Schaefer, M.C.; Krimmer Easton, K.; Carl, A. First Population Study on Winter Breeding Monarch Butterflies, Danaus plexippus (Lepidoptera: Nymphalidae) in the Urban South Bay of San Francisco, California. Insects 2021, 12, 946. [Google Scholar] [CrossRef]
  71. Satterfield, D.A.; Maerz, J.C.; Altizer, S. Loss of migratory behaviour increases infection risk for a butterfly host. Proc. R. Soc. B Biol. Sci. 2015, 282, 20141734. [Google Scholar] [CrossRef] [PubMed]
  72. Satterfield, D.A.; Villablanca, F.; Maerz, J.; Altizer, S. Migratory monarchs wintering in California experience low infection risk compared to monarchs breeding year-round on non-native milkweed. Integr. Comp. Biol. 2016, 56, 343–352. [Google Scholar] [CrossRef]
  73. Wenner, A.M.; Harris, A.M. Do California Monarchs Undergo Long-Distance Directed Migration; Natural History Museum of Los Angeles County: Los Angeles, CA, USA, 1993; Volume 38, pp. 209–218. [Google Scholar]
  74. Bean, W.T.; Prugh, L.R.; Stafford, R.; Butterfield, H.S.; Westphal, M.; Brashares, J.S. Species distribution models of an endangered rodent offer conflicting measures of habitat quality at multiple scales. J. Appl. Ecol. 2014, 51, 1116–1125. [Google Scholar] [CrossRef]
  75. Weber, M.M.; Stevens, R.D.; Diniz-Filho, J.A.F.; Grelle, C.E.V. Is there a correlation between abundance and environmental suitability derived from ecological niche modelling? A meta-analysis. Ecography 2017, 40, 817–828. [Google Scholar] [CrossRef]
  76. Battin, J. When good animals love bad habitats: Ecological traps and the conservation of animal populations. Conserv. Biol. 2004, 18, 1482–1491. [Google Scholar] [CrossRef]
  77. Semerdjian, A.E.; Butterfield, H.S.; Stafford, R.; Westphal, M.F.; Bean, W.T. Combining occurrence and habitat suitability data improve conservation guidance for the giant kangaroo rat. J. Wildl. Manag. 2021, 85, 855–867. [Google Scholar] [CrossRef]
  78. Crone, E.E.; Schultz, C.B. Resilience or catastrophe? A possible state change for monarch butterflies in western North America. Ecol. Lett. 2021, 24, 1533–1538. [Google Scholar] [CrossRef]
  79. Halsch, C.A.; Code, A.; Hoyle, S.M.; Fordyce, J.A.; Baert, N.; Forister, M.L. Pesticide Contamination of Milkweeds Across the Agricultural, Urban, and Open Spaces of Low-Elevation Northern California. Front. Ecol. Evol. 2020, 8, 162. [Google Scholar] [CrossRef]
  80. Schultz, C.B.; Brown, L.M.; Pelton, E.; Crone, E.E. Citizen science monitoring demonstrates dramatic declines of monarch butterflies in western North America. Biol. Conserv. 2017, 214, 343–346. [Google Scholar] [CrossRef]
  81. WAFWA. Western Monarch Butterfly Conservation Plan 2019–2069; Western Association of Fish and Wildlife Agencies (WAFWA): Boise, ID, USA, 2019. [Google Scholar]
Figure 1. FR and COR model extents and Western Monarch presence points. FR, or the population-wide level (hatched green and blue), extended to the land west of the Rocky Mountains, where the minimum temperature of the warmest month during migration (October) is above −6 °C. The COR model (blue) extended to the land west of the Rocky Mountains, where the minimum temperature of the coldest month (December) was above −6 °C. Both extents were constricted to the United States borders due to data limitations. Presence points in FR (both yellow and black) included all observations of adults and overwintering clusters from October to February. Presence points in the COR model (black only) include observations of adults and overwintering clusters during peak overwintering season only (November to February).
Figure 1. FR and COR model extents and Western Monarch presence points. FR, or the population-wide level (hatched green and blue), extended to the land west of the Rocky Mountains, where the minimum temperature of the warmest month during migration (October) is above −6 °C. The COR model (blue) extended to the land west of the Rocky Mountains, where the minimum temperature of the coldest month (December) was above −6 °C. Both extents were constricted to the United States borders due to data limitations. Presence points in FR (both yellow and black) included all observations of adults and overwintering clusters from October to February. Presence points in the COR model (black only) include observations of adults and overwintering clusters during peak overwintering season only (November to February).
Diversity 16 00640 g001
Figure 2. (A) FR (population-wide) model, and (B) COR (peak overwintering habitat) model predictive outputs. Values from 0 to 1 indicate the relative probability of presence for monarch butterfly overwintering sites. A higher-resolution view of Figure 2A is presented in Figure 3. The COR model shown is the output of the LT model. For a higher resolution view of Figure 2B, see Supplemental Material Figure S1. See Supplementary Materials (Table S4) for a view of the LQT output (Figure S9), which performed equally well and is nearly identical.
Figure 2. (A) FR (population-wide) model, and (B) COR (peak overwintering habitat) model predictive outputs. Values from 0 to 1 indicate the relative probability of presence for monarch butterfly overwintering sites. A higher-resolution view of Figure 2A is presented in Figure 3. The COR model shown is the output of the LT model. For a higher resolution view of Figure 2B, see Supplemental Material Figure S1. See Supplementary Materials (Table S4) for a view of the LQT output (Figure S9), which performed equally well and is nearly identical.
Diversity 16 00640 g002
Figure 3. FR predictive output. With (AE) images at a larger scale (~1:4,000,000) for regions with high suitability estimates from the FR model (overwintering and migratory). With subpanels (AE) corresponding with inset boxes (AE).
Figure 3. FR predictive output. With (AE) images at a larger scale (~1:4,000,000) for regions with high suitability estimates from the FR model (overwintering and migratory). With subpanels (AE) corresponding with inset boxes (AE).
Diversity 16 00640 g003
Figure 4. Difference Map. Shows the difference in suitability between FR and COR, which was calculated by subtracting COR from the FR predicted suitability (values from 0 to 1). Subpanels (AC) show different maps at higher resolution, with (AC) corresponding with inset boxes. Areas that had higher suitability in FR than COR have positive values and are in red, orange, and yellow, whereas areas that had higher suitability in COR than FR have negative values and are in light and dark blue. Areas with the largest drop in suitability from FR to COR are Big Sur, the Santa Barbara Coast, and the Channel Islands, which all have low or relatively low human population density. Meanwhile, the suitability of COR increased in the urban areas of Las Vegas, Tucson, Phoenix, and Sacramento, as well as Stockton, Modesto, Fresno, and Bakersfield in the Central Valley of California. Presence points are in black.
Figure 4. Difference Map. Shows the difference in suitability between FR and COR, which was calculated by subtracting COR from the FR predicted suitability (values from 0 to 1). Subpanels (AC) show different maps at higher resolution, with (AC) corresponding with inset boxes. Areas that had higher suitability in FR than COR have positive values and are in red, orange, and yellow, whereas areas that had higher suitability in COR than FR have negative values and are in light and dark blue. Areas with the largest drop in suitability from FR to COR are Big Sur, the Santa Barbara Coast, and the Channel Islands, which all have low or relatively low human population density. Meanwhile, the suitability of COR increased in the urban areas of Las Vegas, Tucson, Phoenix, and Sacramento, as well as Stockton, Modesto, Fresno, and Bakersfield in the Central Valley of California. Presence points are in black.
Diversity 16 00640 g004
Figure 5. Univariate response curves for each of the top predictors in the FR model. The predicted value represents the relative probability of presence (0 to 1) associated with the value of the given predictor. VPD = vapor pressure deficit.
Figure 5. Univariate response curves for each of the top predictors in the FR model. The predicted value represents the relative probability of presence (0 to 1) associated with the value of the given predictor. VPD = vapor pressure deficit.
Diversity 16 00640 g005
Figure 6. Univariate response curves for each of the top predictors in the COR LT model. The LQT model had nearly identical results (see Supplemental Material Figure S4). The predicted value represents the relative probability of the presence of monarch butterflies associated with the value of the given predictor.
Figure 6. Univariate response curves for each of the top predictors in the COR LT model. The LQT model had nearly identical results (see Supplemental Material Figure S4). The predicted value represents the relative probability of the presence of monarch butterflies associated with the value of the given predictor.
Diversity 16 00640 g006
Table 2. Final predictors used in the FR model, along with each predictor’s percent contributions and jackknife training gain.
Table 2. Final predictors used in the FR model, along with each predictor’s percent contributions and jackknife training gain.
PredictorPercent ContributionJacknife Training Gain
Minimum Temperature (December)83.71.4192
Cloud Cover (December)9.10.4655
Precipitation Total (October to February)4.40.2271
Vapor Pressure Deficit (December)2.00.5341
Wind (February)0.80.1435
Table 3. Final predictors used in the COR model, along with each predictor’s percent contributions and jackknife training gain. Values are reported from both the best-performing models (LT and LQT). The variables are ordered based on the relative percent contribution to the LT model result, which differs slightly from the LQT. The numerical values for training gain are nearly identical for each predictor regardless of the model. * Only 4 km resolution was used for human footprint.
Table 3. Final predictors used in the COR model, along with each predictor’s percent contributions and jackknife training gain. Values are reported from both the best-performing models (LT and LQT). The variables are ordered based on the relative percent contribution to the LT model result, which differs slightly from the LQT. The numerical values for training gain are nearly identical for each predictor regardless of the model. * Only 4 km resolution was used for human footprint.
PredictorOptimal ScalePercent Contribution (LT)Jacknife Training Gain (LT)Percent Contribution (LQT)Jacknife Training Gain (LQT)
FR Logistic Output1 km61.31.600163.31.6003
% Cover Medium Development1 km31.81.299530.11.3045
Human Footprint Level4 km *1.80.73340.20.7334
% Cover Low Development1 km1.80.82571.50.8291
% Cover Open Water4 km0.90.5240.20.5342
% Cover Shrubs4 km0.80.35752.10.3587
% Cover Cultivated Land4 km0.50.65351.10.6535
% Cover Wetlands4 km0.30.21770.20.2177
Topographic Rugosity Index4 km0.30.23130.30.2313
% Tree Cover1 km0.20.57690.60.5769
Topographic Position Index1 km0.10.09060.10.0906
Stream Density (Excluding Ephemeral Water Sources)4 km0.10.198800.1988
Stream Density4 km00.2220.10.222
% Cover High Development1 km00.994600.9946
% Cover Barren Land4 km00.221100.2211
% Cover Unclassified Land4 km0000
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Fisher, A.R.; Bean, W.T.; Villablanca, F.X. A Multi-Scale Species Distribution Model for Migrating and Overwintering Western Monarch Butterflies: Climate Is the Best Predictor. Diversity 2024, 16, 640. https://doi.org/10.3390/d16100640

AMA Style

Fisher AR, Bean WT, Villablanca FX. A Multi-Scale Species Distribution Model for Migrating and Overwintering Western Monarch Butterflies: Climate Is the Best Predictor. Diversity. 2024; 16(10):640. https://doi.org/10.3390/d16100640

Chicago/Turabian Style

Fisher, Ashley R., William T. Bean, and Francis X. Villablanca. 2024. "A Multi-Scale Species Distribution Model for Migrating and Overwintering Western Monarch Butterflies: Climate Is the Best Predictor" Diversity 16, no. 10: 640. https://doi.org/10.3390/d16100640

APA Style

Fisher, A. R., Bean, W. T., & Villablanca, F. X. (2024). A Multi-Scale Species Distribution Model for Migrating and Overwintering Western Monarch Butterflies: Climate Is the Best Predictor. Diversity, 16(10), 640. https://doi.org/10.3390/d16100640

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop