Opposing global change drivers counterbalance trends in breeding North American monarch butterflies

Abstract Many insects are in clear decline, with monarch butterflies (Danaus plexippus) drawing particular attention as a flagship species. It is well documented that, among migratory populations, numbers of overwintering monarchs have been falling across several decades, but trends among breeding monarchs are less clear. Here, we compile >135,000 monarch observations between 1993 and 2018 from the North American Butterfly Association's annual butterfly count to examine spatiotemporal patterns and potential drivers of adult monarch relative abundance trends across the entire breeding range in eastern and western North America. While the data revealed declines at some sites, particularly the US Northeast and parts of the Midwest, numbers in other areas, notably the US Southeast and Northwest, were unchanged or increasing, yielding a slightly positive overall trend across the species range. Negative impacts of agricultural glyphosate use appeared to be counterbalanced by positive effects of annual temperature, particularly in the US Midwest. Overall, our results suggest that population growth in summer is compensating for losses during the winter and that changing environmental variables have offsetting effects on mortality and/or reproduction. We suggest that density‐dependent reproductive compensation when lower numbers arrive each spring is currently able to maintain relatively stable breeding monarch numbers. However, we caution against complacency since accelerating climate change may bring growing threats. In addition, increases of summer monarchs in some regions, especially in California and in the south, may reflect replacement of migratory with resident populations. Nonetheless, it is perhaps reassuring that ubiquitous downward trends in summer monarch abundance are not evident.


| 4727
CROSSLEY Et aL. of scales, is often implicated in falling insect numbers (Fox, 2013;Habel, Samways, et al., 2019;Leather, 2018;Sánchez-Bayo & Wyckhuys, 2019). A key local driver has been heavy herbicide and insecticide applications associated with agricultural intensification (Habel, Ulrich, et al., 2019). Urbanization and associated automobile collisions (Baxter-Gilbert et al., 2015;Kantola et al., 2019) and light pollution bring additional challenges (Owens et al., 2020). At global scales, climate change can heighten physiological stress to insects while triggering spatiotemporal misalignment with, or reduced quality of, host plants or other resources (Bale et al., 2002;Jamieson et al., 2012), although even climate change can create variable regions of insect decreases and increases (Crossley et al., 2021;Koltz et al., 2018) such as when increased temperature enables faster population growth. Particularly damaging are cases where local and global drivers both are moving in harmful directions, for example when long-distance migrants must move through increasingly hot and dry regions that also are seeing more intense land use (Saunders et al., 2019).
Monarch butterflies (Danaus plexippus) in North America have become the public face of insect declines (Gustafsson et al., 2015), largely because of the well-publicized diminishing of winter colonies in Mexico and California (Boyle et al., 2019;Pelton et al., 2019).
Monarchs are iconic insects due to their large size, attractive and distinctive coloration, wide range, host association with horticulturally popular milkweeds (Asclepias spp.), and fascinating long-distance seasonal migrations. This has led to the prominent use of monarchs as ambassadors to engage the general public in insect conservation, for example, by facilitating the widespread planting of milkweed in home gardens (Thogmartin, López-hoffman, et al., 2017). However, some of these same traits that make monarchs so charismatic to humans also subject the butterflies to particular risk. Best documented is habitat loss and changing climate at concentrated overwintering sites, which has apparently led to an ongoing, multi-decadal decline of those colonies (Brower et al., 2012;Pelton et al., 2019;Thogmartin, Wiederholt, et al., 2017;Zylstra et al., 2021). A second widely-touted threat is removal of milkweed from agricultural fields within the monarch's core breeding range in the American Midwest, following widespread adoption of glyphosate-tolerant corn and soybean (Stenoien et al., 2018). Thirdly, since migration in the human-dominated world is risky (Wilcove & Wikelski, 2008), their particularly long-distance movements could expose monarchs to multiple threats along the 2 month journey (e.g., deaths from traffic collisions, Kantola et al., 2019;McKenna et al., 2001). Additionally, agricultural and residential pesticides (Olaya-Arenas & Kaplan, 2019) and sensitivity to temperature and precipitation extremes as the climate changes (Lemoine, 2015;Saunders et al., 2018) may be adversely affecting monarchs at various stages of their life cycle.
Altogether, these perceived threats have led to the recent decision by USFWS that federal protection is warranted in the United States (USFWS, 2020). However, evidence is ambiguous whether monarchs continue to be in consistent, recent decline across the annual cycle (i.e., outside of the winter stage), with studies variously reporting steady or falling monarch numbers at different places and seasonal milestones (Brower et al., 2018;Davis & Dyer, 2015;Espeset et al., 2016;Ethier, 2020;Inamine et al., 2016;Ries et al., 2015). Uncertainty about whether breeding populations are continuing to steeply decline, or show some resiliency to overwintering losses in at least some regions or at some stages, complicates efforts to target conservation programs to points in the life-cycle where they will be most effective.
Here, we used the North American Butterfly Association's (NABA) summer citizen-science counts to assess spatiotemporal patterns and drivers of relative abundance of breeding, adult monarchs, and across most of their summer range throughout the United States (east and west) and southern Canada. Prior work with these or similar citizenscience datasets have focused on specific regions of the country, such as the western US , or the Midwest (Zylstra et al., 2021). For a species like the monarch, which has a continental breeding range, it is important to assess the population throughout this large area, so that local or regional hotspots of decline or increase do not bias the interpretation of the entire population's status.
These NABA data are broad in scope, collectively recording 135,705 monarchs at 403 sites across North America, over time periods of 10-26 years from 1993 to 2018. We analyzed NABA data using methods developed for a similar citizen-science program, the Audubon Christmas Bird Count (Meehan et al., 2019), yielding monarch relativeabundance trends that accounted for spatial and temporal variation in sampling effort as well as spatial and temporal autocorrelation among neighboring counts. Our central goals were to (1) quantify trends in monarch relative abundance among NABA sites throughout the United States and southern Canada, and (2) characterize relationships between those trends and two dominant global change factors: agricultural intensification, specifically glyphosate use, and climate change, specifically temperature and precipitation change.

| Butterfly data
We used direct counts of monarch adults from the North American Butterfly Association's summer citizen-science counts (https://www. naba.org/). Butterfly counts are made within a 15-mile (~24 km) diameter circle, typically in July, and are open to participation from the public. For each count event, the abundances of butterfly species are tallied and the sum of associated party hours (a measure of sampling effort that aggregates the number of hours spent by each observer) is recorded. To minimize bias due to differences among sites in the day of year when butterfly counts were conducted, we limited our analysis to butterfly counts that occurred between June 1 and August 31. Prior to estimating trends in abundance, we removed sites that had <5 years of monarch detections and that spanned <10 years (Didham et al., 2020). Lastly, butterfly counts were assigned to 50 × 50 km (2500-km 2 ) cells on a uniform grid covering North America to enable spatial smoothing of estimated relative abundance trends and covariate effects. Grid cells contained an average of 1.21 ± 0.03 and a maximum of three circles. The curated dataset recorded a total of 135,705 monarchs from 403 sites occupying 334 grid cells, over time periods of 10-26 years from 1993 to 2018.

| Modeling relative abundance trends
We modeled monarch counts, y i,k,t , in grid cell i encompassing count circle k during year t as a random variable from a negative binomial distribution. Expected values for counts per grid cell, i,t , were assumed to be a function of spatially structured grid-cell, count-effort, and year effects, plus unstructured variation among count circles.
The linear predictor for i,t took the form Parameters i were modeled as cell-specific random intercepts with an intrinsic conditional autoregressive (iCAR) structure. Parameters i were modeled as spatially structured (iCAR), cell-specific, random slope coefficients for the local effects of effort E i,k,t . Effort was represented by E i,k,t , the number of party hours expended during a count, where a party hour was the count effort of one party of unspecified size for 1 h. Pairing log-transformed expected counts with log-transformed effort in the linear predictor yielded a power function for effort correction, a flexible mathematical form that accommodated a decreasing, linear, or increasing impact of effort on expected counts (Link & Sauer, 1999). Parameters i were modeled as spatially structured, cell-specific, random slope coefficients for a log-linear year effect. Year, represented by T, was transformed before analysis such that max(T) = 0, and each preceding year took an increasingly negative integer value. Given the scaling of effort and year variables, exp( i ) could be interpreted as a cell-specific expected count given one party hour of effort during the final year in the time series. Parameters i,t were modeled as exchangeable random intercepts that accounted for variation in relative abundance per grid cell and year that was not accounted for by the log-linear year effect. The final term in the model, k , was an exchangeable random intercept that accounted for variation in relative abundance among circles, possibly due to differences in habitat conditions or observer experience. This spatially varying coefficient (SVC) model was analyzed within a Bayesian framework using the R-INLA package in R (Lindgren & Rue, 2015;R Core Team, 2021;Rue et al., 2017). For parameters i , i , and i , with iCAR structure (Besag et al., 1991), precision matrices were scaled such that the geometric mean of marginal variances was equal to one, and priors for precision parameters were penalized complexity (PC) priors, with parameter values U PC = 1 and a PC = 0.01 . Precision for the zero-centered, exchangeable, random circle effect, k and grid cell by year effect, i,t , were also assigned a PC prior with parameter values U PC = 1 and a PC = 0.01 . The overdispersion term for the negative binomial count distribution, Φ, was assigned a PC prior with parameter value l = 7. .
Following trend model analysis, posterior medians and symmetric 95% credible intervals were computed per cell for i , i , and i and per cell and year for i,t , by sampling the respective posterior distributions 5000 times. Posterior summaries were then mapped to visualize spatial variation in abundance indices, effort effects, and relative abundance trends.

| Explaining spatiotemporal variation in relative abundance
The North American monarch breeding range spans nearly the entire United States and southern Canada, which includes widely differing landscapes (see Figure S1), including the heavily agricultural region in the Midwest. This region is where 38% of monarchs in Mexico come from (Flockhart et al., 2017), and this is where there has been significant losses of milkweeds due to application of glyphosate in crop fields (Brower et al., 2012). These losses have been proposed as one of the major reasons for the declines in winter colonies in Mexico, because of the temporal synchrony of glyphosate application and colony size decreases (Pleasants & Oberhauser, 2013). As such, determining the impact of glyphosate use on monarch abundance was a priority for us here. In addition, summer climate variables are also known to influence relative abundance of monarchs (Zylstra et al., 2021), and our analyses also incorporated such data.
We used posterior samples along with a subset of the linear predictor to calculate an annual relative abundance index, i,t , per year and grid cell, as i,t = exp i + i + i,t . We then modeled relative abundance indices and their associated uncertainty for grid cell i during year t as a random variable from a gamma distribution.
Expected values for annual abundance indices per grid cell, Ω i,t , were assumed to be a function of spatially structured grid cell, agricultural glyphosate use, average temperature, and cumulative precipitation effects (using data summarized below). The linear predictor for Ω i,t took the form Parameters i were modeled as cell-specific random intercepts with iCAR structure. Parameters i were modeled as spatially structured (iCAR), cell-specific, random slope coefficients for the local effects of glyphosate use, P i,t . Agricultural glyphosate use was calculated as the pounds of active ingredient applied in a county multiplied by the proportion of the county planted in corn or soybean, to account for the expectation that the majority of glyphosate use in a county that is negatively impacting monarch host plants is through applications to corn and soybean acreage (Zylstra et al., 2021). Estimates of pounds glyphosate applied were obtained from the United States Geological Survey-National Water-Quality Assessment Project (USGS, 2022), and corn and soybean acreage were obtained from United States Department of Agriculture National Agricultural Statistics Service (USDA-NASS, 2022) using the "rnassqs" R package (Potter, 2019). To obtain an estimate of agricultural glyphosate use for each grid cell × year, values of glyphosate use offset by the proportion corn or soybean in each county overlapping grid cell i in year t were multiplied by the proportion of overlap with grid cell i. Spatial operations were done in R using functions available from the "rgeos," "raster," and "rgdal" R packages (Bivand & Rundel, 2021;Bivand et al., 2021;Hijmans, 2022).
Maps of glyphosate use (kg active ingredient per acre corn or soybean) in 1993 and 2017 are provided in Figure 1.
Parameters i were modeled as spatially structured (iCAR), cellspecific, random slope coefficients for the local effects of mean annual temperature, Z i,t . Parameters i were modeled as spatially structured, cell-specific, random slope coefficients for the local effects of cumulative precipitation, N i,t . Mean temperature and precipitation data were obtained from CRU TS 4.03 (Harris et al., 2014), which provides monthly gridded estimates at 0.5° latitude/longitude resolution. Mean temperature for grid cell i in year t was calculated as the annual average of monthly mean temperature estimates in year t. Annual cumulative precipitation for grid cell i in year t was calculated as the sum of monthly precipitation estimates in year t.
Maps of mean temperature and cumulative precipitation in 1993 and 2017 are provided in Figure 1.
To propagate uncertainty in relative abundance indices, i,t , during covariate analyses, the analysis was repeated 5000 times using randomly sampled values from the posteriors of α i , τ i , and γ i,t .
Estimates for i , i , and i from each of the 5000 replicates were then used to generate posterior medians and symmetric 95% credible intervals per cell for i , i , and i . Posterior summaries were then mapped to visualize spatial variation in covariate effects.

| RE SULTS
Considering all available NABA data for monarchs across the entire breeding range in eastern and western North America, the median of posterior distributions for relative abundance trends ( i ) pooled across all grid cells suggested an overall annual increase in monarch relative abundance of 1.36% per year. However, there was an 84% chance of the global trend being >0 and a 16% chance of the global trend <0 ( Figure S2). Cell-specific relative abundance trends were generally the most negative in the US Northeast, parts of the Midwest, and in northwest California, and were generally the most positive in the US Southeast and Northwest (Figure 2a). Only 11 of the 334 grid cells exhibited relative abundance trends whose 95% credible intervals did not overlap zero, 10 of which were positive trends in The lack of strong trends in the core breeding range is in contrast to studies that focus on winter colony size as measures of population abundance, where there are clearly multi-decadal declines (Brower et al., 2012;Oberhauser et al., 2017;Semmens et al., 2016;Thogmartin, Wiederholt, et al., 2017;Zylstra et al., 2021), but is in general agreement with various breeding season and fall migration studies that have shown high variability in monarch abundance trends (Table S1) Even though the most recent evidence indicates the monarchs west of the Rocky Mountains should not be considered a separate population (Freedman et al., 2021;Talla et al., 2020), the assessment of monarch abundance in the west has traditionally been via counts of wintering monarchs along the California coast (Espeset et al., 2016;Pelton et al., 2019). However, as we found with the larger cohort of monarchs east of the Rockies, the trend of diminishing wintering colonies in California does not appear to mirror long-term trends in breeding monarch abundance to the north or northeast, either in Oregon or Idaho (Figure 2a). In fact, the (admittedly) few locations with long-term data in that region indicate an overall increasing trend (sampling the posterior distributions of i for grid cells overlapping Oregon and Washington revealed an 86.7% probability that the monarch relative abundance trend was >0). However, we do note that the NABA data we had access to ended in 2018, before a dramatic drop in colony size in the winter of 2020/2021 James, 2021), so we do not know how this may (or may not) have affected breeding monarchs in the northwest. We also note that the NABA data showed a small region in central California where summer numbers are declining, Cyan and pink highlighting denotes estimates whose 95% credible intervals were greater or less than zero, respectively. Maps of (d) glyphosate effect estimates ( i ), (e) cumulative precipitation effect estimates ( i ), and (f) mean temperature effect estimates ( i ) among grid cells. For (d-e), cyan and pink highlighting denotes estimates whose 95% credible intervals were greater or less than zero, respectively.
which is consistent with other long-term surveys from that same region (Espeset et al., 2016), and this consistency provides confidence in the NABA data. Further, reasons for declines in wintering monarchs in California have been the subject of ongoing debate, with some speculating that western monarchs may be transitioning to a less migratory lifestyle in California, which is being fueled by homeowner plantings of non-native milkweed that thrives yearround (Davis, 2022;James et al., 2022). Regardless of the reason, the discrepancy between wintering numbers and breeding abundance in the west, like that of the east, argues that overall population assessments should be based on multiple sources, and from across life stages.
Our analysis indicated that glyphosate use, while an important contributor to local monarch numbers, is significantly affecting only a portion of the summer breeding range (portions of upper Midwest, Figure 2d). The initial rapid increase in glyphosate use in midwestern corn and soybean, which likely devastated weedy milkweeds in those fields, has now leveled (Zylstra et al., 2021), such that harmful indirect effects of herbicides on monarchs may no longer be increasing in magnitude. This suggests that the loss of agricultural milkweed in the US Midwest will not inevitably lead to ongoing drops in summertime abundances (Agrawal & Inamine, 2018). In fact, a recent inventory of the western half of the United States revealed billions of previously uncounted native milkweeds that are available for monarchs (Spaeth et al., 2022), supporting the notion that there are sufficient hostplants to maintain a stable summer population throughout much of the breeding range.
Recent analyses indicate that changing climate is driving increases and decreases in overall butterfly numbers across North America (Crossley et al., 2021;Forister et al., 2021), and, there is evidence that temperature and precipitation in North America is indirectly and positively impacting abundances of overwintering monarchs, via positive effects on breeding monarch population size (Zylstra et al., 2021).
Such counterbalancing effects of seasonal temperature and precipitation appear to be common in butterflies (Davies, 2019;Konvicka et al., 2021;Roland & Matter, 2016). In line with this, we found a pattern of increasing monarch relative abundance with increasing average temperature in the northern US, with the strongest effects evident in the midwestern US (Figure 2f), where glyphosate use appeared to have the strongest negative effect (Figure 2d). Positive and negative effects of precipitation were also evident, but this signal was statistically less robust (Figure 2e). The eastern US and Canada, the area corresponding to the major monarch summer breeding ground for the eastern subpopulation, have generally seen increases in precipitation and only modest increases in summer temperature (IPCC, 2018), conditions that have apparently been providing favorable conditions for many butterfly species (Crossley et al., 2021).
However, Texas and the northern portions of Mexico, a vital corridor region, have seen recent pronounced increases in temperatures (Cuervo-Robayo et al., 2020) which could be negatively affecting survivorship during the arduous southward migration.
Overall, our findings suggest monarch populations may have some ability to recover, on average, from declines at overwintering colonies. Of course, the total loss of overwintering monarchs would make it impossible for any summer rebound to be ignited, and there almost certainly is some inflection point well before total winter extinction where spring migrants would be too few to reliably spark a summer resurgence. This would leave only the year-round resident monarch populations, with the loss of the epic migrations that inspire much human interest in monarchs as conservation icons. For those monarchs that do return northward in the spring, our results argue that following the winter period, monarchs experience high population growth, perhaps facilitated by reduced intraspecific competition among larvae. Indeed, monarch larvae are known to exhibit negative interactions with conspecifics, including egg cannibalization (Brower, 1961), aggression (Collie et al., 2020), and oviposition avoidance on optimal host plants (Jones & Agrawal, 2019), behaviors that are presumably reduced under the smaller population returning from recent years in winter migration. Considering the general lack of widespread breeding season declines found here, our evidence suggests, alongside the ongoing declines at winter colonies, that monarchs must be experiencing increasingly higher levels of mortality during their fall migration. Contrasting evidence of no change in the number of tagged monarchs returning to Mexico in the fall suggested otherwise (Taylor et al., 2020), but that finding remains contested due to difficulties in accounting for changing tagging effort through time (Fordyce et al., 2020). In support of our assessment, a recent study of the monarch parasite, Ophryocystis elektroscirrha, has shown that nation-wide prevalence has increased in the last 15 years, and that this increase is leading to considerable migratory losses and corresponding reductions in winter colony sizes (Majewska et al., 2021). Therefore, conservation attention along the migration routes, and/or actions that reduce parasite transmission, may be more imperative for the monarch's long-term survival compared to efforts directed at the breeding grounds.
Our data were collected by citizen-scientists, a method that requires careful use (Burgess et al., 2017), but that nonetheless enables inquiry at spatiotemporal scales otherwise unachievable by individual research groups (e.g., Herremans et al., 2021). The number of party hours spent monitoring butterflies in the North American Butterfly Association dataset increased on average by 1.2% (±0.3%) per year between 1993 and 2017 ( Figure S3). However, our analyses accounted for annual variation in sampling effort while allowing for a variety of relationships between increasing sampling effort and monarch counts, following methods developed to analyze conceptually similar Audubon Christmas Bird Counts (Meehan et al., 2019). Importantly, we did not find evidence of increasing or decreasing trends in sampling effort around sites dominated by cropland or forest, suggesting that changes in sampling effort have neither masked declines nor exaggerated abundance increases ( Figure S4). Furthermore, we found that the local effects of sampling effort exhibited an increasing impact on Beyond monarchs, the conservation of insects has received far less attention than most other taxa, despite the ubiquity of insects in terrestrial ecosystems. Undoubtedly, citizen-science efforts targeting the charismatic monarch have exposed many non-scientists in North America to the importance of insects and the value of their conservation. Given our results, we suggest that there could be considerable ecological gain from broadening citizen scientists' attention to also consider the many butterfly species that do appear to be experiencing major summer declines across North America.
For example, the summer butterfly count data suggest that Lycaeides melissa is declining across much of its broad range (Figure 3), and even the well-known west coast painted lady, Vanessa annabella, appears to be faring worse than the monarch (Figure 3). In fact, of the 456 butterfly species tracked by NABA, there are 320 species with trends less positive than monarch butterflies (Crossley et al., 2021).
More broadly, our results are consistent with other recent analyses of large-scale insect data that have also revealed complex and heterogeneous spatiotemporal patterns of insect decline. For example, a warming climate in Europe is shifting some moth ranges northward, with species unable to do so declining, but leading to a net range increase overall (Fox et al., 2021). Similarly, recent drops in UK moths seem modest relative to increases seen over prior decades (Macgregor et al., 2019), leading to no net change over time. In North America, close examination of long-term insect counts revealed declines in some taxa, but increases in others (Crossley et al., 2020).
The same is true with butterflies, where species declines in western North America may be at least partially offset by abundance increases elsewhere on the continent (Crossley et al., 2021), again, leading to no net change despite troubling declines in some locations and/or for some taxa. Our analyses show that for monarchs, for now, summer abundance increases appear sufficient to buffer winter declines. It will be increasingly important to understand complex interactions among species traits and mechanistic drivers, in order to understand and successfully predict how an ever-more-rapidly changing environment will impact the future persistence of monarchs and other insects.

ACK N OWLED G M ENTS
We thank Anurag Agrawal for helpful feedback on an early version of the manuscript. We are grateful to the thousands of volunteers who have participated in butterfly surveys for the North American Butterfly Association over the years, which has made research like this possible. This paper was developed following many lengthy and stimulating discussions with colleagues who study monarch butterflies, and who we believe all have the same goal: to understand and identify potential threats to the monarch population in North America.

FU N D I N G I N FO R M ATI O N
We acknowledge funding from USDA- NIFA-OREI 2015-51,300-24,155, USDA-NIFA-SCRI 2015, and USDA-Hatch DEL00774.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data (relative abundances and environmental covariates) and R code supporting the findings of this study are available on Dryad (DOI F I G U R E 3 Monarch abundance trend compared with other common North American butterflies. Histogram depicts median abundance trends (%/year) of >450 species monitored by the North American Butterfly Association. Trend for Danaus plexippus (+0.7%/year) is highlighted and compared with three other well-known species, Lycaeides melissa (−2.0%/year), Vanessa annabella (−7.8%/year), and Vanessa atalanta (−1.1%/year). Trends based on sites where butterflies were recorded at least five times over a span of 10 years. See Crossley et al. (2021) for details on trend estimation. All butterfly species trends are available in Table S6. https://doi.org/10.5061/dryad.7pvmc vdwb). Raw data are available from the North American Butterfly Association (NABA) at https:// www.naba.org/