Current distribution and voltinism of the brown marmorated stink bug, Halyomorpha halys, in Switzerland and its response to climate change using a high-resolution CLIMEX model

Climate change can alter the habitat suitability of invasive species and promote their establishment. The highly polyphagous brown marmorated stinkbug, Halyomorpha halys Stål (Hemiptera: Pentatomidae), is native to East Asia and invasive in Europe and North America, damaging a wide variety of fruit and vegetable crops. In Switzerland, crop damage and increasing populations have been observed since 2017 and related to increasing temperatures. We studied the climatic suitability, population growth, and the number of generations under present and future climate conditions for H. halys in Switzerland, using a modified version of the bioclimatic model package CLIMEX. To address the high topographic variability in Switzerland, model simulations were based on climate data of high spatial resolution (approx. 2 km), which significantly increased their explanatory power, and identified many more climatically suitable areas in comparison to previous models. The validation of the CLIMEX model using observational records collected in a citizen science initiative between 2004 and 2019 revealed that more than 15 years after its accidental introduction, H. halys has colonised nearly all bioclimatic suitable areas in Switzerland and there is limited potential for range expansion into new areas under present climate conditions. Simulations with climate change scenarios suggest an extensive range expansion into higher altitudes, an increase in generations per year, an earlier start of H. halys activity in spring and a prolonged period for nymphs to complete development in autumn. A permanent shift from one to two generations per year and the associated population growth of H. halys may result in increasing crop damages in Switzerland. These results highlight the need for monitoring the spread and population development in the north-western part of Switzerland and higher altitudes of the valleys of the south.


Introduction
Agricultural losses from insect pests are estimated up to 20-50% (Oerke 2006;Yudelman et al. 1998) and are projected to increase under future climate conditions (Battisti and Naylor 2009). Ectotherms are especially sensitive to an increase in temperature as their basic functions strongly depend on the environmental temperature (Deutsch et al. 2018). There is strong evidence that climate change is already modifying species geographical range, phenology, voltinism and interactions (Altermatt 2010;Bebber et al. 2013;Devictor et al. 2012;Hickling et al. 2006;Maclean and Wilson 2011;Parmesan and Yohe 2003;Robinet and Roques 2010;Roth et al. 2014;Walther et al. 2002). The type and strength of the effect of climate change on species depends on taxonomic group or location. The potential effect of future climate change on insect pests based on different emission levels and regional climate scenarios has been quantified in various studies (Caffarra et al. 2012;Harrington et al. 2007;Jonsson et al. 2009;Kistner 2017;Kocmankova et al. 2011;Langille et al. 2017;Meynard et al. 2013;Olfert et al. 2016;Porter et al. 1991;Stoeckli et al. 2012; Thomson et al. 2010;Trnka et al. 2007;Warren et al. 2018;Ziter et al. 2012). These projections indicate changes in distributions with expansions northward and uphill and contractions southward and downhill. Estimated changes in seasonal phenology include earlier spring occurrences, extended flight periods, and a change from uni-voltinism to multi-voltinism. Decreased winter mortality may also lead to higher population densities early in the year. Changes in species interactions (e.g. insect-host plant, host-parasitoid, competition or decoupling of mutualism) are realistic scenarios (Both et al. 2009;Gutierrez et al. 2010;Hance et al. 2007;Tylianakis et al. 2008). Climate change comprises various abiotic stressors, such as warmer temperatures, elevated CO 2 , more frequent droughts or storms leading to complex interactions with insect pests (Jactel et al. 2019).
Besides a shift in the distribution and phenology of native pest species, it is evident that the number of non-native species will increase and that climate change will promote their establishment (Bacon et al. 2013;Seebens et al. 2017). The brown marmorated stinkbug, Halyomorpha halys (Stål; Hemiptera: Pentatomidae) is native to East Asia (China, Japan, Korea and Taiwan) and is a highly polyphagous pest with more than 200 wild and cultivated host plants, mainly fruit trees in the family Rosacea, but also vegetables, berries, vine, maize or soya (Lee et al. 2013). Halyomorpha halys invaded North America in the mid-1990s (Hoebeke and Carter 2003), where it soon became a major pest in apples, peaches, sweet corn and beans in the Mid-Atlantic region (Leskey et al. 2012). The oldest photographic evidence for the presence of H. halys in Europe is from Zurich, Switzerland, in 2004(Haye et al. 2015, and since then it has established in many European countries (Cianferoni et al. 2018) and continues its spread throughout Eurasia (Tytar and Kozynenko 2020). Where the pest is bivoltine, severe damage has been observed in less than 5 years after establishment, e.g. Italy and the Republic of Georgia (Bosco et al. 2018;Maistrello et al. 2017).
Initially, the natural spread of H. halys in Switzerland was relatively slow, and populations were mostly restricted to the cities of Zurich and Basel. Until recently, H. halys has mainly been considered an urban nuisance pest. However, since 2015, significant damage in peach and pear orchards has been reported in the south, and since 2017, high population densities and damage were also noticed in fruit orchards in northeastern Switzerland. In northern Switzerland, H. halys is ostensibly univoltine (Haye et al. 2014). However, in recent years (2018-2019), two generations have been reported (Haye, personal observation), and it is assumed that two generations are only possible in warmer years. In the warmer climate of northern Italy, H. halys appears to be bivoltine .
Observations and projections of the impact of climate change on environment, industry and society are highly relevant to decisions and adaptation measures of governments, political and business sectors and the society (EEA 2017;Hofmann et al. 2011;Wilby et al. 2009). Regional impact initiatives, which are of high importance for developing adaptation measures, require localised understanding of how climate change affects natural and human systems (CH2014-Impacts 2014; Henne et al. 2018). Different modelling approaches have been used to estimate the response of native and non-native insects to climate change (Hill and Thomson 2015). One of them, CLIMEX, is a semi-mechanistic bioclimatic model package with proven reliability to explore the impact of climate on the potential distribution and seasonality of invasive species (Kriticos et al. 2015b). This model has some similarities with simple correlative models (Felber et al. 2018), but it simulates the mechanisms that limit species' geographical distributions and determines their seasonal phenology. CLIMEX models have been developed primarily for risk assessment of invasive pests (Byeon et al. 2018;de Villiers et al. 2016;Kistner-Thomas 2019;Kriticos et al. 2015c;Olfert et al. 2016;Ramos et al. 2019;Vera et al. 2002;Yonow et al. 2017).
CLIMEX models are often run on a spatial resolution of 0.16°(18.5 km) or 0.5°(55.5 km) climate data (Kriticos et al. 2012). However, this low spatial resolution does not allow adequate pest risk analyses for countries like Switzerland that have a complex topographic structure . Although H. halys is fairly widespread and abundant, most regions of Switzerland have been considered climatically unsuitable or only marginally suitable for H. halys in previous CLIMEX models run using the CliMond Climate Dataset (Kistner 2017;Kriticos et al. 2017). Accordingly, in the present study, a new CLIMEX version was developed to run with climate data of higher spatial resolution. Using this adapted CLIMEX model, the aim of this study was to reexamine the climatic suitability, population growth, the number of H. halys generations per year and the influence of different stress indices (e.g. cold, heat) under present and future climate conditions in Switzerland. To further validate the new fine-scale model, we were particularly interested in whether the model projections would actually match with long-term distribution records collected in a citizen science initiative in Switzerland. Finally, we aimed to identify highly suitable regions that have not yet been colonised by H. halys and those regions that may suffer from increasing H. halys populations in the near future in order to adjust current monitoring programmes and develop sustainable management strategies.

Bioclimatic modelling with CLIMEX
CLIMEX is a semi-mechanistic bioclimatic modelling package that can estimate potential, climatically suitable areas and population dynamics for poikilothermic organisms (Kriticos et al. 2015b). It combines weekly growth and stress functions to calculate the annual Ecoclimatic Index (EI), as an estimate of the annual climatic suitability of the location for a given species. EI values <1 indicate that a species cannot survive at the location. The weekly (GI W ) and annual (GI A ) Growth Indices are a function of temperature (TI) and soil moisture (MI) indices. Different stress indices as cold (CS), dry (DS), heat (HS) and wet (WS) stress and their interactions are considered to simulate the mechanisms that limit survival during unfavourable seasons. Furthermore, the minimum length of the growing season (PDD) and obligate diapause can constrain the overall climate suitability. The diapause Index (DI) is based on temperature and day length. The CLIMEX model for H. halys developed by Kriticos et al. (2017) was used to evaluate the potential distribution and population dynamics of this species in Switzerland (Table 1).

Spatial climate data and scenarios for Switzerland
Projecting present and future climate for regions characterised by complex topography is challenging, and substantial spatial variation is expected, especially for precipitation patterns (Jasper et al. 2004). The original H. halys CLIMEX model was calibrated with distribution records from Asia and validated with records from the invaded areas of the USA and Europe using a spatial resolution of 0.16°(18.5 km) . A new CLIMEX version (4.1) was recently developed in collaboration with CSIRO (Darren Kriticos, unpublished), to allow importing climate data of higher spatial resolution. Monthly average of daily minimum and maximum temperatures 2 m above ground level and the monthly precipitation sum interpolated from observed data at a regular latitude-longitude grid of 0.02°(2.2 km) were provided by MeteoSwiss. In total, we included 11,211 grid cells. Relative humidity at 9 am and 3 pm was deduced from the saturated vapour pressure at the corresponding air temperature and averaged at the monthly level (Kriticos et al. 2012).
CH2011 Extension Series provides gridded change signals of near-surface temperature at 2 m above ground and precipitation for the two non-intervention emission scenarios A1B, A2, and the climate stabilisation scenario RCP3PD (CH2011 2011; Zubler et al. 2014). These data cover the mean annual cycle for the periods 2020-2049, 2045-2074, and 2070-2099. Scenario data were created by adding the change signal to the monthly baseline data of the reference period . The scenario data were used to simulate the potential future distribution of H. halys. The same localised interpolated climate data set was used to analyse the climate suitability for H. halys at 10 Swiss locations. We selected the grid cell in which the locations falls based on longitude and latitude.

Potential distribution, phenology and voltinism under present and future climate conditions in Switzerland
The potential distribution of H. halys under present and future climate was evaluated based on the Ecoclimatic Index (EI).
We termed a location with an EI > 5 as climatically suitable and a location with an EI > 15 as climatically highly suitable. The shift in the weekly Growth Index (GI W ) was analysed to evaluate the impact of climate change on the phenology. The number of weeks with a GI W > 0 was calculated as a measure for how long H. halys is able to grow under present and future climate conditions. The number of H. halys generations per year was simulated for present and future climate conditions based on a value of 595 degree days above 12°C .

Validation of the CLIMEX model with observation records from Switzerland
The projected distribution of H. halys in Switzerland matched well with its observed distribution (Fig. 1). In the beginning, most observations were recorded from the cities of Zurich, Basel, Bern and Lugano, but in recent years (2017-2019), observation reports from other areas noticeably increased and confirmed the presence of H. halys throughout the Swiss midlands (Fig. 2a-e). For the same period, our model indicated a remarkable increase in the proportion of suitable areas (

Potential distribution, phenology and voltinism under present and future climate conditions in Switzerland
Compared to the reference period 1981-2010 (Fig. 1), the CLIMEX model projects an extensive expansion of the potential distribution in Switzerland for the A2 scenario (2070-2099) (Fig. 3). The north-western part of Switzerland could become completely suitable for H. halys. Southwards, the projected range expansion would reach the foothills of the Alps, and higher latitudes in the alpine valleys could become suitable under future climate conditions. The A2 model output indicates a potential 2.7-fold increase of the percentage of suitable areas (EI > 5) by 2099 based on the reference period (Table 2). Even for the mitigation scenario RCP3PD, a 1.6-fold increase was simulated. Specifically our simulations illustrate that approximately 40-60% of Switzerland will be suitable for H. halys long-term survival by the end of the century. The increase of suitable areas in Switzerland between 2020 and 2049 is similar for all climate scenarios (Table 2). However, between 2045-2074 and 2070-2099, the percentage of suitable areas under the RCP3PD mitigation scenario remains low compared to the A1B and A2 scenario. With the A2 scenario, the model projects a massive increase in the percentage of the area with a high EI (EI > 15): from 0.3% (reference) to 36.9% (A2) ( Table 2). There is no indication of cold, heat, wet or dry stress accrued under any of the climate scenarios (not presented). In Switzerland, the only factor leading to unfavourable conditions (EI = 0) is the relatively high number of degree days (595) required by H. halys to complete a generation.
For the reference period, we simulated that 33.2% of Switzerland would be suitable to support one or more generation of H. halys compared to 63.5% under the A2 climate change scenario (Table 2). Until only the last 2-3 years, two generations are only possible in exceptional years with unusual warm spring and summer temperatures (Table 2). By 2099, the A2 scenario projects two generations per year in 34.8% of the country compared to <1% for the reference period.    Fig. 1 for reference) A longer growing season under future climate scenarios is also visualised by the GI W (Fig. 4) for four different locations, chosen because they represent four typical climate regions in Switzerland. Spring activity of the bugs is projected to start 1-4 weeks earlier compared to the reference period, depending on the climate scenario and time period. Similarly, the fall activity would stop 1-3 weeks later compared to the reference period. Although overall heat and dry stress under future climate scenarios may not affect the suitability, during some weeks high summer temperatures may exceed the limiting upper temperature threshold (DV3 = 33°C). Consequently, for many locations in the southern and north-western parts of Switzerland (Fig. 4a, b), the projections indicate a decrease in the GI W during the hot summer months (weeks 25-35), especially for the time period 2070-2099. For other locations, mainly in the north-eastern part of Switzerland, the growth potential would increase from 2020 to 2099 (Fig. 4c, d).
The projections for 10 representative locations in Switzerland show the highly variable impact of climate change on the EI, number of generations and number of weeks with positive GI W (Table 3). According to the model, at nine out of ten locations the EI is expected to increase. Only in one location (Sion) a potential decrease is indicated due to high summer temperatures limiting the growth potential.
Remarkably, under future climate conditions, all localities have the potential to host two or more H. halys generations per year, with the number of weeks of growth potential (GI W > 0) projected to increase by 5-9 weeks (Table 3).

Discussion
The validation of the CLIMEX model using the observational records collected in a citizen science initiative revealed that more than 15 years after its accidental introduction, H. halys has colonised nearly all bioclimatically suitable areas in Switzerland. Accordingly, there is no strong potential for range expansion into new areas in Switzerland under present climate conditions. By using local Swiss climate data with a spatial resolution of 0.02°(approx. 2 km), we identified many more climatically suitable areas for H. halys in Switzerland in comparison to previous CLIMEX models with a lower resolution (Kistner 2017;Kriticos et al. 2017). As suggested by Kriticos and Leriche (2010), our study demonstrates the increase in explanatory power when projecting a species niche model to climate data of higher resolution. For future research, we recommend examining the goodness of fit of other mechanistic models that also consider factors such as seasonal abundance of H. halys, diversity of agricultural landscapes, or presence of preferred host plants, which are not reflected in CLIMEX models. In addition, most CLIMEX models use standardised (average) species values, but such a simplification ignores intraspecific (within-species) variation that can be as extreme as the trait variation across species.
Our results are in line with the projection that the global potential range of H. halys will expand polewards and contract from its (warmer) southern range limits under future climate scenarios (Kistner 2017). The estimated climate suitability of H. halys in Switzerland under future climate scenarios suggests an extensive range expansion into higher altitudes of mountain regions (Alpine valleys, Swiss Jura Mountains) that are currently too cold for establishment. The percentage of area suitable for H. halys (EI > 5) could more than double under the more extreme climate change scenario by 2099 (A2 , Table 2). Remarkably, we identified a pronounced increase in the proportion of highly suitable area (EI > 15) from 0.3% (reference) to 36.9% for the A2 climate scenario until the end of the century. More worryingly, climate change is projected to increase the number of weeks with positive GI W (Table 3 and Fig. 4), resulting in earlier H. halys activity in spring and a prolonged period for nymphs to complete development in autumn. This in turn increases the likelihood that in most years, H. halys will have two instead of one generation, as already observed in the unusual warm years 2018 and 2019. These observations were underlined by the simulated increase of the area with two generations from 0.3% (reference) to 34.8% (A2) until the end of the century. However, for some locations, 3-4 generations per year are projected. This may be an overestimation of the number of annual generations due to the parameters used for the current CLIMEX model (Tab. 1), which is based on 595 degree days above the critical temperature threshold (T 0 ) of 12°C per generation (PDD). However, literature values for PDD and T 0 are quite variable, ranging from 471 to 649 degree days (dd) and 11.1 to 13.9°C (Haye et al. 2014;Musolin et al. 2019;Nielsen et al. 2008). In addition, the preoviposition period of overwintering H. halys females (148 dd) (Nielsen et al. 2008) is not reflected in the PDD. Another factor that may have influenced the model output is that the accuracy of the local climate data is not the same for all Swiss regions, e.g. for the temperature, an imprecision of 0.6°C (midlands), 1.0 (alps) and 1.2 (Ticino) has been indicated (www.meteoschweiz.ch). Furthermore, we used the CH2011 climate scenarios based on 30-year mean temperature and precipitation changes, which does not account for potential changes in interannual variability, changes in wet-day frequency or dry-spell length.
Crop damage by H. halys in Switzerland was first reported in 2012 in greenhouse vegetables (Sauer 2012), and since 2016, damage in other crops such as pears, apples and grapes have steadily increased. Accordingly, a national monitoring programme was established in 2018 with the aim to better understand the distribution and seasonal phenology of H. halys in Switzerland. However, since H. halys was considered to still be spreading, the selection of monitoring locations was subject to uncertainty, and thus, the current model was developed to identify agricultural areas that are at higher risks to H. halys under current and future climates.
Our model projections suggest that the future increase in climate suitability and the expected expansion of the activity periods in spring and late autumn will likely result in a significant population growth of H. halys. Whereas warm temperatures in early spring would enable H. halys to develop two generations per year, unusual warm temperatures late in the season would allow more nymphs to complete development (Haye et al. 2014). Our model suggests that after 10 years of low population growth, such a scenario already occurred in the exceptional warm year of 2018. Consequently, a permanent shift from one to two generations per year and the associated population growth of H. halys would result in increasing crop damages in Switzerland, particularly in northern Switzerland and the Swiss midlands. Monitoring the spread and population development in the north-western part of Switzerland and higher altitudes of the valleys in the south, where new establishment of H. halys may occur in the near future, are recommended to protect threatened crops in a timely manner, using a variety of practices such as exclusion netting, bagging fruits, chemical or biological control (Candian et al. 2018;Lee et al. 2013;Leskey and Nielsen 2018).
The knowledge on the biology and ecology of H. halys is quickly growing and may help to further enhance the current model in the near future. Ideally, the precision of projecting  -20459, 2045-2074, 2070-2099) at the four Swiss locations a Changins, b Basel, c Wädenswil, d Güttingen species' distribution and seasonal phenology should not only consider high-resolution climate data, but also habitat factors. For example, Kriticos et al. (2015a) demonstrated that including such factors resulted in a 22-53% reduction in the areas estimated to be endangered by the invasive plant species, Parthenium hysterophorus (Asteracae). Increasing temperatures may also lead to an increase in crop water demand resulting in an expanding area requiring irrigation during summer months, especially in southern and central Europe (EEA 2017;Wada et al. 2013). For the invasive wasp, Vespula germanica (Hymenoptera: Vespidae), it was shown that including an irrigation scenario into CLIMEX significantly improved the model fit between present and estimated distribution (de Villiers et al. 2017). Similarly, the use of irrigation systems may reduce the potential dry stress of H. halys during summer months, particularly in southern Switzerland, where irrigation is common practise. This may explain the observations records of H. halys in locations, where our simulations indicate a season-long dry stress However, for a polyphagous insect like H. halys, a complex modelling effort would be needed to include local and crop specific irrigation levels with daily resolution, differentiated timing and different amounts of irrigation, especially for countries with a complex landscape like Switzerland. Haye et al. (2014) showed that overwintering mortality of H. halys adults in an open wooden shelter in Switzerland was relatively low (39 %). Although temperature refugia provided by human-built structures can be crucial for overwintering survival in some areas (Cira et al. 2016), we do not consider overwintering mortality a major limiting factor for the distribution of H. halys in Switzerland. Many CLIMEX models have focused on simulating the distribution of invasive pests, but only a few investigated the impact of climate change on the distribution of their antagonists (Haye et al. 2018;Olfert et al. 2016). Currently, there are no sufficient tools available for managing H. halys in Switzerland, but a very promising and sustainable approach could be the use of exotic natural enemies. The Asian samurai wasp, Trissolcus japonicus (Hymenoptera: Scelinoidae), has been identified as the most promising classical biological control agent in Asia and adventive populations were recently found in the Switzerland and Italy (Haye et al. 2020;Moraglio et al. 2020;Sabbatini Peverieri et al. 2018;Stahl et al. 2019). It remains unknown though, which areas will be most suitable for the establishment of the parasitoid in Switzerland and how climate change may alter its distribution. Similar to the current study, evaluating the bioclimatic suitability of T. japonicus, using the CLIMEX model of Avila and Charles (2018) and fine-scale climate data from Switzerland, could help to project the performance of the biological agent and identify optimal areas for its potential application. Table 3 Modelled climate suitability (CLIMEX Ecoclimatic Index, EI), number of generations and number of weeks with CLIMEX Growth Index, GI) > 0 for the reference period