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 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.

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.

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).

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.

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).

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).

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).

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.

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).

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).

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).

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.

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.