Data

We used two different time series depicting dynamics of the Eastern, migratory population of monarch butterflies. First we used the log of the total extent (ha) of overwintering forest area occupied in Mexico per year from 1993–2014 measured by the World Wildlife Federation- Mexico and the Monarch Butterfly Biosphere Reserve (MBBR)5. Second, we used the log of the estimated total amount of egg production in the Midwest per year from 1999–20126 and extended this through 2014 for the current analysis (Extended Data Table 1). The egg production per year was based on the average estimated eggs per milkweed stem for that year multiplied by the number of available milkweed stems on the landscape in that year. Eggs per stem estimates come from weekly monitoring data reported to the Monarch Larva Monitoring Project (MLMP - http://www.mlmp.org) from citizen scientists throughout the monarch range. We used the eggs per stem value from the week of peak egg production as an index for the production estimation. The number of available milkweed stems was based on estimates of the density of milkweeds in different habitats based on surveys23,24 and the area on the landscape occupied by those habitats using U.S. Department of Agriculture (USDA) databases25 and a function describing the decline of milkweeds in agricultural fields6.

Estimating Population Parameters

We developed a multivariate first-order auto-regressive state-space model10,11,12 to generate monarch population parameter estimates for quasi-extinction risk forecasting. We chose a Bayesian modeling approach because the resulting parameter posterior distributions provide a complete characterization of parameter uncertainty that can be seamlessly propagated through to our quasi-extinction risk analysis. In so doing, we can account for uncertainty in quasi-extinction risk due to uncertainty in the parameters controlling population change. Moreover, the hierarchical nature of state-space modeling is easily handled through Bayesian estimation.

The log-scale population modeling framework takes the following form:

In the above set of equations, x represents the state process (estimated log of the true size of the overwintering population in Mexico) across all years t for which we have data. The state process evolves from one year to the next according to a mean population growth rate, and associated random yearly deviates to growth, w t , which we assumed to be normally distributed with a mean of 0 and a standard deviation q (process noise). The exception to this is that we must directly estimate the population state in the first year (1993) using an uninformative uniform prior, since there is no prior year to evolve from. Note that , the average annual (non-logged) population growth rate, where λ values of <1 result in population decline, while values of >1 results in population growth.

Values of m t are the log of yearly estimates of Mexican overwintering habitat occupied, which we assumed to deviate from the state x t by v t . 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. We used this parameterization based on the assumption that process noise in the time series is greater than the measurement error associated with the Mexican overwintering data. This assumption is based on a consensus among the authors that measurement error is exceeded by variation in population sizes caused by demographic and environmental variation and because process noise is typically the predominant form of variability in time series’ of wild populations26. As such, measurement error in m is defined to be a proportion (p) of process noise q.

Similarly, the log of annual estimates of Midwestern egg production, e t , are assumed to deviate from the state x t by a, a scaling parameter11 that shifts the egg production index to the same scale as the overwintering habitat 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). Note that, in the context of the linear modeling framework outlined above, the parameter a is essentially an intercept term.

We fit our model using R27 and JAGS28 (Just Another Gibbs Sampler) and assessed convergence by examining parameter trace-plots and calculating Gelman-Rubin diagnostics using the CODA29 package in R. Parameter estimates are provided in Extended Data Tables 1 and 2.

For the purposes of quasi-extinction risk forecasting, we report the following posterior estimates after removing variation caused by measurement error:

x 2014 , the estimate of the true size of the Mexican overwintering population during the winter of 2014–15,

, the monarch population growth rate, and

q, the process noise associated with the monarch population.

The process noise, q, represents year-to-year variability in the population after accounting for overall trend through time (ū, the mean growth rate in log space) and after removing variation caused by measurement error. Note that using these model estimates, we can simulate the population forward in time from its estimated current size (x 2014 ).

Calculating Quasi-extinction Risk

Because each posterior draw from the Bayesian state-space model represents a complete set of likely parameter values, we can use all the posterior draw sets to generate annual probabilistic quasi-extinction risk estimates that account for both uncertainty in population parameters and uncertainty due to the stochastic population process (random yearly growth or decline due to process noise). For each posterior draw i, we simulate the population 20 years into the future 1000 times, starting at x 2014i , using the growth rate and process noise q i . For each of i simulation sets, we subsequently calculated the percent of runs that fall below a given quasi-extinction threshold (described below). Because we carried out this exercise for each i posterior draw, we can subsequently generate median and 95% credible interval estimates of quasi-extinction risk that account for population parameter uncertainty. Using a model run with a burn in of 4e5 iterations and a sample window 10e5 iterations (3 separate chains) thinned by 600, our model achieved satisfactory convergence based on both visual inspection of trace plots and Gelman and Rubin diagnostics30; the Potential Scale Reduction Factors for all parameters were below 1.05.

The quasi-extinction threshold, or population size at which extinction of the Eastern monarch migration becomes inevitable, is unknown. We do know that monarch ecology exhibits at least two characteristics that suggest the likelihood of a strong Allee26 effect and therefore the existence of an extinction threshold: tightly clustered overwintering colonies convey important microclimate advantages that diminish as colony size decreases19,31 and the increased efficiency of locating mates in overwinter aggregations32,33,34. Diminishing colony size can therefore result in higher winter mortality rates and lower fecundity in the spring, which can cause the population growth rate to drop below replacement.

Regarding the size of winter aggregations specifically, expert opinion among the monarch biologist authors of this manuscript (Taylor, Oberhauser, Pleasants) favored an extinction threshold of no less than 0.05 ha, with most favoring a threshold closer to 0.25 ha. For reference, the smallest observed size of a Mexican overwintering colony is 0.01 ha5 and the minimum number of colonies that has been observed in Mexico in any given year is 75. To be sure that we have captured the real quasi-extinction threshold, we opted to consider a range of values from 0.01 ha (1 viable colony, equivalent to approximately 1 occupied tree) to 0.25 ha. Our lowest (most optimistic) quasi-extinction threshold is therefore equivalent to the smallest observed colony size on record, which permits the loss of all but one core colony and the associated redundancy of multiple colonies before quasi-extinction occurs. The presence of multiple colonies provides a buffer against extinction during winter storms because of their variable storm severity within the overwintering region14. Our least optimistic estimate of the threshold suggests that quasi-extinction could occur well before the population declines to a single core colony at the minimum observed size. These values are intended to bookend the plausible range of the extinction threshold based on the best available information.

A mechanistic approach to estimating suitable threshold values was presented by Wells et al.33, who developed a model of monarch butterfly mating in California overwinter clusters that demonstrated a relationship between mating success and overwintering density. They found that reproductive success was highest for aggregations over 250,000 individuals. As aggregation size dropped below 250,000, reproductive success started to decline and the rate of decline increased substantially below 50,000. They further noted that stable overwintering aggregations in California normally fall within this range (50,000–250,000), but acknowledge that they did not account for predation, which is a significant factor at Mexican overwintering sites and would likely require shifting this range higher. In addition, this model does not account for the population-level benefits of having multiple colonies.

The number of individuals present in overwintering colony areas is strongly dependent upon the density of monarchs per hectare. There are five published estimates of monarch overwintering densities, ranging from 6.9–60.9 million monarchs per hectare13,14. At the low end of this range, a quasi-extinction threshold of 0.05 ha would yield ~345,000 individuals, which may be just small enough to impact mating success according to the Wells et al.33 model, once higher Mexican predation rates are accounted for. A threshold of 0.01 ha yields just 69,000 individuals at the lowest density and ~324,000 at the average density, which is more firmly in the realm of reduced mating success according to the Wells et al.33 model.

Estimating a target population size

We ran the quasi-extinction risk simulation across a range of initial starting population sizes and thresholds to develop associated quasi-extinction risks for each combination of values. The quasi-extinction risk simulation was run at initial population sizes ranging from 1 ha–10 ha, in 1 ha increments. We ran these simulations for quasi-extinction thresholds of 0.01, 0.05, 0.15 and 0.25 ha over 10 and 20 years. Because the intent of this exercise was to inform the selection of a recovery goal based on population sizes that confer protection against quasi-extinction risk, we conducted our simulation exercise under the assumption that population declines have been halted and the annual growth rate (λ) of the population is 1. All quasi-extinction events in our simulations are thus exclusively a function of process noise.