Estimating the annual distribution of monarch butterflies in Canada over 16 years using citizen science data

,


Introduction
Monarch butterflies (Danaus plexippus, Linnaeus, 1758) in North America are known for their longdistance movements among the United States, Canada, and Mexico.Monarchs are comprised of one genetic population with two migratory flyways separated by the Rocky Mountains (Brower 1995;Lyons et al. 2012).The eastern North American population overwinters in the highlands of central Mexico and migrates to breeding areas in the midwestern and eastern United States and Canada over successive breeding generations before a final generation migrates back to Mexico (Brower 1995;Flockhart et al. 2013).The western population overwinters in coastal areas of California and some individuals migrate north over multiple breeding generations to all western states and occasionally reaching very southern portions of Canada in British Columbia before a final generation migrates back to coastal California (Brower 1995;Dingle et al. 2005;James et al. 2018).Both migratory populations have declined over the last several decades because of multiple threats (Brower et al. 2012;Schultz et al. 2017), and monitoring their distribution and abundance is a necessary aspect of their ongoing management (Commission for Environmental Cooperation 2008).
Their small population size and large fluctuations in population growth rate (Brower et al. 2012;Flockhart et al. 2015;Semmens et al. 2016;Schultz et al. 2017) were primary reasons why monarchs were designated a species of 'special concern' in Canada in 1997, a status that was reconfirmed in 2010 (COSEWIC 2010).Monarch population trajectories currently meet the threshold to be considered endangered in Canada and the Committee on the Status of Endangered Wildlife in Canada has recommended this listing to the federal government for a classification decision due in 2019 (COSEWIC 2016).Both North American populations are at risk from multiple factors such as disease, extreme weather, and overwinter habitat loss (Brower et al. 2012).There is increasing evidence to support that population viability is three times more sensitive to breeding habitat loss compared with habitat loss on the wintering grounds (Flockhart et al. 2015) and, therefore, the monarch population decline is driven, at least in part, by decreasing abundance of their obligate host plants, milkweeds (Asclepias spp.) across the United States and Canada (Pleasants and Oberhauser 2013;Flockhart et al. 2015;Stenoien et al. 2016;Oberhauser et al. 2017;Thogmartin et al. 2017a).Prioritizing locations based on the greatest probability of occurrence for breeding habitat restoration in Canada, therefore, relies on spatially explicit mapping the distribution of monarchs across space and time that extends beyond simple generalized maps of distribution.
A variety of long-term citizen science monitoring programs across North America provides data to map the distribution and phenology of many organisms.Monarchs are the focus of multiple citizen science programs that document locations, migration patterns, reproduction, and disease rates (e.g., Journey North, Monarch Alert, Monarch Watch; for a more complete list see, Ries and Oberhauser 2015).Citizen science programs that document monarch occurrence data may record first monarch observations of the season (e.g., Journey North; Davis and Howard 2005) or a compilation of observations throughout the spring and summer for monarchs as well as other North American butterflies (e.g., eButterfly; Prudic et al. 2017).All programs record the date and location where monarchs are observed, giving rise to so-called presence-only data, and have been used to document migration routes (Howard and Davis 2009), migration speed (Davis and Howard 2005), and breeding distributions (Batalden et al. 2007;Flockhart et al. 2013).Presence-only citizen science data from different programs, if not biased to only include specific locations or portions of the annual cycle, can be combined to increase sample sizes to document annual variation in distribution patterns inherent in insects.Understanding annual variation in distribution would provide valuable information to aid in management and conservation planning (Warren et al. 2001).
The annual breeding distribution of monarchs in Canada depends on the colonization of individuals that migrate from the United States (Brower 1995) and the conditions (habitat, physiological, and geographic) under which such migrations occur (Brown et al. 1996).The eastern North American monarch butterfly population often reaches into southern parts of Canada as far west as Alberta (Brower 1995;Flockhart et al. 2013;Flockhart et al. 2019).In some extreme cases, individuals observed early in the breeding season have likely migrated directly from their overwintering grounds in Mexico (Miller et al. 2012).In contrast, monarchs that reach the extreme southern portions of British Columbia are likely from the western North American population, which overwinters in California (Environment and Climate Change Canada 2016; Yang et al. 2016;James et al. 2018).Introgression among the eastern and western populations (Lyons et al. 2012) implies that climatic, physiologic, and geographic factors may consistently influence migration and hence colonization of  (Oberhauser et al. 2001), geographic limits of migration and recruitment (McKinnon et al. 2010), or climatic thresholds for flight in insects.Additionally, these factors may not be mutually exclusive if, for instance, host plant availability may be limited by climate (Lemoine 2015).Because Canada legislates species at risk, including monarchs, independently of the US or Mexico, effective conservation planning requires understanding the extent of the breeding range in Canada, the natural fluctuation from year to year (Brown et al. 1996;Prysby and Oberhauser 2004), and what might be the primary drivers of annual variability.
In this study, we use the citizen science programs eButterfly (Prudic et al. 2017) and Journey North (Davis and Howard 2005) to estimate the annual breeding distribution of monarch butterflies in Canada between 2000 and 2015.In doing so, we test predictions from several hypotheses to explain variation in monarch occurrence (Table 1).Our results provide the geographic context for decision-makers to identifying priority areas to engage in restoration and other conservationrelated activities.Observations of monarchs are biased towards areas with higher human population density.

Monarch butterfly observations
We used data from eButterfly (e-butterfly.org/)and Journey North (learner.org/jnorth/monarchs/),which are online data repositories for citizen scientists to record their butterfly observations across North America (Howard and Davis 2004;Prudic et al. 2017).Citizen scientists occasionally report the same observations to both programs and given that both programs produce presence-only data, these data can be combined to increase sample size to better fit species distribution models.Occurrence records (n = 49 098) were vetted to remove observations that were outside of Canada, prior to 2000 or after 2015, or missing either a full date or the latitude and longitude.Our refined data set contained 22 974 observations.

Species distribution models
Our objective was to develop models to predict the annual distribution of monarchs across Canada across 16 years (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015).For each year, we used likelihood-based models to assess the probability of occurrence of monarch observations based on underlying environmental variables to test monarch breeding distribution determinants (Flockhart et al. 2013).We then used the best-supported model to predict a probability of occurrence for a given year across Canada and calculated the breeding distribution occupied (km 2 ) at several probability thresholds.Although linear trends in parameter estimates from our model may indicate directional change and hence increasing importance of certain factors to explain monarch distribution that could be used to predict future distributions, our goal was to examine past breeding distributions, recognizing that colonization of the breeding ground in Canada is an annual event and that 16 years is a short time series to examine directional change in breeding distributions under global change.
We followed the approach of Flockhart et al. (2013), who estimated monthly breeding distributions across eastern North America testing among three non-mutually exclusive hypotheses to explain monarch occurrence (Table 1).If breeding distribution is dependent on habitat suitability (Oberhauser et al. 2001), then monarch occurrence should be associated with percent tree cover and the normalized difference vegetation index (NDVI).Although monarchs typically occur in the eastern Great Plains and eastern broadleaf forest ecoregions, depending on the year and monarch relative abundance, observations may be positively (e.g., if a larger populations occur in the eastern broadleaf forests region) or negatively (e.g., if a larger populations occur in the Great Plains region) correlated with percent tree cover.Monarch occurrence should be positively related to NDVI, which is an index of green vegetation.If monarch distribution is determined by the geographic limits of migration then monarch occurrence should be associated with latitude, longitude, elevation, and slope (Cockrell et al. 1993;Batalden et al. 2007;Zipkin et al. 2012;Flockhart et al. 2013).Latitude and longitude relationships could either be linear (e.g., latitude) or quadratic (e.g., latitude + latitude 2 ) forms, whereas elevation and slope are predicted to be linear.Finally, monarch occurrence could be dependent on physiological constraints on development and movement imposed by weather conditions such as temperature and precipitation (Zalucki 1982;Cockrell et al. 1993;Batalden et al. 2007;Flockhart et al. 2013).Physiological constraints of monarch occurrence could manifest in linear or quadratic relationships with mean temperature and precipitation.As monarch observations reported by citizen scientists are not from a random or systemic sampling design (Yackulic et al. 2013), for all models we included human population density as an explanatory variable to control for sample selection bias, as citizen scientist observations were concentrated in urban centers with large numbers of people (Flockhart et al. 2013).
Environmental variables were available at fine resolution but the small sample size of monarch observations in some years, and the geographic extent of those observations, limited the resolution of the environmental variables that we could consider.Small sample sizes of observations may result in large sampling variance of environmental variables when using high-resolution spatial data.To ensure that true effects of environmental variables explaining monarch distribution were not obscured under these conditions, spatial raster layers for all explanatory variables were downloaded, aligned, and resampled (resolution = 0.5 degrees) at the extent of Canada (Table 1).The explanatory variables were: NDVI, percent tree cover, mean temperature, precipitation, latitude, longitude, elevation, slope, and human population density (log x+1 transformed).NDVI, temperature, and precipitation were calculated at monthly time intervals; percent tree cover and human population density were calculated at annual time intervals; and latitude, longitude, elevation, and slope were constant over the study period.As our analysis was conducted on an annual interval, explanatory variables at monthly intervals were weighted-average to an annual surface based on the sample size of occurrences in each month of the focal year.

Statistical analysis
We used the maxlike package (function maxlike) in program R to conduct our analysis (Royle et al. 2012;R Core Team 2014).We used Akaike Information Criterion (AIC) values to conduct model selection in two stages.In the first stage, for each year, we tested a global model of all additive variables and all combinations of reduced-parameter models in each hypothesis.The model with the lowest AIC value for each hypothesis designated the model that we carried forward to the second stage of analysis.In total for each year, we compared 22 models in the first stage.In the second stage, we conducted a similar model comparison exercise where we considered a global model of three hypotheses (all variables from the three hypotheses carried forward from the first stage) and thereafter all reduced-model combinations.In total for each year, we compared seven models in the second stage.The model from the second stage with the lowest AIC value was considered the best explanatory model for a specific year; model comparisons from the second stage are presented in Table S1.Once we had the best-supported model, we then used the predict function to estimate the probability of occurrence across the 5.77 million km 2 study area which included all areas of Canada south of 60°N.We used the annual mean human population density as a constant value across the study area when making predictions of breeding distribution to smooth for the effect of occurrence records being more commonly made in areas with high human population density but which does not reflect distribution patterns of monarch butterflies (Flockhart et al. 2013).
The predicted probability of occurrence ranged between zero and one so we considered three probability thresholds to present our results: 0.5, 0.25, and 0.05.The value represents the lowest probability value considered.For example, the breeding distribution at the 0.25 level is the area of all pixels with a probability of occurrence of 0.25 or higher.We calculated the breeding distribution area for each year for western and eastern Canada at the Alberta-British Columbia border that represents the continental divide that separates the eastern and western populations in Canada.We then correlated the annual breeding area estimates with the annual western (Schultz et al. 2017) and eastern (Rendón-Salinas et al. 2018) overwintering sizes.
Overall, all variables considered were supported in the top model in at least one year (Table 3) and model selection uncertainty was low in all years except in 2000, 2012, and 2015 (Table S1).
Monarch distribution was best explained by geographic variables followed by habitat variables (Table 3).Latitude (16 of 16 years) and longitude (15 of 16 years) were the most common variables found in top models, whereas the least common variables supported included linear (3 of 16 years) and quadratic (1 of 16 years) effects of mean temperature and a quadratic effect of precipitation (2 of 16 years).Parameter estimates for human population density were positive in all models, which support the notion that human presence needs to be controlled for when using citizen science data to predict monarch distributions (Table 3; Fig. 3d).
There was no correspondence between the preceding or forthcoming overwintering population estimates in California and Mexico and the estimated breeding distribution in western (Fig. S2a) and eastern Canada (Fig. S2b).However, there was a relationship between a larger estimated breeding distribution and the number of observations of monarchs in Canada (Fig. S3).

Discussion
The core breeding distribution of monarchs in Canada includes portions of southern Ontario (south of Sudbury and Ottawa to the USA border), Quebec (south of Québec City including Montreal and Sherbrooke), and Manitoba (south of Winnipeg).In some years, monarchs were predicted to range over more than a million square kilometers that included southern portions of all provinces.Monarchs in eastern Canada occupied an area that was an order of magnitude larger than the area occupied by monarchs in western Canada.Patterns of occurrence in Canada were best described by geographic limits that are presumably related to limits on northward migration in the late summer or limits in the availability of milkweeds (Asclepias spp.), the obligate host plants (Lemoine 2015).Given the inter-annual variation in breeding distribution, all portions of southern Canada should be considered when developing long-term management strategies for this species; however, the population in eastern Canada occurs each year over a much larger area.
The distribution of monarchs was best explained by geographic limits of migration described with relationships between geographic variables of latitude, longitude, elevation, and slope.These are fixed attributes compared with regional weather and vegetation patterns that vary annually or seasonally and are known predictors of monarch abundance and distribution (Zipkin et al. 2012;Flockhart et al. 2013;Flockhart et al. 2017;Saunders et al. 2018).Therefore, variation in migration timing and geographic distribution in Canada could be dependent on environmental conditions in the United States experienced by monarchs during earlier portions of the annual cycle (Zipkin et al. 2012;Saunders et al. 2018).Monarch butterfly development and survival is known to be impacted by local temperature and precipitation.Although we did not see their effect at the broad spatial scale of our modeling efforts, both temperature and precipitation may be important within the breeding population and warrant future research at a finer geographic scale than performed in this study.For monarchs in eastern Canada that migrate north in the spring, the first observations that occur in southern Ontario are of individuals moving north through Michigan and Ohio that largely originate from the central and southern United States (Miller et al. 2012).In some years, first observations may occur in late April and increase over the next two months suggesting colonization from areas  farther south until approximately July when monarchs have reached their breeding distribution limit (Flockhart et al. 2013).In June, the first observations also occur in southern Manitoba and these butterflies likely originate from the US Midwest via Minnesota and North Dakota.Therefore, while monarchs colonize Canada at two different locations the origins of these individuals also occur from the midwestern and southern regions of the United States where regional weather and habitat conditions will have influenced population processes.In western Canada, monarchs occur in the Fraser Note: To account for the observed higher probability of occurrence in areas with high human density (urban areas) we controlled for human population density (Human) by using the mean human population density across the study area when deriving predicted breeding occurrences for each year.Locations for habitat restoration should be prioritized based on where monarchs have a high probability of occurrence annually, occur regularly each year, and contribute to maximum monarch population growth.Despite interannual variation in the extent of monarch breeding distribution in Canada, the core of the breeding distribution is in eastern Canada, primarily southern Ontario, and to a less extent southern Quebec and Manitoba (Brower 1995;Wassenaar and Hobson 1998;Flockhart et al. 2013).Southern Ontario is densely populated and all areas in the core of the breeding range are found in land uses that are predominantly agricultural.Given the tendency for monarchs to lay more eggs per host plant in agricultural areas compared with other land uses in the US Midwest (Oberhauser et al. 2001;Pleasants and Oberhauser 2013) and in southern Ontario (Pitman et al. 2018), management strategies that engage with the agricultural sector will almost certainly be required (Thogmartin et al. 2017b;Pitman et al. 2018).For example, one strategy may include compensation for delayed harvest to ensure host plant establishment or monarch recruitment (Drechsler et al. 2007).Alternatively, another option in agricultural landscapes would be to promote habitat restoration in roadsides (Kasten et al. 2016;Pitman et al. 2018).Although the immediate focus of habitat restoration for monarch butterflies in Canada should prioritize productive habitats in southern Ontario, it must be recognized that those potential gains may only contribute marginally to the viability of monarchs in eastern North America (Flockhart et al. 2015).Given the large-scale restoration of host plants necessary to recover monarchs, a diverse portfolio of strategies must be considered across an appropriate geographic extent (Thogmartin et al. 2017b).
There are a few considerations to note when interpreting the results of our work.Our data come from two citizen science programs and their use assumes that sampling occurs randomly throughout the landscape and that detection probability is constant across sites.Citizen science observations are rarely collected from random locations or uniformly across a landscape and, therefore, should be controlled in some manner during analysis (Yackulic et al. 2013).Several approaches have been proposed to account for biased sampling (Phillips et al. 2009).For monarchs, variation in sampling intensity is correlated with the number of potential observers in a given area then controlling for human population density is one manner to account for this bias when mapping breeding distributions (Flockhart et al. 2013).Additionally, the distribution of citizen science observations some years invokes limits to the geographic extent or resolution of the environmental variables that could be considered to explain monarch distributions, especially when considering species distributions at large geographic scales.Detection probability is rarely perfect for any organism but as we are using presence-only observation, any detection probability below unity results in a reduced sample size of occurrence records.However, if detection probability is heterogeneous across the landscape, then it likely correlates with our explanatory variables, which would confound our interpretation of the probability of occurrence with detection probability (Yackulic et al. 2013).Systematic surveys for monarchs in other citizen science projects have focused on monarchs that allow for effort to be controlled, presence-absence data are recorded, and surveys that occur repeatedly over the breeding season (e.g., Monarch Larvae Monitoring Program, Mission Monarch; Ries and Oberhauser 2015) would allow the concurrent analysis of the probability of occurrence accounting for variation in sampling intensity and detection probability (Yackulic et al. 2013).The effects of these limitations in our study, if they exist, would be to lower the probabilities of occurrence for monarchs across Canada and thereby reduce the annual breeding distribution area, but these effects would presumably occur across all years and, therefore, would not have a large impact on variability of the breeding distribution area over time.Future analyses should consider a range of distribution models of higher resolution using measures of monarch density based on counts conducted at precisely georeferenced locations.Our quantitative assessment of the probability of occurrence of monarchs in Canada provides important information to develop successful management strategies because it sets the geographic context for management decisions.The next step is to prioritize management actions while accounting for variation in use of habitat in different land use by monarchs (Pitman et al. 2018), differential recruitment of offspring in these habitats (Oberhauser et al. 2001), and the associated costs of restoring habitat that will maximize population growth of monarchs.Our study supports efforts to identify strategic locations for restoration that will be necessary in regional and national conservation plans for monarch butterflies in Canada.

Fig. 1 .
Fig. 1.The cumulative estimated breeding area (km 2 ) based on the probability of occurrence of monarch butterflies in Canada between 2000 and 2015.The distributions are plotted to represent (a) western and (b) eastern populations based on the Alberta-British Columbia border and for (c) all of Canada.The area estimates are based on the top statistical model of the probability of occurrence testing the geographic, physiological, and habitat hypotheses of distribution based on monarch observations reported in eButterfly (e-butterfly.org/).The probabilities are the lowest values considered in the estimate; for example, a probability of 0.05 represents all areas with a probability >0.05.

Fig. 2 .
Fig. 2. Examples of the predicted probability of occurrence of monarch butterflies in Canada for (a) 2004, a year with below-average breeding distribution; (b) 2010, an average year; and (c) 2012, a year with above-average breeding distribution.In all panels the blue dots represent the occurrence locations in that year used to fit the species distribution model.Map produced in QGIS 2.6 using map data from GADM (gadm.org/).

Fig. 3 .
Fig. 3. Counts of the number of years between 2000 and 2015 where the modeled monarch probability of occurrence was (a) >0.5, (b) >0.25, and (c) >0.05.These values represent the lowest probability value considered; for example, for each pixel, the count is the number of years where the probability of occurrence was 0.25 or greater (in the case of panel b).High counts depict the core of the monarch distribution in Canada that includes southern Ontario and Quebec; monarchs rarely occur in western Canada and only at low probabilities.The locations of all monarch observations (black dots) in Canada used in the study are shown in panel (d).Map produced in QGIS 2.6 using map data from GADM (gadm.org/).

Table 1 .
Attributes of the variables used in the species distribution models.

Table 2 .
Akaike Information Criterion (AIC) values of the top models.All models included the effect of human population density, which added one parameter to each model above those from the explanatory variables listed.See text for details and definitions of the explanatory variables.FACETS Downloaded from www.facetsjournal.comby 70.64.17.116 on 06/28/19 FACETS | 2019 | 4: 238-253 | DOI: 10.1139/facets-2018-0011 243 facetsjournal.com

Table 3 .
Top model parameter estimates (with standard errors in parentheses).See text for details and definitions of the explanatory variables.
FACETS Downloaded from www.facetsjournal.comby70.64.17.116 on 06/28/19River Valley in the Lower mainland and in Douglas-fir/bunchgrass ecosystems in the Central Interior.Monarchs in the west, therefore, likely fly north up the Okanagan River valley in the summer.