The mechanics of DDRP are detailed in Barker et al. (55). Briefly, DDRP uses a process-based modeling approach in which degree-days and climate stress are calculated daily and accumulate over time to model phenology and climatic suitability, respectively ( Figure 1 ). Required inputs include gridded daily minimum and maximum temperature data (T_min and T_max , respectively) for a time frame of interest and a species parameter file (summarized in Table 1 ). Models include four separate life stages (egg, larva, pupa, and adult) plus a separately parameterized overwintering stage.

Phenology models often incorporate variation in phenology using a user-defined number of time-distributed cohorts that reflect the distribution of times in which the overwintering stage completes development. Thus, the daily accumulation of degree-days and accompanying life cycle transitions are carried out in separate bins starting from the completion of development of the overwintering stage of each cohort. For computational simplicity, individuals within each cohort are assumed to develop in synchrony. We typically apply seven cohorts because they approximate a normal distribution of development completion times and allow for reasonable model run times since the phenology model is run separately for each cohort.

DDRP uses accumulated stress unit totals for the year to delimit two levels of exclusion (i.e., moderate vs. severe stress exclusion) at any given location ( Figure 1 ) (55). Areas not excluded by either moderate or severe temperature stress are part of the potential distribution and therefore at highest risk of establishment. Moderate stress may inhibit long-term establishment, with a high probability that short-term establishment could occur during favorable years, whereas areas under severe stress would likely prevent even short-term (one complete year) establishment.

Over the daily time-step, DDRP produces outputs that include maps of the potential distribution, number of generations, life stages present, and dates of phenological events. Additionally, it integrates maps produced by the phenology and climatic suitability model to provide insight into phenology (i.e., life stage present or the date of specific event) only in areas at risk of establishment.

Parameter values used in the final phenology model are provided in Table 1 .

We developed and validated a climatic suitability model for EAB using a combination of eco-physiological information and georeferenced presence records from China (N = 148), European Russia and Ukraine (hereafter Europe; N = 139), and North America (N = 2947). Presence records were derived from peer-reviewed literature, theses, reports, the NAPIS database (https://cmr.earthdata.nasa.gov/search/concepts/C1214608223-SCIOPS, accessed 29 Aug 2022), the Global Biodiversity Information Facility (91) (http://gbif.org, accessed 27 Dec 2022), and CERIS Pest Tracker (https://pest.ceris.purdue.edu, accessed 27 Dec 2022). Of the 3234 presence records, 1068 were spatially resolved only to the level of a county or city, which is typically a coarser scale than climate data. However, predictive uncertainty was not an issue for these locations because none occurred near the predicted range limits for EAB (see ‘3.2 Results’). Presence records were deposited in Zenodo (https://doi.org/10.5281/zenodo.7493142).

Climatic suitability models used daily estimates of T_min and T_max for the most recent 20-year time frame available for each focal region. These included (1) the CDAT dataset for China for 1999 to 2018 at a spatial resolution of 0.1° (92)(accessed on 28 Jun 2022; https://zenodo.org/record/5513811#.YtnLNXbMIuU), (2) the Daymet dataset for North America for 2002 to 2021 at a spatial resolution of 1 km² (93, 94) (https://daymet.ornl.gov/getdata, accessed 13 Jul 2022), and (3) the E-OBS dataset for Europe for 1999 to 2018. E-OBS data for 2019 and 2020 were excluded because of a large amount of missing data in European Russia. Daymet data were chosen over PRISM data because they include climate estimates for Canada; however, they were cropped to a maximum latitude of 60°N and aggregated to a spatial resolution of 4 km² to decrease model run times.

We produced phenological event maps of first adult emergence and first egg hatch for EAB in 2021 for North America and Europe to provide insight into its phenology and potential distribution during a particularly warm year. The summer of 2021 was tied with 1936 for the warmest on record in CONUS (NOAA website https://www.ncei.noaa.gov/news/national-climate-202108, last accessed 22 Jan 2023) and was the hottest on record for Europe (European State of the Climate 2021, https://climate.copernicus.eu/esotc/2021, last accessed 22 Jan 2023). The model used Daymet data for North America and E-OBS data for Europe.

We tested for trends in the DOY of first adult emergence and climate stress exclusions over a 20-year period to explore potential changes in EAB’s phenology and potential distribution due to climate warming. Tests used the ‘Mann-Kendall’ function in the Kendall R package v. 2.2.1 (103) and excluded grid cells with missing data (Europe only). The Mann-Kendall statistical test (104, 105) is a rank-based non-parametric method that is widely used for time series of remote sensing data (106, 107). Analyses of climate stress exclusions were only conducted for North America because Europe was climatically suitable for most years (see ‘3.2 Results’). We used absolute values for climate stress exclusion values in DDRP (0 = no stress exclusion, −1 = moderate stress exclusion, and −2 = severe stress exclusion) to make trends easier to interpret (e.g., a decreasing trend would indicate decreasing climate stress over time). Trends were considered statistically significant at the P = 0.1 level because Mann-Kendall tests may have low power for time-series data (107).

The process of filtering combinations of five cohort parameter values to determine overwintering cohort settings resulted in 22 parameter combinations ( Supplementary Table S3 ). The lowest combined MAE for first and peak adult emergence (i.e., sum of MAE for both events) was obtained for the model that used a low bound = 60, mean = 200, high bound = 350, and variance = 15000 of development completion times [in degree-days (°C), DDC] for J-larvae. Accordingly, adult emergence was predicted to occur between 217 and 460 DDCs ( Supplementary Table S6 and Figure S4 ). Predicted DOYs exhibited a linear relationship with observed DOYs ( Figure 3 ), although the line of fit for peak adult emergence curved downward at higher DOYs. The MAE between predicted and observed DOYs for first and peak adult emergence were 3.7 and 4.6 days, respectively ( Table 3 ). On average, first adult emergence was overpredicted by 1.2 days [standard deviation (SD) = 4.8] whereas peak adult emergence was underpredicted by 0.8 days (SD = 6.1).

According to model validation analyses, predicted DOYs for phenological events in adults exhibited mostly linear relationships with observed DOYs ( Table 3 ; Figure 4 ). However, a single observation from Washington DC (34) resulted in a slight upward curve at lower DOYs for peak adult activity, whereas a single observation from Louisiana (80) resulted in a downward curve for peak adult emergence at lower DOYs. MAE values for first adult emergence, peak adult emergence, and peak adult activity were 7.4, 8.4, and 7.2 days, respectively. On average, first adult emergence was overpredicted by 1.6 days (SD = 9.2) whereas peak adult emergence and peak adult activity were underpredicted by 0.1 (SD = 10.2) and 0.9 days (SD = 8.2), respectively. Predicted DOYs for first pupation and first egg hatch exhibited lower correspondence with observed DOYs ( Figure 4 ) and were less accurate (MAE = 12.0 and 13.7 days, respectively) than predictions of adult events. However, sample sizes for these two events were very small (N = 3 observations each).

According to equivalence tests for first adult emergence and peak adult activity, the two one-sided 95% CIs for difference between the means of predicted and observed DOYs were completely contained in the equivalence interval of 14 days ( Figure 5 ), indicating that the means were statistically equivalent (no bias). However, the CIs for the difference between the means of predicted and observed DOYs for first adult emergence exceeded the upper equivalence bound of seven days, indicating overprediction of up to seven days. The CIs for the difference between the means of predicted and observed DOYs for peak adult emergence exceeded both the lower and upper bounds of the equivalence interval of 14 days, indicating under- and over-prediction of up to 14 days.

The climatic suitability model correctly predicted presence of EAB for 99.9% (3024/3028) of presence records used for model validation, which provides evidence for high model sensitivity. The potential distribution based on climate data for 20 years included most of China ( Figure 6A ), and it overlapped with the ranges of all native Fraxinus species in North America ( Figure 6B ) and Europe ( Figure 6C ). In North America, areas predicted to be excluded by cold stress for at least one modeled year corresponded roughly with the northern range edges (ca. north of 49 °N in continental areas) of green ash and black ash, F. pennsylvanica Marsh and F. nigra Marsh, respectively ( Figure 6B ). However, a few areas along these range edges were excluded by moderate cold stress during an extreme year in terms of cold stress accumulation (1999; Figure 7A ). Conversely, cold stress accumulation in Europe did not exceed EAB’s limits in any region, except for some small areas of the northernmost latitudes of Europe and European Russia in 1999 ( Figure 8A ).

Heat stress played a minor role in shaping EAB’s potential distribution within the ranges of native Fraxinus species in North America and Europe ( Figures 7B , 8B ). However, parts of the ranges of Mexican ash (F. berlandieriana DC) and littleleaf ash (F. greggii Gray) in northeastern Mexico as well as California ash (F. dipetala Hook. & Arn.) in Baja California were excluded by moderate heat stress ( Figure 7B ). Areas excluded by moderate and severe heat stresses expanded throughout this area and in southern United States (Texas and Oklahoma) during an extreme year in terms of heat stress accumulation (2011).

Predictions of first adult emergence and egg hatch for EAB in 2021 varied substantially by latitude. For most of Europe and northern parts of North America, first adult emergence occurred between June and July ( Figures 9A , 10A ) whereas egg hatch occurred between late June and August ( Figures 9B , 10B ). The earliest dates for these events in North America were in southwestern areas (Mexico, southern Texas, and Florida) and latest in southeastern Canada, whereas in Europe the earliest dates were in Spain and latest in the United Kingdom and Scandinavia. There was insufficient degree-day accumulation for egg hatch to occur in parts of the ranges of Fraxinus species in southern Canada and northern Europe (e.g., in Ireland, Great Britain, and Scandinavia), which suggests that EAB may be unable to complete its life cycle in these areas. Within the host ranges, moderate cold stress was predicted for a small area of southern Saskatchewan and Manitoba in Canada.

Trends in the DOY of first adult emergence over 20 recent years were significant in only a few areas. Significant trends towards earlier emergence in North America occurred in both southern areas (parts of Mexico and the southeastern United States) and northern areas (parts of New England and the Midwest; Figure 11A ), whereas in Europe they occurred predominantly in central-eastern regions (Ukraine and European Russia; Figure 11B ). A significant trend towards later adult emergence were predicted in a few small areas of the midwestern United States and eastern Turkey.

Significant trends in range-limiting climate stress for EAB in North America between 2002 and 2021 occurred predominantly in areas outside of the range of native Fraxinus species. Areas excluded by cold stress declined in parts of Quebec, the Hudson Bay area in central Canada, and northern Alberta, whereas they increased in parts of western Canada ( Figure 12A ). Areas excluded by heat stress significantly increased in parts of northern Mexico and the southwestern United States, including in parts of the ranges of F. berlandieriana and F. greggii in northeastern Mexico and F. dipetala in Baja California ( Figure 12B ).

According to model validation analyses, predicted DOYs for phenological events in adults exhibited strong concordance with observed DOYs, with mean absolute errors of ca. 7 days and low estimates of bias. However, an equivalence test indicated that the mean predicted DOY of first adult emergence was overpredicted by as much as seven days, which is concerning because surveillance and management tactics that target adults are often employed several days before adult emergence (5, 30). As a precautionary measure, activities such as installing detection devices and applying insecticide treatments could be implemented at least a week before the predicted date of first adult emergence to help avoid the potential consequences of model overprediction (e.g., failing to detect beetles that had already emerged).

Additional phenological observations collected from across a latitudinal range are needed to conduct a more robust evaluation of model performance. For example, high levels of under- and over-prediction of the DOY of peak adult emergence (±14 days) according to equivalence tests may be partly due to low statistical power owing to a small sample size (N = 8 observations). Sample sizes were particularly small for first pupation and first egg hatch (N = 3 observations each), so we refrained from interpreting results of analyses for these events. Additionally, observations of all life stages are needed from southern populations to better evaluate model performance in warmer climates. EAB’s activity and development in areas with warm climates are not well understood because the pest began spreading into the southern United States only over the past decade (14, 108, 109). To date, monitoring data collected from the southern United States indicate that EAB overwinters primarily as J-larvae and that populations are univoltine (78–80, 84). Thus, model assumptions of a one-year life cycle in which the J-larva is the overwintering stage are likely valid for this region.

The model could potentially be improved by incorporating stage-specific temperature thresholds and a low temperature requirement for the completion of diapause development in the J-larval stage. Using different thresholds for immature stages (10°C for larvae and pupae, and 13.7°C for eggs) did not significantly affect prediction error, but the analysis was based on a limited observation dataset. Another area needing further study is whether termination or completion of diapause development in J-larvae requires a specific low temperature regime or other cue (110, 111). Recent work indicated that J-larvae may require at least two months at moderately low temperatures (ca. 12.8°C) to complete diapause development before pupation can begin (36, 64). However, future work should investigate whether J-larvae in southern populations have different thermal requirements to terminate diapause. A low temperature requirement for diapause completion was not added to the model in part because of DDRP’s single-year modeling structure, which prevents it from using climate data for the entire winter (e.g., November–February). The single-year modeling structure also prevents modeling a two-year life cycle for EAB.

In CONUS, trends towards earlier adult emergence over 20 recent years were predicted in some areas that experienced marked increases in average winter (December−February) temperatures (≥ 1.7°C) and spring (March−May) temperatures (0.5−3.4°C) since 1970, such as parts of New England and the South (112, 113). Similarly, trends in Europe occurred predominantly in areas (European Russia and Ukraine) that experienced rapid increases in winter and spring temperatures (0.07−0.08°C/year) between 1985 and 2020 (50). However, temporal trends in the DOY of first adult emergence were not significant for most areas of North America and Europe, which may be due to our use of a rather short (20 years) time-series dataset that likely exhibited serial correlation (i.e., the variable is correlated with itself over periods of time) (107). Alternatively, climate warming may not be impacting the timing of adult emergence as much as other phenological events. For example, a degree-day modeling study of EAB in Winnipeg, Canada concluded that climate change between 1970 and 2010 significantly impacted the date of peak activity but not the date of emergence in adults (73).

The potential distribution of EAB based on 20-year climate averages overlapped with the ranges of all native Fraxinus species in North America and Europe, which is broadly consistent with previous climatic suitability models for EAB based on climate normals (41–47). Very high estimates of sensitivity of the climatic suitability model (99% of known presence records correctly modeled) suggest that it does not underpredict the potential distribution of EAB. However, we were unable to quantify model specificity (i.e., the proportion of known absence records correctly modeled) owing to a lack of data on absences of EAB in relation to climate. Thus, absence records are needed to assess whether the model overpredicts the potential distribution and to quantify overall model accuracy.

For certain years, cold stress excluded populations at the northern range of Fraxinus species in North America (i.e., green ash and black ash), which may be explained by severe cold events that can reduce overwintering survival rates in EAB (48, 54). While North America and Europe are experiencing warming winter temperatures and an overall decrease in the frequency and duration of extreme cold events (52, 106, 114), the frequency of high intensity cold events increased in some areas between 1950 and 2020 (115). Spatiotemporal variation in extreme cold events may explain why range-limiting cold stress for EAB between 2002 and 2021 in North America significantly decreased in some areas (e.g., in Quebec) but increased in others (e.g., in British Columbia). Thus, the impacts of climate change on EAB’s survival in the coldest parts of its distribution will likely vary across space and time.

Models based on 20-year climate averages indicated a minor role for heat stress in shaping EAB’s potential distribution in areas with native Fraxinus species. However, summer temperatures and the intensity, frequency, and duration of heat waves are increasing in western North America and are projected to worsen under enhanced global warming (49, 106, 116). Thus, the ranges of certain host species in southwestern North America such as F. berlandieriana, F. greggii, and F. dipetala could potentially become increasingly unsuitable under enhanced warming. Conversely, range-limiting heat stress for EAB did not significantly increase in the southeastern United States. Previous studies likely underestimated climatic suitability for EAB in this region because they used correlative algorithms and did not completely sample the species’ climatic distribution (42, 43).

Climatic suitability modeling in DDRP is based only on survival-limiting climate stresses, but the ability for a pest to complete its life cycle in a new climate also influences establishment risk [e.g.,(40)]. Insufficient degree-day accumulation for egg hatch to occur in parts of southern Canada and northern Europe during 2021, which had record-breaking summer temperatures, suggests that EAB is incapable of completing its life cycle in these areas. Similarly, models based on annual growing degree-days for Europe predicted that EAB could not complete development in many areas of Norway, Sweden, Finland, Ireland, and Great Britain (7, 46). Overlaying DDRP predictions of egg hatch produced for multiple years may provide greater insight into spatiotemporal variability in establishment risk based on life cycle completion. Future work could also investigate the relative roles of cold stress vs. incomplete life cycle progression in potentially excluding EAB from parts of southern Canada.

Differences in monitoring methods used across studies likely explain at least some prediction error for adult phenological events. For example, some studies counted the number of exit holes to monitor adult emergence whereas others used trapping data. Adults are capable of immediate flight upon emergence (16); however, the efficiency in counting adults with traps and visual observation of adult exit holes may differ. While trap captures of EAB adults are affected by trap deployment location and many other factors such as beetle age and weather, visual observation of adult exit holes on infested ash trees may be less affected by those factors and thus may provide more accurate counts of emerging adults (J. Duan, unpublished data).

An additional explanation for prediction error for adult phenological events is our use of seven cohorts and a single set of cohort parameter values for all populations. These settings produced a distribution in adult emergence times that corresponded well with field data used for model calibration; however, distributions likely vary across the range of EAB. For example, decreasing the low bound for the timing of the completion of J-larvae (xdist1) may be appropriate for populations that have cohorts which require fewer degree days to complete J-larval development. The phenology model could be modified to accept a raster of values for cohort numbers and cohort parameter values; however, significantly more observations of pupal development and adult emergence of overwintered insects across latitudinal gradients are needed to estimate how J-larval development varies within populations across EAB’s range.

Differences between air temperatures and under-bark microclimates experienced by immature stages of EAB are a likely source of model error (32, 41, 117). Under-bark temperatures may be 1−7°C higher than ambient air temperatures during winter owing to the buffering capacity of snow and tree bark, potentially increasing rates of development and overwintering survival for EAB (117, 118). However, accounting for these temperature differences in our model is difficult because making assumptions of a constant level of thermal buffering are not valid owing to variability both within and between trees (117). For example, urban heating (45, 117) as well as variation in solar insolation within a forest patch (e.g., density of canopy cover) and within individual trees (e.g., the south side of a tree experiences warmer temperatures because it faces sunlight) can affect under-bark microclimates experienced by EAB (37, 66, 117, 118).

Our model assumes that tolerance to climate stress is constant across EAB’s range, but cold tolerance in EAB is a phenotypically plastic trait that varies both temporally and geographically (118–120). Mid-winter warm spells followed by an extreme cold event may have lethal effects on non-diapaused overwintering larvae owing to the breakdown of cryoprotectants during deacclimation, which decreases their cold tolerance (120). Populations in the northernmost (coldest) parts of EAB’s range may have higher survival rates during extreme cold events compared to southern populations because they acclimate to colder winter temperatures (119). A mechanistic model of overwintering mortality for EAB incorporated an equation to model the temperature dose-response relationship for overwintering J-larvae (48). However, the single-year modeling structure of DDRP would likely hinder using a similar method for the EAB model.

Biological factors not included in the model that may affect EAB’s development and survival include host tree health and nutritional quality of host tissues (66, 83, 121), presence of parasitoids and predators (122, 123), levels of host plant resistance (66, 124), and density of infestations (83). Phenology within a site may vary over time owing to host composition. Upon establishment, EAB kills the most susceptible ash species (green and black ash) first and then infests more resistant ash species (e.g., white ash, Fraxinus americana L.) and even white fringe trees (Chionanthus virginicus L.), which may result in delays in adult emergence in later years because of poorer host tree nutrition or other detrimental factors on larval development (125–127).

The DDRP model for EAB has been operationalized for real-time decision support for CONUS at USPest.org (https://uspest.org/CAPS) and at the USA National Phenology Network (https://usanpn.org/data/forecasts/EAB). USPest.org provides all model outputs including phenological event maps for all life stages in raster (GeoTIFF) and summary map (PNG) formats. Conversely, the USA National Phenology Network presents predictions of first adult emergence and egg hatch as Pheno Forecasts (38) in summary map and interactive formats, and end users can sign up to receive e-mail notifications that provide advanced warnings (3, 2, and 1-week) of when these events will occur in their area. Interactive forecasts allow end users to zoom, pan, and interrogate maps (e.g., click on individual raster pixels in an area of interest). Pheno Forecasts for EAB show predictions only in climatically suitable areas and convert dates of phenological events to time relative to the map issue date (e.g., “This Week”, “Next Week”, etc.), which is a preferred format of many Pheno Forecast end users (38). Model outputs at both web sites are updated every three days. Climate datasets used for real-time modeling include PRISM daily data and daily-downscaled NMME (North American Multi-Model Ensemble) 7-month forecasts (128) at a spatial resolution of 4 km² (55).

Real-time forecasts of adult emergence can support timely surveillance of adults, which are the most visible and therefore most easily detected life stage of the pest. For example, notifications that adults will emerge in three weeks can alert surveillance teams to finish installing detection devices such as sticky prism traps or funnel traps before adults begin flying in the canopy. Early detection of EAB is critical because infestations in ash trees are usually fatal, and failure to detect adults early may result in dispersal to new locations and increase the cost of ash treatments and removal (31, 129). Thus, early detection may help slow the spread of EAB in new regions such as the west coast of North America, where the pest was detected for the first time in 2022 in Forest Grove, Oregon (26). Forecasts of adult emergence may also improve the timing of systemic insecticide treatments and cover sprays that target adults feeding on host foliage (5, 130). For example, a manager can begin applying systemic insecticides when they know to expect adults in 2−3 weeks because females must feed on leaves for ca. 1−2 weeks before they begin laying eggs (5).

Real-time forecasts of egg hatch may help with the timing of insecticide applications that target EAB larvae, which are most effective when newly hatched larvae encounter the insecticide as they chew through the bark and into the cambial tissue (5, 130). Additionally, forecasts of adult oviposition and egg hatch may help biocontrol practitioners identify the critical periods for releases of biological control agents against EAB, which may increase the likelihood that agents establish and grow in the target area. Releases of the egg parasitoid Oobius agrili Zhang and Huang (Hymenoptera: Encyrtidae) must be timed during the EAB oviposition period because eggs are only suitable for parasitism up to the development of the neonate host larva (30, 131). Conversely, the larval parasitoids such as Spathius agrili Yang (Hymenoptera: Braconidae) and S. galinae Belokobylskij & Strazanac (Hymenoptera: Braconidae) must be released when 3^rd and 4^th instar larvae are present (131, 132).

While the DDRP model for EAB exhibited overall good performance, it should still be used conservatively for decision-making, particularly given that model validation analyses were based on a small observation dataset. As mentioned previously, surveillance and management activities could be implemented several days (e.g., a week) prior to the date of a predicted event. With additional observations of EAB, one could assess how to interpret forecasts in order to maximize the probability that adults are captured upon first emergence or that insecticides or biocontrol treatments are employed within the optimal window of time. To potentially avoid under-predicting the risk of establishment, the potential distribution could be defined as areas not under severe climate stress as opposed to defining it using both stress levels.

A lack of near real-time daily climate datasets for Canada and Europe would likely hinder using the model for real-time decision support in these regions despite a need for additional tools to manage EAB (12, 133). Continued spread of EAB in Canada is predicted to cause substantial economic impacts primarily owing to its destruction of urban forests (17). In Europe, the pest has spread to 16 regions of European Russia and to the east of Ukraine (19, 21, 28) and is expected to expand into neighboring countries within 5−20 years (27, 134). Declines in Europe’s most widespread ash species, European ash (F. excelsior L.), would have severe ecological and economic consequences because this species is a key component of many forests and has been widely planted in cities, parks and along roads as shade or ornamental trees (12, 21, 134, 135).

Future work on the EAB model could add capabilities that would allow it to predict EAB’s spread across the landscape and forecast the timing of its arrival (136, 137), which would provide a more comprehensive assessment of both where and when to expect this pest. For example, surveillance could be focused on areas that are climatically suitable for EAB and have a high likelihood of being colonized in the next five years. Within areas at highest risk of establishment, forecasts of adult emergence over subsequent years may help ensure that surveillance teams install detection devices on time.