Introduction

The migratory population of monarch butterflies (Danaus plexippus) in eastern North America underwent a long period of decline from the mid-1970s to early 2000s (Zylstra and others, 2020) as habitat in the core of their breeding range deteriorated (Pleasants and Oberhauser, 2013; Pleasants and others, 2015; Zaya and others, 2017), leading to increased concerns regarding species imperilment (Semmens and others, 2016). Deterioration of breeding habitat, mainly through removal of the monarch’s larval host plant, milkweed (Asclepias spp.), is associated with the application of glyphosate to corn and soy (Pleasants and Oberhauser, 2013), which have been genetically modified to tolerate herbicide application (Stenoien and others, 2018). Adoption of herbicidally resistant corn and soy now exceeds 90 percent, turning what was once important habitat for monarchs into a nonhabitat matrix through which they must find comparatively sparser remnant milkweed.

Given that milkweed from corn and soy fields is nearly eradicated, continuing decline in population attributable to this threat should have abated about 2006–08 (Thogmartin and others, 2017c; Saunders and others, 2018; Pleasants and others, 2023). It seems the carrying capacity for monarch butterflies is now substantially lower (Pleasants, 2017) with monarch butterflies now varying stochastically between this new upper bound and some lower absorbing boundary. These upper and lower bounds act as basins of attraction (that is, alternative states or conditions leading to long-term behavior in population dynamics). The carrying capacity, or upper bound, has been estimated at 3.2 hectares (ha) of overwinter habitat occupied in central Mexico (Pleasants, 2017); this area estimate serves as a proxy to the species abundance (Thogmartin and others, 2017a). The lower bound is generally labeled an extinction basin and may exist at some level of population size above zero but below a level from which the population cannot resume growth (Semmens and others, 2016). We do not know what this lower-level abundance is, but it is almost certainly smaller than the smallest population size observed to date (2022), 0.67 ha (Rendón-Salinas and Tavera-Olonso, 2014). A population below this threshold level of abundance is identified as quasi-extinct (Semmens and others, 2016), and continental-scale migration of monarch butterflies is not expected to continue at this population size.

Semmens and others (2016) previously assessed extinction risk for the eastern migratory population of monarch butterflies, determining that risk at the 10-year mark was 42 percent (7–88 percent at the 95-percent credible interval [CI]) and at the 20-year mark was 62 percent (9–98 percent at the 95-percent CI). This study was completed before it was suspected monarch butterflies may have entered into a relatively stationary dynamic. Given this apparent abatement in the decline (Thogmartin and others, 2017c), re-assessing the risk faced by the population at its lower but perhaps stationary abundance is appropriate. At this stationary abundance, we can hypothesize that the current risk for the population derives principally from interannual variation in abundance relative to the population’s lowered abundance. The current risk faced by monarch butterflies during this current period of abundance (after 2006 to 2022) was calculated. This time series of abundance was examined for signs that the population may move from the current basin of attraction, as defined by the carrying capacity, to the extinction basin from which the population presumably cannot recover.

The use of generic early warning signals for detecting the proximity of a system moving from one stable state to another is increasingly common (Scheffer and others, 2009; Dakos and others, 2012; Dai and others, 2013). These studies indicate advanced warning can be revealed based on fluctuations in time series, such as temporal variation and autocorrelation. The predictive power of these leading indicators is, of course, limited by the need for long-term observations, which are lacking here. Nevertheless, skewness, kurtosis, autocorrelation, and variance were calculated to determine from these limited data any insight into possible slippage into extinction. If the population were “flickering” ahead of extinction, one would expect variance, autocorrelation, and skewness to increase over time.

Methods

A multivariate first-order auto-regressive state-space model (Holmes and Fagan, 2002; Holmes and others, 2012) described by Semmens and others (2016) was used to replicate their analytical approach. This Bayesian modeling approach permitted separation of observation and process error in these data to allow for generation of probabilistic extinction risk. Explicitly estimating observation and process error is needed to avoid unnecessarily inflating risk deriving from measurement error. Measurement error at the overwintering sites is non-negligible, originating primarily from difficulty in measuring the variable density of monarchs within each colony; published density estimates range greater than about an order of magnitude, from 6.9 to 60.9 million monarchs per hectare (Calvert, 2004; Brower and others, 2004, Thogmartin and others, 2017a). Process error can markedly affect extinction risk because environmental variability can cause populations to stochastically hit extinction thresholds well before a deterministic decline would indicate; this risk from environmentally induced variance, or stochasticity, in the species dynamics is particularly true for small populations (Dennis and others, 1991). Parsing observation error from process error avoids inappropriate amplification of risk from measurement error and allows estimation of risk solely from that relating to environmental stochasticity.

Population data used to inform the model included the annual area (in hectares) of forest occupied by overwintering colonies in Mexico (2006–22; Rendón-Salinas and others, 2023) and the annual number of butterflies counted in summer by the Illinois Butterfly Monitoring Network (2006–21; Illinois Butterfly Monitoring Network, 2022). Semmens and others (2016) used counts of eggs and larvae as reported by the Monarch Larva Monitoring Program to parse the observation process from the state process. The counts of adult monarchs in Illinois were used because they covered a greater range of years than the data available from the Monarch Larva Monitoring Program and were more highly correlated to the overwinter estimate of population size (correlation coefficient [r]=0.77, probability [p]<0.0001; Pleasants and others, 2023).

As in Semmens and others (2016), the log-scale population model takes the following standard form:

Values of v_t follow a normal distribution with a mean of 0 and a standard deviation of q×p (measurement error), where q is the process noise and p is a proportion parameter. As in Semmens and others (2016), this parameterization was used based on the assumption that process noise in the time series is greater than the measurement error associated with the Mexican overwintering data. Process noise is typically the predominant form of variability in time series of insect populations (Fox and others, 2019).

Similarly, the log of annual estimates of Illinois adult monarch counts, c_t, is assumed to deviate from the state, x_t, by a, a scaling parameter shifting summer counts to the same scale as the overwintering habitat area index, and f_t, where values of f are assumed to be normally distributed with a mean of 0 and a standard deviation of r (measurement error):

Using code amended from Semmens and others (2016), this model was fit using R (version 4.2.2, R Core Team, 2023) and Just Another Gibbs Sampler (known as “JAGS;” Plummer, 2003) and checked convergence using the Brooks–Gelman–Rubin potential scale reduction factor where the coefficient of determination (R) is less than or equal to 1.05 was considered indicative of model convergence (Brooks and Gelman, 1998). To forestall concern that the chosen time series was biased by a high initial abundance and to assess the sensitivity of the extinction results, all analyses were repeated covering 2 fewer years (2008–22, n=15 years); results were indistinguishable in almost all cases and are, therefore, not reported. For example, extinction risk at year 10 was 0.287 (0.019–0.810 at the 95-percent CI) for the 17-year set versus 0.274 (0.014–0.853) for the 15-year set, with the difference a result of a slightly lower rate of decline (estimated annual rate of growth [λ]=0.95 [0.64–1.42]) for the shorter set. Stationarity in the time series was tested for with Kwiatkowski–Phillips–Schmidt–Shin (KPSS) and Dickey-Fuller tests.

Because of the limited data, metric-based early warning signals relating to autocorrelation, standard deviation, skewness, and kurtosis were calculated (Dakos and others, 2012). To examine the robustness of these warnings, varying degrees of smoothing were examined (3, 5, and 7 years). Data and analysis code are available at Thogmartin (2023).

Results

Mean population size of the eastern migratory population of monarch butterflies was 3.16 ha between 2006 and 2022 (standard deviation [SD]=1.73), declining 67.8 percent from 6.87 to 2.21 ha (fig. 1; 52.1-percent decline in 2007–22; 56.3-percent decline in 2008–22).

The estimated annual rate of growth (λ) for 2006–22 was 0.93 (0.67–1.30; hereafter, median and 95-percent CI; fig. 2). About two-thirds (68.3 percent) of the posterior distribution fell below λ=1 (fig. 2), indicating about 2:1 odds in favor of there being a declining population. Tests for stationarity indicated a similar degree of uncertainty in trend with the KPSS test indicating stationarity (H₀=time series is stationary; KPSS trend=0.1202, p=0.098) and an augmented Dickey-Fuller test indicating non-stationarity (H₀=time series is non-stationary; Dickey-Fuller=−2.12, p=0.526).

Autocorrelation in the time series was relatively low (fig. 3, left panel), indicating a rather moderate return toward equilibrium with each move away. However, relating patterns in annual change to population size indicated that most decreases in abundance from one year to the next occurred when the population size was less than 3 ha (9 of 10 decreases), whereas most increases (4 of 6) occurred when the population size was greater than or equal to 4 ha.

Process noise (in units of standard deviation) was estimated at 0.43 SD (0.22–1.02 CI) and overwintering habitat area observation error at 0.18 SD (0.01–0.47 CI). The estimated probability of extinction over 10- and 20-year periods based on an extinction threshold of 0.25 ha was 28.7 percent (1.9–81.0 CI) and 52.0 percent (3.2–97.7 CI), respectively (fig. 4); note, however, the wide uncertainty, reflecting uncertainty in the estimate of trend.

In contradiction to the hypothesis of an impending state change, autocorrelation declined over the period (Kendall τ_{7 year}=−0.6), whereas skewness (Kendall τ_{7 year}=0.644) and kurtosis (Kendall τ_{7 year}=0.778) increased, in comportment with the hypothesis (fig. 5). Warning signals at finer degrees of smoothing (3 and 5 years) were similar but generally of a lower magnitude than the coarsest smoothing. Variance as measured by standard deviation provided contradictory evidence across the levels of smoothing.

Discussion

Populations can proceed to extinction in many ways. For small populations, irregular and unpredictable perturbations become increasingly important in determining a population’s risk to extinction (Shaffer, 1987; Lande, 1993; Fagan and Holmes, 2006). At low abundance, risk from annual fluctuations in abundance may accumulate, even if the trend in underlying abundance is nonnegative. So too for the eastern migratory population of monarch butterflies. At its reduced abundance, risk for this species accumulates year over year so that by year 20 (about 2042) there is a predicted 50+ percent probability of dropping below a putative extinction threshold of 0.25 ha. Much of this risk comes from the stochastic nature of its population dynamics. Other early warning signals of impending extinction provide a mixed picture; autocorrelation in abundance seems to be declining with time, whereas skewness and kurtosis (measures of extremes) are increasing. Together these metrics indicate a system increasingly dominated by stochastic variation, heightening risk (Fagan and Holmes, 2006; Melbourne and Hastings, 2008; refer to Appendix 1).

Mean population size for 2006–22 was 3.16 ha, in accordance with the carrying capacity estimated by Pleasants (2017). The λ_2006–2022 of 0.93 (0.67–1.30) is indistinguishable from that reported by Semmens and others (2016) for the 1993–2014 period (λ_1993–2014=0.94, 0.69–1.30 CI). Process error was similar among periods, q–σ_2006–2022=0.43 versus q–σ_1993–2014=0.49; measurement error, however, was nearly half that (m–σ_2006–2022=0.18, from surveys of adult butterflies, versus m–σ_1993–2014=0.44, from egg production data). Median estimates of extinction risk at 10- and 20-years for a quasi-extinction threshold of 0.25 ha were slightly lower for 2006–22 (0.29 and 0.52, respectively) compared to those reported by Semmens and others (2016) for 1993–2014 (0.42 and 0.62, respectively); uncertainty in risk, however, was large in both periods, precluding differentiation.

Perhaps the most concerning observation is that small populations (less than 3 ha) more frequently continued their decrease whereas most increases in abundance only came about when the monarch butterfly population was larger (greater than or equal to 4 ha). Decreasing abundance when the population is particularly small may signal the pull of the extinction vortex (Gilpin and Soulé, 1986) via Allee effects (Courchamp and others, 1999), perhaps by lowering mating success (Kuussaari and others, 1998) or delaying migration to breeding grounds (Edwards and others, 2021).

These analyses are, however, based on a time series of only 17 years. A series of years longer than 17 would have included years in which the monarch butterfly population was still declining from milkweed eradication in row crops. Shorter series would have made risk assessment increasingly difficult (Foley, 1994; Lotts and others, 2004; Fox and others, 2019). Based on these limited data, though, it is reasonable to conclude the eastern migratory population of monarch butterflies remains at substantial near-term risk of dropping below quasi-extinction thresholds from which the possibility of recovery is uncertain.

Risk to this comparatively small population would be mitigated by either an increase in abundance or diminishment of the stochastic variation experienced by the species. Considerable evidence associates the annual variation in monarch butterfly abundance to weather (Thogmartin and others, 2017c; Saunders and others, 2018), a stressor that is likely irreducible given current management tools. Thus, the most promising means of mitigating risk is increasing population size (Semmens and others, 2016), principally through restoration of milkweed to regions where it has been eradicated by herbicidal application to corn and soy agriculture (Pleasants, 2017; Thogmartin and others, 2017b; Solis-Sosa and others, 2021). Strategies are in place for exactly this purpose (Midwest Association of Fish and Wildlife Agencies, 2023), largely focused on milkweed restoration on private agriculture and conservation lands (Janke and others, 2021; Thakur and Hurley, 2023), protected natural lands (Kurzejeski and others, 2020; Lukens and others, 2020), transportation and utility rights of way (Kasten and others, 2016), and urban and developed lands (Johnston and others, 2019). Restoring milkweed provides monarch butterflies the breeding substrate it requires to loosen the grip of stochastic variation as it grows away from the absorbing boundary of extinction.