Ips typographus and Dendroctonus ponderosae Models Project Thermal Suitability for Intra- and Inter-Continental Establishment in a Changing Climate

Introduction

Living forests are critical components of the earth's climate system. Growing trees absorb carbon dioxide (CO₂) and can store it for hundreds to thousands of years. Increasing atmospheric concentrations of CO₂, largely caused by industrial activities, are contributing to global climate change (IPCC, 2013) and forest carbon sequestration and storage have become important mitigation strategies (Dugan et al., 2018; Xu et al., 2018). Multiple factors affect in situ forest carbon storage capacity including tree growth and decay rates, management, and tree mortality due to deforestation and natural disturbances (Breshears and Allen, 2002 Li et al., 2003; Canadell and Raupach, 2008). Insects are significant tree mortality agents globally, and bark beetles (Coleoptera: Curculionidae, Scolytinae) are among the most important in forest ecosystems (Hicke et al., 2013; Mezei et al., 2017). Bark beetles are integral ecosystem components that promote heterogeneity and resilience across multiple scales (Kulakowski et al., 2017), yet can also impact economic, societal (Grégoire et al., 2015) and short-term carbon sequestration (Pfeifer et al., 2011; Hansen et al., 2014) goals. The bark beetle and forest carbon feedback to climate change highlights the importance of understanding future climate influences on global distributions of bark beetle populations (Kurz et al., 2008; Raffa et al., 2008; Adams et al., 2010). Climatic changes have already influenced patterns and timing of bark beetle population outbreaks in their native ecosystems (Buotte et al., 2016; Senf and Seidl, 2018), in addition to intra- and inter-continental range expansion and invasion of novel habitats by some species (Liu et al., 2014; Cooke and Carroll, 2017). Predicting tree-killing bark beetle population responses, including range expansion and potential contraction, to future climate projections is an important component of forest management and carbon sequestration strategies globally.

The European Ips typographus (L.) and North American (NA) Dendroctonus ponderosae (Hopkins) are considered among the most important disturbances affecting tree mortality on each continent. Due to thermally-dependent fitness traits and positive responses to drought-induced tree stress, population activity has increased recently with a changing climate (Weed et al., 2015; Kolb et al., 2016; Marini et al., 2017; Senf and Seidl, 2018). In western Canada and the United States (US) recent tree mortality due to D. ponderosae exceeded 28 Mha (Hicke et al., 2015; Cooke and Carroll, 2017). In Europe tree mortality caused by bark beetles, mainly I. typographus, exceeded 2.88 million m³ per year during the period 1958–2001 and has increased in recent decades (Schelhaas et al., 2003; Seidl et al., 2011).

Both species feed and reproduce in tree phloem, typically killing the host when at epidemic levels, and they select mature trees within the appropriate genera that are of sufficient size and bark/phloem thickness for offspring production. Both species have relatively large geographic distributions that follow that of their primary host trees, Picea abies for I. typographus and Pinus species for D. ponderosae. Across their vast ranges, phenological plasticity in thermal- and photoperiod-regulated traits allow for adjustments in lifecycle, and hence generation timing, depending on year-round available heat, day length, and timing of cold temperatures that can induce offspring death (Bentz et al., 2014; Schebeck et al., 2017; Schroeder and Dalin, 2017). Many of these physiological traits are directly influenced by climatic changes and can differ between the species (Bentz and Jönsson, 2015). Moreover, although both species have adapted to feed and breed in living tree tissues, their tolerance to active plant defenses and density dependent traits also differ (Kärvemo and Schroeder, 2010).

The distribution of I. typographus follows the entirety of its host tree Picea abies which has a continuous range in Scandinavia, north-eastern Europe and western Russia, and at higher latitudes in central Europe (Farjon and Filer, 2013). The current distribution of D. ponderosae, in contrast, is not as extensive as species within its host genus Pinus, which are contiguous further north and east in Canada, further east in the US, and further south in Mexico than the current distribution of D. ponderosae. Canada was covered in ice 21,000 years BP, and as glaciers retreated several Pinus species colonized the new habitat (Godbout et al., 2008; Roberts and Hamann, 2015). Despite the availability of Pinus host trees, the northernmost areas remained thermally unsuitable for D. ponderosae (Safranyik, 1978). Recent warming in western Canada resulted in a rapid northward expansion of D. ponderosae into Pinus habitat of northern British Columbia and western Alberta (de le Giroday et al., 2012). Climate change-induced range expansion of D. ponderosae southward has not occurred, despite the availability of known Pinus hosts. Multiple Pinus species are found at low elevations in Arizona, Mexico, and Central America, yet the most southern locations where D. ponderosae has been observed is high elevation forests of southern Arizona, US, northern Chihuahua, Mexico, and northern Baja California Norte, Mexico (Armendáriz-Toledano et al., 2017; Bentz and Hansen, 2017; Dowle et al., 2017).

Given the demand and volume of softwood timber products in the global market there is a risk of inter-continental introduction of both species, and fitness traits that influence population success are known to be useful in predicting invasion potential (Ward and Masters, 2007). Following introduction, establishment is dependent on (1) propagule pressure, which may be estimated as frequency of individual interception at ports, (2) host tree availability, and (3) climate suitability. Ips typographus has been intercepted almost 500 times in NA ports, although establishment has not been recorded (Liebhold et al., 2017). There are no records of D. ponderosae interception in Europe. Ips typographus and D. ponderosae can successfully reproduce in Picea and Pinus species not found on their native continents (Økland et al., 2011; Fries, 2017; Rosenberger et al., 2017; Flø et al., 2018), suggesting that if introductions were to occur, suitable trees for reproduction on either continent would not be a barrier to establishment and spread. A significant unknown is the climate suitability for I. typographus in NA, and D. ponderosae in Europe. Evolved, thermally-dependent traits in both species affect population success and voltinism (i.e., number of lifecycles or generations in a year), and thereby the potential for establishment and spread.

Our goal was to describe potential changes in the distribution and voltinism of I. typographus and D. ponderosae given current and projected temperatures across NA and Europe. We used previously published phenology models for predicting lifecycle timing and number of annual generations based on thermally-influenced fitness traits specific to each species. By running our models for each species on each continent we assessed potential shifts in voltinism and distribution in their native continents, in addition to thermal suitability for D. ponderosae in Europe and I. typographus in NA. Model results provide measures for predicting establishment and persistence of each species if invasion should occur in non-native habitats and continents.

Materials and Methods

Study System: A Comparison of I. typographus and D. ponderosae Life History Strategies

Ips typographus and D. ponderosae use aggregation pheromones during the colonization of live host material which enables a rapid mass attack (i.e., within days) that can overwhelm host defenses (Christiansen and Bakke, 1988; Raffa et al., 2008). Male I. typographus initiate attacks and they are generally accompanied by 1–3 females (i.e., polygamous), whereas in D. ponderosae females initiate attacks and both males and females can have multiple mates (i.e., polyandry and polygyny), at least during epidemic phases (Janes et al., 2016). Both species have fungal associates that can assist in compromising host tree defenses and for D. ponderosae they also provide nutrition for developing brood (Hofstetter et al., 2015; Krokene, 2015; Zhao et al., 2018). Dendroctonus ponderosae is obligatory dependent on several fungal associates that are carried in mycangia, highly specialized structures of the exoskeleton (Six and Bracewell, 2015). Mycangia have not been identified in I. typographus and mutualistic fungal associates are unclear. Although both species are capable of infesting and reproducing in live standing trees, I. typographus is commonly found in wind-felled trees with no or weak defenses (Christiansen and Bakke, 1988; Schroeder, 2010) and storms are therefore considered an important factor in population outbreaks (Kärvemo et al., 2014; Marini et al., 2017; Potterf and Bone, 2017). Dendroctonus ponderosae rarely colonizes downed trees and has evolved semiochemical feedback mechanisms whereby high population density facilitates high population growth (i.e., positive density dependence) when live standing trees are attacked (Raffa et al., 2008). In live standing trees the average reproductive success of D. ponderosae is approximately 3 times higher than for I. typographus, in part due to significantly lower attack densities and therefore less intraspecific competition in D. ponderosae relative to I. typographus (Anderbrant et al., 1985; Kärvemo and Schroeder, 2010; Komonen et al., 2011). Outbreaks of both species can be initiated when drought stresses large landscapes of host trees, particularly when associated with warm temperatures (Kolb et al., 2016; Marini et al., 2017), allowing smaller numbers of attacking beetles to overcome the weakened defenses (Mulock and Christiansen, 1986; Raffa et al., 2008).

Following attack and oviposition parent beetles of both species can emerge to establish an additional cohort, known as a sister brood, in nearby trees although this trait is more common in I. typographus (Reid, 1962; Davídková and Doležal, 2017). Ips typographus is generally univoltine (i.e., one generation annually) across all of Scandinavia except Denmark, bivoltine (i.e., two generations annually) in Denmark, and bi- or multivoltine in areas further south except at high elevations (Christiansen and Bakke, 1988; Wermelinger, 2004). Dendroctonus ponderosae is univoltine across the majority of its range in western NA with semivoltinism (i.e., 2 years for a single generation) at the highest elevations and rare to no occurrences of bivoltinism (Bentz and Powell, 2014; Bentz et al., 2014).

Larval stages of I. typographus are not cold tolerant and the adult stage is the typical overwintering lifestage (Koštál et al., 2011; Dworschak et al., 2014). In D. ponderosae the last larval stage is most cold tolerant and also the typical overwintering lifestage, although other lifestages can survive warm winters (Bentz and Mullins, 1999; Lester and Irwin, 2012; Rosenberger et al., 2017). Both species have an arrested development, known as diapause, although the lifestage where it occurs differs. Diapause in I. typographus occurs in the adult stage and serves to prevent reproduction in late summer or fall that would result in larvae during winter (Doležal and Sehnal, 2007). The critical day length for induction of adult reproductive diapause varies latitudinally with southern populations entering diapause at shorter day lengths than northern populations (Schroeder and Dalin, 2017). The most northern populations also include a considerable proportion of beetles with an obligatory diapause wherein diapause occurs regardless of environmental conditions. Dendroctonus ponderosae has a facultative diapause that occurs in the last larval stage prior to pupation (i.e., prepupal) that is induced predominately by low temperatures and serves to prevent molting to the cold intolerant pupal stage. Diapause in D. ponderosae also varies latitudinally with a greater proportion of northern US populations entering diapause at warmer temperatures compared to southern US populations (Bentz and Hansen, 2017).

Climate Data

Climate data used in this study were selected from regionally downscaled global climate model (GCM) scenarios for NA and Europe. The data were produced within the CORDEX project (Giorgi and Gutowski, 2015, http://www.cordex.org/) on a grid-spacing of 0.44° × 0.44°(~50 × 50 km), covering the period of 1971–2100 on a daily temporal resolution. To account for uncertainties in climate model simulations, data of two GCMs were selected, CCCma-CanESM2 and ICHEC-EC-EARTH. Projections from these two GCMs are within the mean temperature and precipitation responses of other GCMs within CMIP5. Both data sets were downscaled with the regional model SMHI-RCA4 and bias adjusted using an empirical quantile mapping method (Themeßl et al., 2011, Wilcke et al., 2013). For the European domain, bias adjustment was based on high quality observation data from MESAN-EURO4M (Bärring et al., 2014; Landelius et al., 2016) for the period 1989–2010. For NA, re-analysis data from WFDEI for the 1979–2012 period (Weedon et al., 2014) was used.

We selected climate model runs based on the representative concentration pathway, RCP 8.5, that characterizes a future with relatively high greenhouse gas emissions (Moss et al., 2010; IPCC, 2013). RCP 8.5 represents the highest increase in global mean temperature by the end of the century. Temperatures expected under RCP 2.6 at the end of the century are similar to present day climate and RCP 4.5 and 6.0 expected temperatures by the end of the century are similar to mid-century RCP 8.5 temperatures. Patterns resulting from all RCPs can therefore be extracted from runs based on RCP 8.5 by examining output at multiple time frames. Daily maximum and minimum temperature climate simulation data from 1981 to 2100 were used to run I. typographus and D. ponderosae physiological models, and output from multiple time periods throughout the century therefore represent multiple RCP scenarios.

Ips typographus Physiological Model

Thermal response of I. typographus was modeled using a previously published phenology model that was developed to simulate large-scale trends of swarming (i.e., adult host-seeking and dispersal) and lifecycle development based on gridded climate data with low spatial resolution (50 × 50 km) (Jönsson et al., 2011; Figure 1). Therefore, the model is more generalized than concepts developed to simulate local and fine-scale topoclimatic conditions (Baier et al., 2007). Thermal sums are expressed as degree-days (dd) above a developmental threshold of +5°C for all developmental stages (Annila, 1969), and emergence from winter hibernation was set to 120 dd. Flight in search of suitable breeding material is related to a temperature threshold, calibrated to +16°C using data on pheromone trap catches in combination with gridded temperature data (Jönsson et al., 2011). Egg development was assumed to start on the seventh day after swarming, accounting for a pre-oviposition period and sufficient time to oviposit 50% of eggs (Jönsson et al., 2009). Variation in microclimate within forest stands, in particular differences caused by sun exposure on tree boles, creates a continuum of breeding conditions. Because I. typographus prefers sun-exposed over shaded substrates, model parameters included in this study were originally derived from field data using exposed boles, with developmental time from egg to mature adult corresponding to a thermal sum of 625 dd (Harding and Ravn, 1985; Jönsson et al., 2007). Populations are strongly regulated by availability of weakened and downed host trees following storms and sanitary and salvage logging, and the phenology model does not explicitly consider factors such as availability of brood material and reproductive success (Jönsson et al., 2012). Furthermore, as I. typographus adults are cold hardy to −20 to −22°C (Koštál et al., 2011), cold-induced mortality is not considered to be a factor that could influence I. typographus distribution.

Figure 1. (A) Three Dendroctonus ponderosae models were integrated to simulate population response to temperature for individual climate grid cells. The phenology model (MPB-Phen) (Régnière et al., 2012) describes temperature-dependent lifestage development and adult emergence timing from a brood tree. The cold-induced mortality model (MPB-Cold) (Régnière and Bentz, 2007) predicts survival as a function of annual temperatures and their effects on larval supercooling. Both models were driven by daily maximum and minimum temperatures projected by two Global Climate Models for each year between 1981 and 2099. MPB-Phen produces a temporal distribution of adult emergence that was used by the demographic model (MPB-R) (Powell and Bentz, 2009) to predict univoltine population growth (R). MPB-R assumes that adult emergence timing is critical in the acquisition of breeding material and that emergence must be synchronous and appropriately timed in the summer. Synchronous emergence is quantified as the number of emerged adults that exceed a threshold (A) required to overwhelm host tree defenses. Adults that emerge on a given day between 30 June and 1 September in sufficient numbers to exceed A are considered effective (E). Parameters α, β, and H were estimated using field-collected temperatures and data on hectares killed by D. ponderosae (Powell and Bentz, 2009, 2013). Univoltine population growth (R) was modified by predicted larval survival based on MPB-Cold. MPB-Phen was also used to evaluate if temperatures across a given year were sufficient to produce two generations in a single year (bivoltine). (B) The Ips typographus phenology model (Jönsson et al., 2012) uses daily maximum and mean temperatures to predict the occurrence of adult diapause, which dictates either the univoltine or bivoltine lifecycle. Ips typographus population growth is highly dependent on the stochastic occurrence of storm felled and drought-stressed trees, which are also regulated by forest management. Population growth of I. typographus was therefore not modelled for this analyses.

Reproductive diapause of the first generation was modeled to occur after fulfillment of day length and temperature (daily mean<15°C) requirements (Jönsson et al., 2011). The day length setting was based on threshold values that vary geographically and accounts for local adaptation to cold autumn temperatures. The gridcell-specific threshold corresponds to the day length of the earliest date during a climate reference period (1961–1990) when I. typographus were not able to reach the adult lifestage required for winter survival. This parameterization, evaluated against independent field monitoring data (Jönsson et al., 2011), follows a latitudinal gradient from approximately 15–16 h in central Europe to 20–22 h in the northern part of Scandinavia. The model captures the day length threshold below which individuals are responsive to lower temperatures that trigger diapause induction. We did not evaluate the potential for a third generation, and therefore reproductive diapause of the second generation was not modeled. Ips typographus model output is summarized as the number of years projected to have one and two generations annually for 30 year periods (1981–2010, 2011–2040, 2041–2070, 2071–2100).

Dendroctonus ponderosae Physiological Models

Thermal response of D. ponderosae was simulated with the integration of three models (Figure 1). A demographic model, MPB-R (Powell and Bentz, 2009), describes annual univoltine population growth (R) as a function of thermally-driven phenology, MPB-Phen (Régnière et al., 2012), which is then modified by predicted larval survival based on a cold tolerance model, MPB-Cold (Régnière and Bentz, 2007). MPB-R indirectly incorporates the role of tree defense in population success by connecting brood adult emergence timing, predicted using MPB-Phen, with population consequences based on adult emergence synchrony required for successful mass attacks on a tree. Live standing trees are the primary host type and defenses are overwhelmed by large numbers of conspecifics attacking in a short period of time (Logan and Bentz, 1999). Thermal regimes across a generation result in a predicted brood adult emergence distribution that is “effective” when emergence between 30 June and 1 September exceeds a daily threshold of 250 attacking beetles required to successfully mass attack non-stressed host trees (Powell and Bentz, 2009). Emergence during this time window is considered seasonally optimal for development to an appropriate overwintering stage (Logan and Bentz, 1999). High R values >1 signify thermal conditions that support highly synchronized and seasonally appropriate adult emergence. Low R values <1 occur when simulated adult emergence occurs either outside the window of 30 June and 1 September or is unsynchronized.

MPB-Phen predicts development and adult emergence timing as a non-linear function of temperature with upper and lower developmental thresholds that vary among lifestages and lognormal genetic variability among individuals. Mechanistically, MPB-Phen is based on the cohort approach and predicts the probability that an individual in a given lifestage will complete development within a given time interval. Prepupal diapause is currently implemented as a high temperature threshold for pupation. Following initiation with a distribution of attacks on a tree, described by a normal distribution with a mean on 24 July and standard deviation of 2 days, hourly temperature inputs drive development through eight lifestages, culminating with an emergence distribution of brood adults which is then used by MPB-R to predict population growth for that year (Figure 1).

As described in Powell and Bentz (2009), unknown parameters in MPB-R were estimated based on non-linear regression between annual aerially-observed D. ponderosae-caused tree mortality in central Idaho, US, and brood adult emergence predictions from MPB-Phen. Observed hourly air temperatures taken from within the impacted stands were used to drive MPB-Phen. MPB-R successfully described 92% of the variance in observed population growth rates when applied to the same stands in central Idaho where parameters were derived, and >90% of variance in novel stands in northern Washington, US (Powell and Bentz, 2013). In a separate study in Montana, US, model-predicted increasing trends in population growth were significant in explaining observed increasing trends in D. ponderosae-caused tree mortality (Bentz et al., 2016).

Successful bivoltinism has not been observed for D. ponderosae in the field (Bentz and Powell, 2014; Bentz et al., 2014; but see Mitton and Ferrenberg, 2012), and conditions that would support this phenological strategy are unclear. Here, we assume that completion of two generations in a year requires similar conditions as for a univoltine lifecycle, and we used MPB-Phen to evaluate potential bivoltinism in a future climate and in Pinus habitat outside the current range of D. ponderosae. MPB-Phen was run over sufficient time to complete two generations for each grid cell of climate simulation data and year, with tests for temperatures that cause mortality of the cold-sensitive egg and pupal lifestages (Logan and Bentz, 1999; Bleiker et al., 2017), resulting in a binary response of thermally unsuitable (0) or thermally suitable (1) for bivoltinism. Because bivoltine populations have not been observed, an MPB-R model for bivoltine populations has not been developed and model predictions have not been field-evaluated.

Cold temperature is considered a significant mortality factor in D. ponderosae population success and a critical limiting factor of the species' distribution (Weed et al., 2015). We used a mechanistic model (MPB-Cold) (Régnière and Bentz, 2007) derived from laboratory estimates of supercooling points of individual larvae that were field-collected during fall, winter and spring (Bentz and Mullins, 1999). Annual probability of population survival was predicted as a function of daily maximum and minimum temperatures at each grid cell. We assumed population growth to be most affected when a cold event results in high mortality and very low offspring survival. If the predicted annual probability of larval cold survival from the MPB-Cold model was <20%, MPB-R predicted population growth was multiplied by the probability of survival, effectively reducing population growth. No effect due to cold was implemented when predicted survival was >20%.

For each grid cell of climate simulation data in NA and Europe, MPB-Phen and MPB-Cold were simulated using hourly temperatures derived from sinusoidal interpolation of daily extrema and Simpson's rule for rate curve integration using five steps per day (Bentz et al., 2016). Daily extrema came from two GCM temperature simulations (CCCma-CanESM2 and ICHEC-EC-EARTH; see Climate Data) for each year from 1981 to 2100. Each D. ponderosae model was run separately for each GCM, and annual MPB-R output was as an average of model output based on the two GCMs. Based on observed and model-predicted R values we considered a threshold value of R > 1 as indicative of positive trends in univoltine population growth (Powell and Bentz, 2009). The number of years in each 30-years period (1981–2010, 2011–2040, 2041–2070, 2071–2100) with R > 1 and bivoltinism = 1 were summed to describe thermal suitability for univoltine population growth and potential for a successful 2nd generation, respectively.

Study Areas and Vegetation Masks

The NA study area included northern Mexico and all of Canada and the US (~20 to 68°N latitude). In Europe the study boundary extended west of Russia from ~35 to 69°N latitude. Outputs of the physiological models were restricted to grid cells of the climate simulations (~50 × 50 km) on both continents estimated as having at least 5% of the grid cell containing either Pinus species for D. ponderosae or Picea species for I. typographus. Areas (1 km resolution) with host tree species were extracted from the North America Land Cover Characteristics spatial database (https://lta.cr.usgs.gov/glcc/na_int.), and from the European Atlas of Forest Tree Species (de Rigo et al., 2016; 1 km resolution). The rotated climate simulation grid cells for the NA and Europe domains (www.cordex.org) were converted into vector polygons, un-rotated and subsequently transformed into the same spatial reference as the two land cover datasets, respectively. For each polygon representing one climate simulation grid cell, the number of land cover grid cells classified as Pinus or Picea were extracted and recalculated into a percentage of the respective polygon area. Elevation was extracted for each grid cell for both climate models by host tree species ranging from ~2 to 2,019 m in Europe (1,634 grid cells) and ~1 to 3,281 m in NA (3,447 grid cells) for Pinus species, and ~2 to 2,170 m in Europe (1,040 grid cells) and ~1 to 3,281 m in NA (3,794 grid cells) for Picea species (Figure 2).

Figure 2. Spatial variation in elevation (m) of (A) Pinus and (B) Picea habitats in North America and the Europe study area based on the North America Land Cover Characteristics spatial database and European Atlas of Forest Tree Species.

Results

Trends in Climate

In the historical 30-years period (1981–2010), Pinus habitats in the Europe study area were on average 2.6 and 6.6°C warmer in summer and winter, respectively, than Pinus habitats within the current range of D. ponderosae in NA (Figures 3C,D). During this same time period, Picea habitats in Europe were on average 2.5°C warmer in summer and 14.2°C warmer in winter than Picea in NA (Figures 3A,B). To assess habitat on each continent projected to experience the greatest thermal change, we calculated climate anomalies that describe the change in mean summer and winter temperature from historical (1981–2010) to the near future (2011–2040) (Figure 4) and end of the century (Figure S1). According to RCP 8.5, by 2011–2040 summer temperatures in Pinus and Picea habitats on both continents will warm the greatest at the highest elevations, and winter temperatures at the highest latitudes (Figure 4; see also Figure 2). In the near future, warming in Picea habitat is projected to be on average 1°C greater in NA than Europe in the winter and 0.3°C greater in the summer months (Figures 4A,B). Warming in Pinus habitats will also be slightly greater in NA compared to Europe (Figures 4C,D). Trends in warming at the end of the century (2071–2100) are spatially similar to the near future (2011–40) but a greater magnitude (Figure S1). Degree days from 1 January, calculated using 6°C as a lower development threshold, also highlights warming in Pinus and Picea habitats on both continents by the end of the century that is specific to development of each beetle species (Figure S2).

Figure 3. Distribution and mean (A) summer temperatures (°C) in Picea habitats (i.e., grid cells of climate simulation data), (B) winter temperatures in Picea habitats, (C) summer temperatures in Pinus habitats, and (D) winter temperatures in Pinus habitats in North America and the Europe study area during 1981–2010. Temperature projections are an average of two Global Climate Models for RCP 8.5. Overlaps between North America and Europe are represented by the darkest color.

Figure 4. Projected mean temperature change between the historical period (1981–2010) and near future (2011–2040) (°C) in North America and the Europe study area according to two Global Climate Models representing RCP 8.5 for (A) Picea habitats in summer, (B) Picea habitats in winter, (C) Pinus habitats in summer, and (D) Pinus habitats in winter.

Trends in I. typographus population success

During the historical period (1981–2010) the calculated pattern of thermal suitability for I. typographus univoltinism in European Picea habitats was similar to the bimodal pattern across latitude predicted for D. ponderosae univoltinism in European Pinus habitats (Figures 5, 6A,C). Picea habitats in the Europe study area between ~42 and 48°N (representing high elevation areas with cool climates) and 55 and 65°N were predicted to have the greatest proportion of years with univoltinism (Figure 6C). This estimate is in agreement with observed historical patterns of voltinism (Christiansen and Bakke, 1988; Wermelinger, 2004). In a warming climate, thermal suitability for univoltinism was projected to shift northward (Figures 5a–d, 6C), and the proportion years with univolitism at latitudes lower than ~65°N was less in 2071–2100 than in all earlier time periods (Figure 6C). Between 2041–2070 and 2071–2100 simulated univoltinism decreased even at the highest latitudes. The projected decline of univoltinism in a warming climate was associated with an increasing proportion of a 2nd generation (Figures 5e,f). 1981–2010 model simulations suggest that 45% of Picea habitat in the Europe study area was suitable for a 2nd I. typographus generation, and this was projected to increase to 69% by the end of the century. The majority of Picea habitat in the Europe study area suitable for I. typographus in 1981–2010 remained suitable throughout the century for either one or two generations annually Figure 5).

Figure 5. Proportion of years in each 30-years period projected to have one I. typographus generation annually in (a) 1981–2010, (b) 2011–2040, (c) 2041–2070, and (d) 2071–2100 for Picea habitats in North America and the Europe study area. Proportion of years where thermal conditions were suitable for a 2nd generation in (e) 1981–2010, (f) 2011–2040, (g) 2041–2070, and (h) 2071–2100 are also shown.

Figure 6. Proportion of years with projected univoltine population growth (R > 1) for D. ponderosae in Pinus habitats in (A) Europe study area and (B) North America. We assumed that model output of R > 1 was indicative of positive D. ponderosae univoltine population growth.Projected proportion of years with I. typographus univoltinism in Picea habitats in (C) Europe study area and (D) North America as a function of latitude. Model projections were averaged for 30-years periods between 1981 and 2100, and are based on D. ponderosae and I. typographus physiological models and two Global Climate Models representing RCP 8.5. A Loess regression was applied to 30-years average model projections across latitudes.

Model simulations suggest high thermal suitability for I. typographus establishment in NA. During 1981–2010, Picea habitats in NA between ~40 and 60°N were projected to have relatively high probability of univoltinism (> 50% of years in the 30-years period) with much of the boreal forest in Canada exceeding 80% of years with univoltinism (Figures 5a, 6D). As projected temperatures warmed in the near future and through mid-century, the proportion years with univoltinism decreased at lower latitudes and increased at higher latitudes including Picea forests in Alaska, US and eastern Canada (Figures 5b,c 6D). During the historical period (1981–2010), a relatively low probability of a 2nd I. typographus generation was predicted for Picea habitats in NA between ~45 and 55°N (Figure 5e). By the end of the century, thermal suitability for a 2nd generation was predicted to be the highest at any time, including Picea habitats in Alaska and much of the boreal forest in Canada (Figures 5h, 6D). In the near future (2011–2040) and throughout the century, all Picea habitats in NA will be thermally suitable for I. typographus.

Trends in D. ponderosae Population Success

Aerial observations during the historical (1981–2010) period indicate that D. ponderosae population activity occurred in the western US and Canada generally from ~31 to 60°N (Meddens et al., 2012). Based on ground observations, population activity was recorded to be low at the lowest latitudes in this range, except at the highest elevations (Bentz et al., 2014). Our coupled MPB-R, MPB-Phen, and MPB-Cold models estimated positive univoltine population growth (R > 1) in the same general area between ~40 and 60°N latitude in the western US and Canada (Figure 7a). The overlap in distributions of observed D. ponderoase population activity during 1981–2010 and modeled univoltine population growth during the same time period highlights the predictive ability of our models, particularly at the southern US range extent (Figure S3). During 1981–2010 the highest thermal suitability was modeled to be along the Rocky and Cascade Mountains from ~41 to 52°N (Figures 6B, 7a). Although the range of D. ponderosae does not currently extend across the boreal forest in central and eastern Canada, Pinus does occur and model simulations suggest positive univoltine suitability in these areas during 1981–2010. Although thermal conditions were conducive to a univoltine lifecycle, low projected population growth in central Canada was due to low cold survival (<20% survival predicted by the MPB-Cold model) (Figure S4a). By the middle of the century (2041–2070), >50% overwintering survival of D. ponderosae was projected across the majority of Pinus habitat in Canada (Figures S4b,c), and by the end of the century even the highest elevations and latitudes in NA were projected to have high overwintering survival (>70%) (Figure S4d). Despite increasing overwintering survival, predicted population growth of univoltine cohorts declined throughout the century across Canada and the western US, and by the end of the century univoltine cohort probability was the greatest at the highest latitudes (> ~60°N) (Figures 6B, 7d). In the near future (2011–2040) 20% of the habitat that supported univoltine cohorts in 1981–2010 was estimated to no longer be thermally suitable for this strategy (Figure 7b), and by the end of the century 77% of habitat that previously supported high univoltine population growth (R) was no longer thermally suitable (Figure 7d). Based on MPB-R model assumptions (Figure 1), simulated univoltine population growth rates (R) are low when thermal conditions are either too warm or too cool for a phenological pathway that leads to synchronous adult emergence that occurs between 30 June and 1 September. The shift north and loss of thermal suitability for univoltinism occured as summer temperatures in NA increased up to 9°C by the end of the century (Figures 4C,D, Figure S1c). During all time periods, Pinus habitats in the eastern US, which are lower elevation than western US forests at the same latitudes (Figure 2A) and where D. ponderosae is not currently found, were estimated to be thermally unsuitable for successful growth of univoltine cohorts (Figures 7a–d).

Figure 7. Proportion of years in each 30-years period with projected D. ponderosae univoltine (1 year annually) population growth (R > 1) in (a) 1981–2010, (b) 2011–2040, (c) 2041–2070, and (d) 2071–2100 for Pinus habitats in NA and the Europe study area according to two Global Climate Models representing RCP 8.5. Predictions were based on the integrated MPB-R and MPB-Cold models. Dark grey areas represent range contraction, where thermal conditions are not suitable for either univoltinsim or bivoltinism (a 2nd generation occurs annually), relative to thermal suitability in 1981–2010. Proportion of years projected to have thermal suitability for a 2nd generation in (e) 1981–2010, (f) 2011–2040, (g) 2041–2070, and (h) 2071–2100 are also shown.

Although D. ponderosae is not currently found in Europe, thermal suitability for univoltine population growth during the 1981–2010 period would have occurred in Pinus habitats > ~41°N, with the highest thermal suitability between ~55 and 65°N in parts of Sweden and Finland (Figures 6A, 7a). MPB-Cold predicted moderate to high probability of survival across Europe during 1981–2010, with very few locations experiencing low (i.e., <20%) survival events, and by 2071–2100 very high winter survival was simulated in even the most northern locations (Figure S4d). In the near future (2011–2040) thermal conditions that promote univoltine population growth (R) was projected to be elevated at the highest elevations and latitudes in the Europe study area (Figures 6A, 7b). A continual shift northward was projected, and by the end of the century D. ponderosae univoltinism in Pinus habitat in the Europe study area diminished across all latitudes < ~65°N. Between 1981–2010 and 2071–2100, climatic conditions in Pinus habitat in Europe that was conducive to univoltine population growth (R) was projected to decline by 85% (Figure 7d).

Because bivoltine D. ponderosae populations have been rarely observed, an R model that is based on field observations and comparable to that developed for univoltine poulations (MPB-R) has not been developed for bivoltine populations. We used MPB-Phen, developed for univoltine populations, to describe years and habitats where thermal conditions support two generations in a single year. Modeling suggests that thermal suitability for seasonally appropriate 2nd generations (bivoltine) of D. ponderosae occurred during the 1981–2010 period in Pinus habitat further south into Mexico and further east in the southern US than the current D. ponderosae range in NA (Figure 7e). During this same time period in Europe, a high probability of a 2nd D. ponderosae generation was predicted for Pinus in Spain and Portugal and other low latitudes of the study area (Figure 7e). As temperatures warm in the near future, and throughout the century, the probability of a 2nd generation was projected to increase, but remained low in Pinus habitats located above ~45°N on both continents (Figures 7f,g). As indicated above, 77% of NA areas with univoltine population growth (R > 1) in 1981–2010 were projected to be no longer thermally suitable for this strategy by the end of the century. Only 17% of the NA areas where univoltinism was lost were projected to be suitable for a successful 2nd generation (Figures 7d,h).

Discussion

Intra- and Inter-Continental Establishment Potential

Differences between I. typographus and D. ponderosae life history traits are reflected in the phenology models used for simulating thermal suitability for population success of each species. We modeled the role of positive density dependence in D. ponderosae tree attacks by restricting population success to those years when thermal conditions for development resulted in adult emergence that occurs at a seasonally appropriate time in summer and is synchronized among individuals to overwhelm defenses of live standing host trees (Powell and Bentz, 2009). Ips typographus attack success, in contrast, was assumed to depend on the availability of wind-felled trees with no or weak defenses (Christiansen and Bakke, 1988; Schroeder, 2010), and the I. typographus model determined whether populations undergo one or two generations between spring emergence and the following winter as a function of temperature (Jönsson et al., 2011). Moreover, two or more generations annually are common for I. typographus in warm European locations (Wermelinger, 2004), but not for D. ponderosae in warm NA locations (Bentz and Powell, 2014; Bentz et al., 2014). Our model simulations reflect and highlight these differences in the responses of both species to projected temperatures in their native and non-native habitats in NA and the Europe study area. Generally, I. typographus was projected to benefit from warming on both continents whereas thermal suitability for D. ponderosae was projected to decline, a result also found for D. ponderosae in NA based on correlative niche modeling (Sidder et al., 2016). In NA, 60% of the area thermally suitable for D. ponderosae in 1981–2010 was projected to be unsuitable by the end of the century, whereas very little I. typographus habitat was simulated to be lost over the same time period due to changing thermal conditions. The capacity of I. typographus to switch from univoltinism to a 2nd generation, and the requirement for D. ponderosae adult emergence synchrony, highlight how thermally dependent fitness traits can differ among species and also have positive and negative effects on population success in a changing climate.

Our results corroborate observed recent northward intra-continental range expansion of D. ponderosae in Pinus habitats of Canada (Cooke and Carroll, 2017) where at >60°N univoltine population growth was projected to increase 30% by mid-century as a result of warming in Pinus habitats in summer and winter. Thermal suitability for I. typographus univoltine populations was also projected to shift north in native habitats in the Europe study area where by mid-century the proportion years with univoltinism at 65°N will be higher than at any other time period. As the occurrence of the phenological strategy to produce a single generation annually shifts north for both species, our model simulations suggest response at lower latitudes will differ. A common strategy for I. typographus in Denmark and areas further south in Europe, except at the highest altitudes, has been to produce two or more generations in a year (Christiansen and Bakke, 1988; Wermelinger, 2004). Warming summer temperatures throughout this century are projected to result in an earlier completion of the first I. typographus generation and a northward expansion into Sweden and Finland of a successful 2nd generation. By the end of the century the majority (~70%) of Picea habitat in the Europe study area is projected to support a 2nd I. typographus generation.

Producing a 2nd successful generation is not a common strategy for D. ponderosae in its native NA habitat, and the probability of bivoltinism in the current D. ponderosae range and climate was projected to be low. Simulations suggest, however, that habitats further south of the current D. ponderosae range, including Mexico and the southeast US, were thermally suitable for D. ponderosae bivoltinism, highlighting the potential for intra-continental expansion southward. As temperatures warm, a northward expansion in thermal suitability for bivoltinism was projected for D. ponderosae in NA, although not to the extent as for I. typographus due to differing life history strategies that are reflected in the models. Because bivoltinism has historically been rare, conditions that will support bivoltinism in D. ponderosae, and alternatively conditions that might select against a 2nd generation, remain unclear. We model bivoltinism based on our knowledge of univoltine populations and the requirement for synchronous adult emergence. Two generations have historically been limited by cool fall and winter temperatures that can invoke a facultative prepupal diapause (Bentz and Powell, 2014; Bentz and Hansen, 2017). As winter temperatures warm, simulations suggest that in the warmest southerly habitats this constraint is relaxed, but potentially remains a barrier in northern latitudes. The capacity for predator/prey communities and symbiotic fungal associates to maintain synchrony with multiple D. ponderosae generations is also unknown. To fully understand the potential for D. ponderosae bivoltinism, and hence population success in a changing climate, additional research on community associates, host tree constraints for tree attacks at times other than the middle of summer, and limits of developmental plasticity are needed. As shown in our predictions, it is also likely there will be a lag period in thermal suitability of some Pinus habitats where neither the univoltine or bivoltine strategy are successful as habitats move through thermal transition zones (Roff, 1980). As the area supporting current adaptations for a univoltine D. ponderosae lifecycle is reduced and limits in the plasticity of thermally-regulated and behavioral traits are exceeded, adaptive changes in life history strategies will be needed in order to persist (Bale et al., 2002). Our models do not include the potential for adaptive capacity in life history strategies to changing temperatures, which when accounted for in modeling can significantly reduce projected range losses (Bush et al., 2016).

Ongoing range margin shifts northward are occurring in a wide range of taxa (Weed et al., 2013; Mason et al., 2015). Intra-continental range expansion of D. ponderosae northward and eastward in Canada is occurring through climate-assisted pathways (Cooke and Carroll, 2017). Ips typographus expansion southward and westward within Europe has been a result of extensive planting of P. abies beyond its natural range (Mayer et al., 2015). Inter-continental movement of non-native insects, however, involves human-assisted movement and with ever-increasing global trade the risk of invasion of new continents for both species will increase (Ramsfield et al., 2016). Wood-boring and ambrosia beetles have been a predominate guild intercepted at ports globally (Aukema et al., 2010), although introduction and establishment of non-native phloem feeders, including D. valens into China (Liu et al., 2014) and I. grandicollis in Australia (Neumann, 1987) have also been devastating. Many factors influence establishment and spread including host, habitat and environmental suitability (Liebhold and Tobin, 2008). Previous research indicates that suitability of breeding material should not be a barrier for either species, as I. typographus and D. ponderosae can both successfully reproduce in Picea and Pinus species, respectively, found outside their native ranges (Økland et al., 2011; Fries, 2017; Rosenberger et al., 2017; Flø et al., 2018). Our model projections suggest that thermal suitability for population growth will also not be a barrier to establishment for either species.

In the historical period (1981–2010), thermal suitability for D. ponderosae univoltine population growth was projected to be greater in non-native European Pinus habitat than in NA, a result of warmer summer and winter temperatures in the Europe study area than in comparable native habitats in NA. In the near future (2011–2040), thermal conditions in the Europe study area were projected to remain favorable for D. ponderosae univoltinism at high elevations and latitudes, and for a 2nd generation in lower latitudes. As temperatures warm throughout the century, however, thermal suitability for D. ponderosae is predicted to decrease except at the highest latitudes in Europe, with continued high probability of a 2nd generation at lower latitudes. During 1981–2010, I. typographus thermal suitability for univoltinism was also projected to be high across the majority of non-native Picea habitats in NA. In contrast to D. ponderosae, thermal suitability for one and two I. typographus generations annually is projected to increase throughout the century in NA. Although winter temperatures were much colder in NA than European Picea habitat during 1981–2010, warming winter temperatures throughout the century at the northernmost latitudes of NA will favor I. typographus overwintering survival.

In addition to habitat thermal suitability and availability of breeding material, the potential for inter-continental establishment and spread of both species could be hampered by life history strategies that influence population persistence at low levels including propagule size at invasion and dispersal capacity. Both I. typographus and D. ponderosae have sexual mating and population growth can be influenced by Allee effects. Allee effects cause a decline in per capita population growth when population density decreases below a critical threshold and can limit establishment when sufficient mates are not available (Liebhold and Tobin, 2008). Both species, however, have evolved strategies that allow for population persistence at low levels. Although at epidemic levels D. ponderosae favors large standing trees with active defenses, thereby requiring conspecifics to overwhelm defenses and create suitable breeding sites (Boone et al., 2011), very low levels of D. ponderosae can survive and reproduce in weakened trees that are often colonized by other beetle species and stem and root diseases (Bartos and Schmitz, 1998; Smith et al., 2011). Similarly, a common strategy for I. typographus is to attack wind-felled trees that have no or weak defenses, and reproductive success in these trees can be greater than in live standing trees (Komonen et al., 2011). Establishment potential for both species will therefore require forest stands with weakened or downed host trees. Although, strong Allee effects can reduce invasion speed (Lewis and Kareiva, 1993), the strength of the effect can vary across landscapes (Tobin et al., 2007) and may be overcome by the density of host plant material that is specifically suitable for low level populations (Powell et al., 2018), in addition to dispersal capacity (Chase et al., 2017). Passive or wind-assisted dispersal above forest canopies was responsible for D. ponderosae invasion into new habitats in Alberta, Canada (de la Giroday et al., 2012) and I. typographus is also known to have strong dispersal capacity (Forsse and Solbreck, 1985).

Model Limitations

We acknowledge several aspects related to phenology of I. typographus and D. ponderosae that could influence response to future climatic change and establishment potential in non-native habitats that are not currently included in our models, including intraspecific lineage differences in thermal traits. Induction of a facultative diapause in adult I. typographus and prepupal D. ponderosae varies with latitude in both species due to local adaptations to photoperiod (Schroeder and Dalin, 2017) and temperature (Bentz and Hansen, 2017). Diapause induction varies not only among populations but also within populations such that a propagule entering a new region could include a mixture of individuals with different responses to day length and temperature thereby increasing the probability that some individuals will be successful. Moreover, predicted phenological shifts can expose lifestages to novel day lengths and temperatures that will influence diapause induction and generation time. Intraspecific differences are currently not implemented in the D. ponderosae model, although geographic variability in day length regulation of autumn diapause is included for I. typographus.

The D. ponderosae model runs on an annual basis with univoltine adult emergence required to occur in mid-summer. The consequences and capacity for adaptation to emergence at other times is not currently included, nor is the contribution of sister broods to population growth. Thermal suitability for population growth of D. ponderosae is modeled based on what we know about successful univoltine populations historically, and requirements for the rare D. ponderosae bivoltine population success are unclear. Supraoptimal temperatures can reduce population growth not only through disruption of seasonality, which is included in the D. ponderosae model, but also through death (Tobin et al., 2014), a factor not included in our models.

We acknowledge that climate projections are uncertain, and that GCMs may overestimate mean annual temperature (Miao et al., 2014). In addition to the direct effects of temperature on population dynamics described by our models, drought can indirectly influence bark beetles through stress on host trees (Anderegg et al., 2015), and storm events through the creation of felled trees and easily acquired brood material (Kärvemo et al., 2014). These effects are not included in our projections, nor are interactions with fungal associates that could become desynchronized (Addison et al., 2013) or altered as novel associations form (Wingfield et al., 2016). Finally, although tree regeneration times are significantly longer than insect generations, tree distributions are continually changing (Aitken et al., 2008). We only considered current host tree ranges in predicting future distributions and voltinism patterns for D. ponderosae and I. typographus.

Conclusion

Forest management and carbon mitigation strategies in an era of climate change must encompass a range of uncertainties including periodic tree mortality events caused by native and invasive insects. As ectotherms, insect life history traits are thermally-dependent making their population success highly sensitive to climatic changes. Our process-based model projections for two important tree mortality agents in NA and Europe, D. ponderosae and I. typographus, highlight how future temperature changes across landscapes could result in positive and negative population consequences in part due to regional variation in climatic conditions and evolved adaptations for seasonal timing of overwintering lifestages. Population consequences could also differ between the species due to differences in the role of density dependence and traits that govern acquisition of breeding material. Simulations suggest that warming throughout the century will generally favor I. typographus in its native habitat in Europe and non-native Picea habitats in NA as thermal suitability for a 2nd generation (bivoltinism) increases. Long-term warming is not projected to be as favorable for D. ponderosae as thermal conditions for the historically optimal strategy of univoltinism declines. Bivoltinism will be a new required phenological pathway for D. ponderosae population success and it is imperative that adaptive potential and thermal, community, and host tree conditions supporting a 2nd generation are investigated.

Our models provide tools for evaluating temporal and spatial shifts in potential bark beetle-caused tree mortality in native and novel habitats, critical aspects of forest management in an era of climate change and increasing global trade. Our simulation results show extensive thermal suitability for D. ponderosae and I. typographus intra- and inter-continental establishment and population success in the near future (i.e., 2011–2040). Moreover, reproductive suitability of non-native hosts on each continent has previously been demonstrated (Økland et al., 2011; Fries, 2017). Collectively, these results highlight the importance of effective legislation to minimize risks of importing forest insects with potentially devastating effects on global conifer markets, and that continued monitoring should be remain a priority at ports on both continents.

Data Availability

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

BB, AJ, and MS conceived study. BB and AJ developed code for model analyses. RW provided climate data. KL and AW developed vegetation masks. All authors contributed to manuscript preparation, review, and revision.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This research was a collaborative effort as part of the International Union of Forest Research Organizations Task Force on Climate Change and Forest Health. We thank Jim Powell for helpful discussions on modeling and comments on a previous version of the manuscript. Guen Grosklos assisted with MATLAB coding. MS was financially supported by Bo Rydins Foundation (grant no. F 07/15), Swedish Research Council for Environment, Agricultural Sciences and Spatial Planning (Formas, grant no. 942–2016–11), BB by the Rocky Mountain Research Station, and AJ by Agricultural Sciences and Spatial Planning (Formas, and grant no. 2010–822). We thank editor Sigrid Netherer and the reviewers for helpful comments and suggestions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2019.00001/full#supplementary-material

References