Biotic predictors with phenological information improve range estimates for migrating monarch butterflies in Mexico

341 –––––––––––––––––––––––––––––––––––––––– © 2019 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Tim Bonebrake Editor-in-Chief: Hanna Tuomisto Accepted 1 November 2019 43: 341–352, 2020 doi: 10.1111/ecog.04886 doi: 10.1111/ecog.04886 43 341–352


Introduction
Coarse-resolution species distribution models (SDMs), which estimate the environmental response of a species based on occurrence data and gridded predictor variables, have classically been fit solely with abiotic variables (e.g. temperature, precipitation). This has been due in a large part to prevailing thought that biotic variables are relevant to species distributions only at local extents (Prinzing et al. 2002, Pearson andDawson 2003), adding mere random noise at larger extents (i.e. the Eltonian Noise Hypothesis sensu Soberón and Nakamura 2009), and that SDMs are expected to result in accurate estimates without the inclusion of biotic variables (Pearson and Dawson 2003). Additionally, the strength of biotic interactions has been shown in some cases to vary considerably over space and time (Thompson 1999, Meier et al. 2011, potentially complicating the underlying assumption that they are stationary for most SDM applications (Wisz et al. 2013). Biotic variables have thus been traditionally relegated to more local-scale applications such as population demographic models (Peterson et al. 2011). However, as both abiotic and biotic factors often contribute to how a species' distribution is determined, abiotic variables alone may not provide the information sufficient for an accurate estimation of a species' distribution at coarser scales. Multiple studies have demonstrated that biotic predictor variables related to the ecology of the focal species have led to increased model accuracy (Araújo and Luoto 2007, Pellissier et al. 2010, Wisz et al. 2013, De Araújo et al. 2014, Belmaker et al. 2015. After all, the spatial patterns of species occurrence data are not only a product of physiological constraints linked to climate or topography, but also of local biotic interactions that govern the probability of presence at a site. Further, failing to consider biotic interactions can result in incorrect predictions of species' distributions under climate change scenarios (Gilman et al. 2010). Although there are many potential methods for integrating interacting species into single-species SDMs (Wisz et al. 2013), a common methodology is including them as predictor variables, whether as occurrence localities or as models of their distributions (Araújo and Luoto 2007, Heikkinen et al. 2007, Bateman et al. 2012, Freeman and Mason 2015, Lemoine 2015. Joint species distribution models can also account for species interactions, as well as separate shared environmental responses from the effects of co-occurrence (Warton et al. 2015), but for studies of single species with unidirectional interactions, using predictor variables based on interacting species reduces the possibility of confounding shared responses with co-occurrence when clear interactions have already been established. This methodology most often assumes that interactions are stationary over space and time, but for studies over broad regions or time scales, it may be important to examine how interactions may change with the environment (Wisz et al. 2013).
In the case that the focal species' distribution is hypothesized to be affected by interactions with multiple species, estimates of species richness can be employed as biotic predictor variables (Koenig and Haydock 1999). Species distribution model predictions of individual interacting species, which should have an equivalent type of interaction with the focal species (i.e. providing a resource), can be combined (or 'stacked') in various ways to estimate richness (Ferrier and Guisan 2006). Other variables such as abundance or total biomass of interacting species may be more closely related to the focal species' biology, but as data on both species interactions and abundance are often lacking (Eltonian and Prestonian shortfalls, respectively; Hortal et al. 2015), estimated richness may be a suitable proxy for some systems. An important caveat is that if populations of the interacting species are affected by the presence of the focal species, a key assumption in single-species SDMs is violated, namely that all predictor variables are 'unlinked' to (e.g. not affected by) the focal species (scenopoetic variables sensu Hutchinson 1978, Peterson et al. 2011. Thus, in order to use stacked SDMs as predictor variables in single-species models, they must be derived from interacting species unaffected by (or 'unlinked' to) fluctuations in the population level of the focal species (Anderson 2017). For example, generalist nectar-feeders and the plants on which they feed (resource plants) can represent such a system, as generalists are likely to have only weak links with any one resource plant species (Anderson 2017). In this case, flowering phenology should be related to the strength of these unidirectional interactions, as only plants in flower can provide resources. Other SDM studies that have used information about resource plants in SDMs for nectar-feeders have focused on specialists, which may violate the 'unlinked' assumption, and they do not consider flowering phenology (Giannini et al. 2013, Silva et al. 2014. Here, we focus on a wide-ranging butterfly that is a generalist nectarivore (at the adult life stage) and its associations with broad assemblages of resource plants with which it has unidirectional interactions during its migration.
The monarch butterfly (Danaus plexippus L. 1758, hereafter 'monarch') is a charismatic and wide-ranging insect that has become established as a conservation icon (Gustafsson et al. 2015). The eastern North American monarch population has a spectacular multi-generational annual migration, traveling from breeding grounds in the U.S. and Canada to small overwintering areas in the highlands of central Mexico (Oberhauser et al. 2015). In recent years, the total abundance of overwintering colonies in Mexico has seen several precipitous drops, due to a combination of deforestation at Mexican overwintering sites, severe weather events, and loss of milkweed in the U.S. (Brower et al. 2012). Developing monarch larvae are obligate to milkweed species (Asclepias spp.), and those species are an important nectar resource for generalist adults. The increasing use of herbicide-resistant crops coupled with an increase in herbicide application in the midwestern U.S. has resulted in widespread reduction in monarch breeding habitat and is widely acknowledged to be a major driver of this decline (Pleasants and Oberhauser 2013). Monarchs migrating through Mexico to and from the overwintering sites also rely on a variety of tree species that provide structure for roosting and suitable microclimates for shelter (Brower et al. 2009, Howard andDavis 2009), which have been impacted by deforestation (Sáenz-Romero et al. 2012). Great interest in monarch conservation has generated a considerable amount of research, resulting in the establishment of the trinational North American Monarch Conservation Plan (Oberhauser et al. 2008) and citizen science monitoring datasets (Howard and Davis 2015).
Monarchs exhibit three discrete behaviors during the multi-generational migration -breeding, migrating and overwintering -and range estimations have only been made in the breeding (Batalden et al. 2007, Flockhart et al. 2013, Lemoine 2015 and overwintering Peterson 2003, Sáenz-Romero et al. 2012) ranges. Although it represents a significant portion of the full range, the migration route through Mexico has remained understudied. During the fall migration, monarchs migrate south from the plains of southern Texas and northeastern Mexico along the eastern slopes of the Sierra Madre Oriental to overwinter in highland conifer forests in the central state of Michoacán (Oberhauser et al. 2015). In the spring, they take the same route back north to breeding grounds in the U.S. and Canada. Due to a scarcity of monarch occurrence data for the Mexican migration route, however, distributional estimates for both the fall and spring are lacking. Fortunately, a new and systematic sampling effort by the National Commission of Protected Natural Areas of Mexico (CONANP) is beginning to close the occurrence data gap between the breeding and overwintering ranges (Botello et al. unpubl.). As the range limits of migrating monarchs in Mexico currently remain unresolved, and the urgency for accurate distributional estimates for monarchs is particularly acute given the recent population declines, we built the first SDMs for monarchs actively migrating through Mexico using occurrence data from CONANP and citizen science datasets.
An important question is whether migrating monarchs are close to equilibrium with their environment, as most SDM implementations assume equilibrium between the focal species and the environmental predictor variables employed (Guisan and Thuiller 2005). For a migrating species moving through a landscape, a potential concern is that presence localities may not adequately represent suitable environmental conditions, as the species has relatively transient relationships with these environments. In the case of the monarch, however, there are multiple reasons why this might not be an issue. The monarch seems to possess a distinct migratory niche differentiated from those associated with its other discrete behaviors (Batalden et al. 2007), which is evidence that it at least partially tracks suitable environmental conditions during the migration. Further, as the monarch makes annual visits to similar areas during the migration season, it is likely that monarchs have been exposed to the suite of environments available in this region over evolutionary time (in contrast to newly colonized invasive species; Jiménez-Valverde et al. 2011). Hence, the possibility of environmental equilibrium for areas seasonally occupied by a migratory species with unique migratory niche characteristics seems reasonable for designing a SDM methodology for this species, as well as beginning with the assumption that interactions are stationary for this migration region and duration. This line of thinking underlies more extensive research done on the seasonal distributions of migratory birds (Martínez-Meyer et al. 2004, Nakazawa et al. 2004, Marini et al. 2010, Laube et al. 2015. Specifically, it follows that monarch SDMs trained on monthly environmental data that span the time period from the initiation of migratory behavior until the transition to overwintering behavior should capture changing seasonal niche characteristics better than models that ignore seasonal differences. With dual goals of advancing techniques to integrate biotic interactions into single-species SDMs and pioneering distributional estimates for migrating monarchs, we used Maxent ) to build SDMs for each month of the migration through Mexico and compared them to determine whether abiotic, biotic or a combined set of variables were optimal ( Fig. 1), and also examined whether accounting for flowering phenology made a difference. We created biotic predictor variables for migrating monarchs based on two estimates of species richness: plants that provide a nectar resource (Asclepias spp.) and highland trees that provide shelter and roosting sites. Although migrating monarchs feed on other nectar plants besides Asclepias spp., we lacked reliable phenology data for other species in Mexico, and thus used Asclepias richness as a proxy for nectar plant richness in this study. In a novel implementation for SDMs, for each month we compared between models using estimated flowering plant richness that accounted for phenology and those that did not. We selected optimal monarch SDMs using information criteria and additionally determined if they performed better than random on spatially withheld data by extending a recently described null SDM approach. We hypothesized that 1) high richness of nectar plants and roosting trees should indicate suitable biotic conditions for migrating monarchs, 2) a combination of abiotic and biotic variables should result in the most accurate monarch SDMs for each month and 3) among these models, those which accounted for phenology of nectar plants should perform the best.

Species occurrence and phenology data
We acquired occurrence data for monarchs from both a monitoring program and a long-term citizen science dataset. We used CONANP monitoring data for the fall season of 2015 (currently proprietary; Commission for Environmental Cooperation 2017), and the Journey North citizen science dataset for the years 2000-2015 (<https://journeynorth. org/>; accessed March 2018). These datasets were combined, and occurrence locality coordinates were extracted for the months of September (n = 168), October (n = 1122) and November (n = 1301), which encompass the period of peak fall migration through Mexico. To reduce artifactual clustering of records due to biased sampling, we spatially thinned the occurrence localities by 10 km for each month (Fig. 2), resulting in more even sample sizes across months (September: n = 63, October: n = 188, November: n = 117; Aiello-Lammens et al. 2015). We did not include models for August (n = 31) and December (n = 26) as the occurrence data were heavily clustered, the predictive ability of preliminary models was very poor, and very few data from CONANP were available for these months.
In order to construct the biotic variables, we selected species of both nectar plants and roosting trees found in the southwestern U.S. and the migration range in Mexico based on available natural history information demonstrating associations between these species and migrating monarchs. We thus acquired occurrence records of 25 species of Asclepias and 27 species of trees across different genera (Abies, Carya, Cupressus, Juglans, Juniperus, Pinus, Quercus and Taxodium) from the Global Biodiversity Information Facility (GBIF, downloaded August 2014). We removed occurrence records that lacked coordinate information and those we determined had erroneous coordinates (i.e. likely spatial outliers) or were incorrectly classified (taxonomic misidentification) based on expert opinion of these species' ranges and identifications. This process included checking questionable occurrence records against herbarium specimens primarily from the Instituto de Biología, UNAM, which were also used to derive flowering months for Asclepias spp. We also relied on a number of other institutional collections and databases: Escuela Nacional de Ciencias Biológicas, Museo Metropolitano de Monterrey, Inst. de Ecología, A.C., Missouri Botanical Garden, and the Tropicos electronic database (Missouri Botanical Garden 2018). To reduce artifactual clustering as with the monarch data, we then spatially thinned each plant species' occurrence records by the same distance, resulting in reductions ranging from 28% to 63% per species (see Supplementary material Appendix 2 Table A1 for sample sizes before and after thinning, and Supplementary material Appendix 1 A1 for methodological details). We retained for analysis only 20 Asclepias species and 24 tree species with 15 or more occurrences in Mexico. In addition to occurrence data, we derived flowering times by month for each species of Asclepias based on phenological information from museum specimens to create a plant phenology database.

Abiotic predictor variables
We acquired abiotic environmental rasters at 30 arc-second resolution (~1 km at the equator) at different temporal resolutions. From the WorldClim 2.0 database, we downloaded long-term  monthly and annual averages of temperature, precipitation and solar radiation (Fick and Hijmans 2017). Although the monarch occurrence data was from 2000 to 2015, we considered that the broad patterns represented by these long-term datasets could still accurately represent environmental conditions for these monarch data. For the plant models, we considered long-term annual averages: 19 bioclimatic variables based on temperature and precipitation, as well as the mean and coefficient of variation of annual solar radiation. For the monarch models, we considered long-term averages for each month (September-November): mean temperature, maximum temperature, minimum temperature and mean precipitation. As we had 345 little a priori knowledge of the climatic variables most associated with either the distributions of the focal plant species or for monarchs throughout the migration route in Mexico, we opted to build models with all considered variables instead of choosing which ones we thought to be most appropriate. Further, as we were solely interested in the predictive accuracy of models and did not transfer models to other times or places, instead of removing correlated variables before modeling, we followed the machine learning approach and allowed the Maxent algorithm to remove those with poor explanatory power via regularization through a tuning exercise Dudík 2008, Merow et al. 2013).

Plant richness estimates
We built SDMs with Maxent 3.4.1 ) for all plant species and stacked them to make richness estimates that became biotic predictor variables for the monarch SDMs. We first built candidate models with the cleaned and spatially thinned occurrence datasets using different combinations of simple and complex model settings. We then selected optimal models based on spatial cross-validation performance by sequentially considering measures of omission rate (10 percentile) first and model discrimination (AUC test ) second (Muscarella et al. 2014; Supplementary material Appendix 1 A2). Finally, we projected these models onto geographic space and summed their continuous predicted suitability values (cloglog transformation) to make 'richness' estimates for roosting trees and Asclepias spp. (Calabrese et al. 2014; Supplementary material Appendix 1 A3). These richness estimates were considered stationary over time, but we additionally made monthly estimates for Asclepias spp. by including only those species that were flowering in that month according to our phenology database. In this study, flowering times per species were assumed to be the same across space for the study extent. This methodology allowed us to compare the performance of models that considered nectar plant phenology (using the monthly richness estimate) with those that did not (using the stationary estimate). Our expectation was that predicted monarch suitability would have a positive relationship with estimated richness for both plant classes, and that the best models would include biotic variables that account for phenology. We also conducted an analysis with low-abundance Asclepias species removed to check if this would bias richness estimates (Supplementary material Appendix 1 A4).

Monarch SDMs
Next, we built SDMs with Maxent for monarchs for each month of the fall migration using different combinations of abiotic and biotic variables (Table 1). We followed the data preparation and modeling methodologies for plants  (Supplementary material Appendix 1 A1, A2), but selected models based on information criteria (AICc) instead of cross validation (for details on this decision, see Supplementary material Appendix 1 A5). Thus, we built abiotic (monthly climate), biotic (estimated plant richness) and combined models, where biotic and combined were either stationary (all Asclepias spp.) or monthly (flowering Asclepias spp.; Fig. 1). For each combined model, we documented permutation importance of variables and plotted response curves. Maxent derives permutation importance by randomly changing the values of each variable in turn and calculating the resulting decrease in the training AUC -large decreases are interpreted to mean that the model depends strongly on the variable (Phillips 2017). As SDMs selected by minimizing AICc may not perform well on withheld data (Galante et al. 2018), we determined whether these models also performed well for cross validation. We extended a recently developed SDM null model approach to determine if the same measures of discrimination and omission used to select optimal plant SDMs (AUC test and OR10, respectively) were higher than random for selected monarch SDM settings when evaluated on spatially withheld occurrence data (Bohl et al. 2019; Supplementary material Appendix 1 A6). See Supplementary material Appendix 3 for the function used to create the null SDMs with spatial block partitioning.

Plants
Optimal model settings for plants differed considerably among species, spanning the full range of potential model complexity considered (Supplementary material Appendix 2 Table  A2). For details on plant SDM results, see Supplementary material Appendix 1 A7. Estimated species richness for roosting trees was greatest in the high elevation region of the Trans-Mexican Volcanic Belt that includes the overwintering sites, with some moderate predictions farther north into the Sierra Madre Oriental (Fig. 3). Estimated richness for flowering Asclepias was also highest in the Volcanic Belt region, but high predictions extended north of Michoácan to the Central Plateau (Fig. 3). Fewer species of Asclepias flowered as winter approached: 12 species for September, dropping to 7 for October and 4 for November. The estimated richness patterns across space for September and October were similar to the stationary Asclepias variable, with moderate estimated richness occurring in gaps north into the Central Plateau. In contrast, November had the highest estimated richness at the southern tip of the Sierra Madre Oriental, with moderate values in the Volcanic Belt and towards the Gulf Coast.

Monarchs
For monarchs, the model settings chosen as optimal by AICc for sets of variables by month ranged from simple (fewer feature classes, higher regularization, fewer parameters) to complex (more feature classes, lower regularization, more parameters). Optimal settings were uniformly simple for September, but covered a wide range of complexity for October and November (Table 2). Overall, the simplest models were consistently the biotic-only ones, whereas the most complex models were the combined ones, those months with higher sample size, and those variable combinations that led to models with a higher number of predictor variables. Removing low-abundance Asclepias species led to very similar results overall (Supplementary material Appendix 1 A4, Appendix 2 Table A3). Model selection among sets of variables by AICc resulted in the combined, phenologyinformed monthly model for all three months ( Table 2). The permutation importance of the biotic variables for these combined monthly models was as high or higher than that of the abiotic ones, with the Asclepias estimated richness receiving a very high percentage for September and November (Supplementary material Appendix 2 Fig. A4).
The combined monthly models varied in both complexity and geographic prediction. The September model used only the variables mean precipitation, maximum temperature, Asclepias and trees, and had smooth and simple response curves (Supplementary material Appendix 2 Fig. A1). For September, monarch suitability was highest in drier areas without high temperature extremes, and also in areas with predicted richness that was low for Asclepias spp. but high for roosting trees. This corresponded mainly to the northern plains and northern Sierra Madre Oriental, but also to spots within the southern Volcanic Belt region (Fig. 4). The October model was more complex, using all input variables, and had responses curves that were somewhat more jagged. Predicted monarch suitability in October was highest in drier areas with lower temperature extremes, and also in areas with predicted richness that was low for Asclepias spp. but high for roosting trees (Supplementary material Appendix 2 Fig. A2), corresponding mainly to the northern Sierra Madre Oriental but also to areas south in the Central Plateau and Volcanic Belt region. The November combined monthly model was even more complex than that of October and also used all input variables. Although monarch suitability in November was highest under similar abiotic conditions to October (although now with a negative relationship to minimum temperature), it was positively related to predicted richness for both Asclepias and roosting trees (Supplementary material Appendix 2 Fig. A3). High suitability was predicted farther south than previous months and focused in the southern Sierra Madre Oriental and Volcanic Belt region.
Most real combined models selected via AICc performed consistently better than real abiotic or biotic models when compared with null distributions for discriminatory ability (AUC test ), with results either close to or exceeding the significance level at α = 0.05, but the same was not true for omission rate (OR10; Supplementary material Appendix 2 Fig. A5, A6). Null distributions of both statistics had a broad range of variation (Supplementary material Appendix 2 Table A4). For more details on the null modeling results, see Supplementary material Appendix 1 A8.

Combined models with flowering phenology performed best
As occurrence data for species of concern become increasingly available, it is crucial that SDMs be fit using predictor variables that reflect not only physiological constraints, but also biotic interactions that are important in shaping their distributions (Wisz et al. 2013). The results of this study add to the growing evidence that biotic predictor variables are important across large spatial extents. From earlier studies demonstrating the advantages of biotic variables at coarse resolutions on butterflies and their host plants (Araújo and Luoto 2007) and owls and their facilitators (Heikkinen et al. 2007), a wealth of studies focusing on different interactions have emerged proposing that biotic interactions be considered more in SDMs ( Van der Putten et al. 2010, Belmaker et al. 2015. Although many examples on animals and food/host plants have focused on specialist-host associations (Giannini et al. 2013, Raath et al. 2018, as these interactions should be particularly strong, a number have shown that predicted ranges for generalists are also improved with biotic predictors (Araújo and Luoto 2007, De Araújo et al. Figure 3. Estimated species richness from stacked SDMs for trees and Asclepias spp. The top row shows the stationary stacks that include all species for each plant class, while the bottom row shows the monthly stacks for Asclepias with species chosen by phenology. The monthly stacks have that month's monarch occurrence localities overlaid to show the relationships between predicted flowering plant richness and observed monarch occurrences. 2014, Lemoine 2015). Supporting these findings, we found that a combination of predictor variables based on both climate and richness estimates of interacting species can lead to improved SDMs for migrating monarchs with generalist associations. The combined models had the best performance every month, and those that accounted for flowering phenology (combined monthly) had better performance than those that did not (combined stationary), although they had the same number of predictor variables. The importance of biotic variables in combined models was high (often higher than abiotic variables) regardless of whether or not flowering phenology was considered. Lastly, the combined monthly models performed best (though not always significantly so) compared with null distributions for discriminatory ability.
Importantly, the combined monthly models featured in this study produce suitability predictions over space that are realistic when considering the monarch's migration route through Mexico. Predicted monarch suitability early in the migration season in September was highest in the north with lower values on the eastern slopes of the Sierra Madre Oriental where migrating monarch abundance should be relatively low. High suitability was predicted, however, farther south in the Volcanic Belt region where few monarch occurrences are found yet estimated richness was high for both classes of interacting plants. In October, most monarchs are moving further south and expanding westward into highland conifer forests that are especially dense near the overwintering grounds (the asterisk in Fig. 4), which is reflected in the model prediction. In November towards the end of the migration, monarch suitability dropped in the northern plains and was highest where most monarchs are beginning to settle at overwintering sites in conifer forests of the Volcanic Belt region. Thus, we think these maps can contribute to management efforts for monarch monitoring and conservation (see 'Conservation implications'). It is also important to note that our results agree with those of Lemoine (2015), who found that the inclusion of a biotic variable derived from nectar plant occurrence data improved range estimates for monarchs in the northern breeding range. Our study presents some methodological advancements, such as using predictors based on the estimated richness of different classes of interactors, as well as incorporating phenological information into biotic predictors, which can be used to produce future range estimates for monarchs throughout their range.

Interpreting model responses of biotic variables
Although the combined monthly models performed best by information criteria, the directionality of model responses for the biotic variables proved difficult to interpret for some months, as they did not align with our expectations regarding monarch ecology. It is important to inspect model response curves to ensure they match ecological expectations (Guevara et al. 2018), especially when using biotic variables chosen based on natural history information. To our knowledge, no studies on biotic predictors and SDMs to date have focused on the directionality of model responses. This is crucial because responses that are not ecologically realistic would compound prediction errors with model transfer to other times and places, particularly for future predictions as climate change can have complex effects on existing interactions (Blois et al. 2013). We expected positive model responses for both biotic variables over all months, as monarchs depend on both classes of plants during the migration. Although model responses were positive for roosting trees for all three months and for Asclepias spp. for November, the Asclepias responses for September and October were negative. The magnitude of this difference is likely exacerbated because of the higher variable importance for Asclepias estimated richness. As both stationary and monthly combined models had these relationships, the differences in estimated richness of November's monthly Asclepias variable compared with the other months cannot explain this discrepancy. Table 2. Maxent SDM settings chosen as optimal for monarchs by month and variable set with associated performance statistics. Models with the lowest AICc score for the month are highlighted in light gray. Settings shown are feature class combinations (features) and regularization multiplier (rm). Statistics shown are AUC calculated on training data (AUC train ) and averaged over testing data (AUC test ), omission rates for 10 percentile (OR10) training values, delta AICc (based on the lowest AICc among optimal models), and the number of non-zero coefficients (nparam).  Compared with the breeding and overwintering ranges, the associations with climate and resources in the migratory range when monarchs are very mobile are likely weaker with higher uncertainty. However, several possibilities may explain these results. The November response could be positive because a greater proportion of the total November monarch occurrences were located farther south in the Volcanic Belt region where estimated stationary and monthly Asclepias richness was highest. Alternatively, a variable missing from the model with a clear ecological explanation could be negatively correlated with estimated Asclepias richness for September and October. It is also possible that Asclepias estimated richness only has a positive response at the end of the migration when most monarchs have reached the overwintering areas, become more sedentary, and engage more in targeted selection of sites with a diversity of nectar plants (with the negative modeled responses for the other two months reflecting either correlation with a missing variable, or a spurious response related to spatial patterns when the monarchs are more mobile). If true, this last possibility entails that interactions between monarchs and nectar plants would indeed be variable across space and/or time, but field studies would be needed to confirm this.
Revising the biotic variables included in future studies may result in responses more aligned with ecological understanding of this system. Some examples include adding more nectar plant species into richness estimates, using predicted abundance instead of richness (although these data do not currently exist for Mexico), testing other combinations of biotic variables (for example, nectar-providing and shelter variables both separately and in combination), exploring finer temporal scales to better account for phenology, or temporal matching of all monarch occurrences with monthly climate data for a single model. Regardless, we interpret that the combined monthly models resulted in reasonable distributional estimates for migrating monarchs in Mexico that can serve as starting points in developing better models for the future.

Conservation implications
These SDMs built with both abiotic and biotic predictor variables help improve distributional estimates for monarchs in the Mexican migration corridor, and thus can serve as important tools for conservation. Prior to the present study, biogeographical information for monarchs in this region has been lacking, but these results can lead to a first approximation of which protected areas are important for the monarch's migratory route through Mexico. As model development progresses with ongoing CONANP data production, these models can be improved over time to benefit prioritization for monitoring sites and habitat conservation efforts. Accurate monarch SDMs in the migration corridor could help local conservation practitioners and managers select areas for mitigating herbicide use (which can decimate nectar food sources) and deforestation (which can reduce habitat for roosting during migration or overwintering). Importantly, under this framework temporal change in predicted monarch distributional area can also be estimated, and SDMs for different years could serve as independent data sets for model evaluation and selection to help improve individualyear model accuracy. Future studies should use these distributional estimates to develop a conservation connectivity model to help optimize conservation efforts, and also investigate the potential effects of climate change to the distributions of monarchs and the plants they use for resources and shelter. Indeed, conditions in the overwintering range are projected to be increasing unsuitable, even lethal, for monarch survival (Barve et al. 2012, Sáenz-Romero et al. 2012, and knowledge of possible climate refuges migrating monarchs may use in the future will be valuable. Certainly, better estimates of the monarch's migratory distribution in Mexico are long overdue (considering the bulk of existing research in the breeding and overwintering ranges), and recent population declines make them all the more urgent. The models presented in this study, and the methodology used to build them, can be foundations for future research that help move monarch conservation forward in Mexico, as well as further strengthen the trinational data-sharing that is essential for the development of high-quality models.