Population genetics of the brown marmorated stink bug Halyomorpha halys in the early phase of invasion in South Tyrol (Northern Italy)

Abstract The brown marmorated stink bug Halyomorpha halys is one of the most harmful invasive species in the world. Native to East Asia, this insect was introduced into North America in the 1990s and into Europe in the 2000s where it subsequently established and spread across the continent. Previous population genetic studies determined the invasion pathways at continental and national levels. However, information on the dynamics on a small-scale is currently scarce. Here we study the genetic diversity and population dynamics of H. halys in South Tyrol, a region in Northern Italy, since its arrival to its widespread establishment over a period of four years. By haplotyping 162 individuals from ten populations (including six previously published individuals) we found a high haplotype diversity in most populations with an increasing diversity across the different years. Most haplotypes were previously found in other regions of Northern Italy, providing evidence for migration from neighboring regions. However, the presence of four previously undescribed haplotypes as well as a haplotype previously found exclusively in Greece highlights additional long-distance dispersal across the continent. Phylogenetic analysis of the haplotypes found in South Tyrol showed that the majority of haplotypes clustered with haplotypes predominantly found in Japan. This suggests a potential recent introduction of H. halys individuals from Japan into Europe, and thus an additional invasion pathway that was previously unidentified.


Introduction
The brown marmorated stink bug Halyomorpha halys (Stål, 1855) is an invasive insect pest introduced from Asia into several countries in the world (Cianferoni et al., 2018). Outside its native range it was first described in Pennsylvania, USA in 1996. It subsequently spread across the continent, invading most states in the USA before it was intercepted in Canada in 2010 Xu et al., 2014). In Europe, H. halys was introduced initially into Switzerland and Liechtenstein in 2004 (Arnold, 2009;Haye et al., 2014) and was subsequently recorded in neighboring countries including Germany, France, Italy and Greece in the following years (Heckmann, 2012;Callot and Brua, 2013;Maistrello et al., 2014;Milonas and Partsinevelos, 2014). H. halys has invaded Russia, Abkhazia, Georgia (Gapon, 2016) and Chile (Faundez and Rider, 2017), and was intercepted multiple times in New Zealand and Australia but no populations have established in these latter countries EPPO 2019).
Detailed knowledge of invasion routes is of considerable interest to prevent and control introductions and establishments of novel invasive species (Estoup and Guillemaud, 2010). In particular, the use of molecular and population genetic tools has shown to provide substantial information about the reconstruction of invasion routes. Several studies have described the invasion history and genetic diversity of H. halys using molecular tools with a particular focus on various mitochondrial regions (Cesari et al., 2014Gariepy et al., 2014;Xu et al., 2014;Zhu et al., 2016;Valentin et al., 2017;Lee et al., 2018;Kapantaidaki et al., 2019).
The comprehensive characterization of the genetic diversity of H. halys showed a high diversity within Asian populations and a significant genetic structure between populations from China, Korea and Japan (Xu et al., 2014). In contrast, invasive populations in the USA showed much lower genetic diversity than in the native range with only two haplotypes present before 2008 (Xu et al., 2014;Valentin et al., 2017) and three additional haplotypes described in 2016 (Valentin et al., 2017). In particular, the presence of a single haplotype (H1) in the east coast suggests that the invasion resulted from a single introduction event of a few founder individuals (Xu et al., 2014). Since in Asia H1 is exclusively present in China, it is assumed that the first individuals on the east coast were introduced from China (Valentin et al., 2017). In contrast, higher levels of haplotype diversity in the west coast indicate that multiple introduction events have occurred (Xu et al., 2014;Valentin et al., 2017). Europe, in contrast shows a different picture: while populations in Switzerland belong mainly to haplotype H3 and the less common haplotype H8, individuals in Northern Italy across different populations in the Emilia Romagna region are mainly H1 (Cesari et al., 2014;Gariepy et al., 2014Gariepy et al., , 2015. This suggests that while the first introduction of H. halys into Switzerland resulted from a direct interception from China, a second introduction occurred in Northern Italy, probably from an eastern US population (Valentin et al., 2017). Moreover, a high haplotype diversity was found in Greece. Since several haplotypesincluding the most common haplotype H33were outside of their native range only described in Greece, it is assumed that populations in Southeastern Europe resulted from an independent introduction event from China (Cesari et al., 2014;Gariepy et al., 2014Gariepy et al., , 2015Valentin et al., 2017).
A recent characterization of populations across Italy compared the haplotype frequency across different years and provided new insight into invasion dynamics on a smaller scale . In 2013 just two haplotypes were identified in Italy, H1 in the Emilia Romagna region and H3 in the Lombardy region (Cesari et al., 2014. While the population in Lombardy likely resulted from an introduction from Switzerland, H. halys in Emilia Romagna was probably introduced from Eastern USA (Valentin et al., 2017;Cesari et al., 2018). The number of haplotypes increased in just a few years with a total of 22 haplotypes recorded in 2016, likely as a result of multiple introduction events and range expansions from other European populations .
Although the invasion history of H. halys on a global (Xu et al., 2014;Valentin et al., 2017;Lee et al., 2018;Kapantaidaki et al., 2019), continental (Cesari et al., 2014;Gariepy et al., 2015) and national  scale has been well-investigated, studies describing the invasion dynamics on a smaller scale are currently scarce. Here we focus on the genetic diversity of H. halys in South Tyrol (Northern Italy). This region is of special interest as it is one of the largest contiguous areas of apple production in Europe and therefore threatened by this emerging agriculture pest. The first exemplars of H. halys were found in South Tyrol in 2016 on a cargo from Bergamo (Lombardy, Italy) (Unterthurner, 2016). Since 2016 H. halys was frequently found in South Tyrol suggesting a successful establishment in this region (Unterthurner et al., 2017). Moreover, from 2016 the frequency of H. halys detected in a monitoring program increased across the region, confirming the successful establishment and rapid spatial spread in this region (Fischnaller et al., unpublished). Whether the current population in South Tyrol resulted from a single introduction in 2016 and a subsequent spread or from multiple introduction events across different years is currently not known. We therefore aim to study the genetic diversity across different populations sampled across 4 years and compare our data to previously published global, continental and national data to understand the population dynamics of H. halys in this area.

Insect collection, DNA extraction and sequencing
A total of 156 individuals of H. halys collected from 2016 to 2019 across different localities in South Tyrol were analyzed. Localities within an area of 3 km were clustered into populations resulting in a total of nine different geographic populations (table 1). The number of individuals per population ranged on average from 14 to 16 individuals with the exceptions of Klausen (population 9, n = 18), Montan (population 9, n = 5) and Laimburg (population 7, n = 43). In total 21 individuals were sampled in 2016, 17 in Table 1. After taxonomic determination all individuals were stored in absolute ethanol at −20°C. DNA was extracted from a single hindleg using the GenElute Mammalian Genomic DNA miniprep kit (Sigma-Aldrich, Steinheim, Germany). DNA was eluted in 100 μl elution solution (10 mM Tris, 1 mM EDTA) and stored at 4°C. A 615 bp fragment of the mitochondrial COI was polymerase chain reaction (PCR) amplified using primers LCO1490 and HCO2198 (Folmer et al., 1994). This region has been characterized already in other studies (Gariepy et al., , 2015Xu et al., 2014;Zhu et al., 2016;Cesari et al., 2018;Lee et al., 2018) allowing a direct comparison with previously published data.
Each PCR reaction was performed in a total volume of 25 μl containing 7.5 μl H 2 O, 12.5 μl UCP HiFidelity PCR Master Mix (Qiagen, Hilden, Germany), 0.25 μM of each primer and 2.5 μl template DNA. All reactions were performed in an Eppendorf Mastercycler Gradient (Eppendorf; Hamburg, Germany) with the following conditions: 35 cycles with 30 s at 98°C, 10 s at 98°C , 10 s at 51°C, 15 s at 72°C followed by a final extension at 72°C for 2 min. All PCR products were run on a 1% agarose gel stained with SYBRsafe (Thermo Fisher Scientific) and Sangersequenced by Eurofins MWG Operon (Ebersberg, Germany).

Population genetic analysis
Sequences were edited manually and aligned with CodonCode Aligner v8.0 (CodonCode Dedham, MA, USA). Our sequences were compared to sequences available in Genbank using the BLAST algorithm in NCBI. To exclude PCR artifacts, sequences with ambiguous polymorphisms and individuals with unique haplotypes (see 'Results' section) were confirmed by sequencing the product of two independent PCR runs. New haplotypes were deposited in the Genbank under the accession numbers MT199091-MT199094. The number of haplotypes (H ), haplotype diversity (h) and nucleotide diversity ( p) were determined using the program DNAsp v6 (Librado and Rozas, 2009). Neighbor-joining tree analysis was performed with Mega X (Kumar et al., 2018) using the Tamura 3-parameter model (T92+G).

Results
We found a total of 15 different haplotypes across all populations of H. halys in South Tyrol. Eleven haplotypes have been described   . The other 12 haplotypes were found in low frequencies in one (0.61%; H155 and H167) to seven (4.29%; H78) individuals.
Overall, we found a high genetic diversity of H. halys in South Tyrol with a haplotype diversity of 0.684 ± 0.00128 and a nucleotide diversity of 0.00455 ± 0.00037(table 2). The number of haplotypes per population ranges from two in Montan (population 8) to eight in Klausen (population 9). Apart from the low haplotype and nucleotide diversity in Brixen in 2016  table 2), the lowest haplotype diversity in our study was found in Montan (population 8) with a haplotype diversity of h = 0.4 ± 0.237. The lowest nucleotide diversity was found in Laimburg (population 7) with p = 0.0024 ± 0.00065 (table 2). In contrast, the highest haplotype diversity was found in Meran (population 2) and Terlan (population 5) with h = 0.85 ± 0.057, while the highest nucleotide diversity was found in Meran with p = 0.00694 ± 0.0008 (table 2).
The comparison between individuals collected in different years shows an increasing number of haplotypes across the different years: while in 2016 six haplotypes and in 2017 four haplotypes were found, the number of haplotypes increased to nine in 2018 and 14 in 2019 ( fig. 2). This is reflected by a lower haplotype diversity (0.453 ± 0.116 and 0.419 ± 0.141) and nucleotide diversity (0.00321 ± 0.00101 and 0.00341 ± 0.00137) in 2016 and 2017, respectively, compared to 2018 (h = 0.716 ± 0.067 and p = 0.00368 ± 0.00073) and 2019 (h = 0.768 ± 0.039 and p = 0.00538 ± 0.00044; table 2).
The most widespread haplotype across all years and regions in South Tyrol is H1. This haplotype is the most common haplotype in the native range in Asia but also in most countries of the world. The haplotype is present in China and Japan but not in Korea, it is the most widespread in the USA, Canada and most European countries (Gariepy et al., 2015;Cesari et al., 2018) (table 3). In 2016 five other haplotypes were already present in South Tyrol. Haplotype H3 has been described in its native range in China (Gariepy et al., 2015) and is present in Hungary, Greece and Italy, and is described as the most prevalent haplotype in Switzerland and France (Gariepy et al., 2015;Cesari et al., 2018;Lee et al., 2018). In Italy this haplotype is frequent in Lombardy and has been recorded also in Piedmont, Emilia Romagna, Lazio and Calabria . Haplotype H23 was described in Switzerland, France  and in Italy where it was exclusively found in Veneto (Morrison et al., 2017) while H78 was described in China (Zhu et al., 2016) and is especially present in Liguria and southern Piedmont. In contrast, H153 has been solely found in Italy where it was exclusively recorded in Veneto table 2). One individual in 2016 belonged to a haplotype (H164) which has not been described yet.
Haplotypes H1, H3, H151 and H153 were still present in 2017 and all six haplotypes described in 2016 were found again in 2018. Three additional haplotypes were found in 2018: haplotype H8 has been described in low frequency in France and Switzerland as well as in Italy in low frequencies in Lombardy, Piedmont, Emilia-Romagna and Veneto  table 3). Haplotype H152 was exclusively found in Piedmont in low frequencies, whereas H41, a haplotype that originates in Japan, has been described outside its native range exclusively in Veneto   In 2019, besides H164 all haplotypes have been found again, whereas six haplotypes appeared for the first time. Haplotype H33 was described already in China and outside its native range exclusively in Greece (Gariepy et al., 2015). In 2019 this haplotype was found in Meran (population 2) and Klausen (population 9). Haplotype H180 occurs in low frequencies in Italy, in Lombardy, Piedmont, Veneto and has been found in Meran (population 2) and Bozen (population 6), whereas H155 outside its native range occurs only in low frequencies in Veneto, Trentino  and has been found in Klausen (population 9). Finally, three haplotypes (H165, H166 and H167) that have not been previously characterized in any other region in the world are scattered across South Tyrol.
Phylogeny and network analysis of the different COI haplotypes described in previous studies showed that haplotypes described in China and Korea formed a different clade than those found in Japan (Zhu et al., 2016;Valentin et al., 2017;Lee et al., 2018). The incorporation of the 15 haplotypes in South Tyrol revealed that six haplotypes (H1, H3, H8, H78, H80 and H152) clustered with haplotypes described in China,

Discussion
The invasion history of the brown marmorated stink bug H. halys in North America and Europe has been investigated in several recent studies (Gariepy et al., , 2015Valentin et al., 2017;Cesari et al., 2018;Lee et al., 2018). By comparing mitochondrial haplotypes of native populations from Asia with populations outside its native range, the invasion history was studied on a global Valentin et al., 2017;Lee et al., 2018), continental (Gariepy et al., 2015;Cesari et al., 2018) and national (Cesari et al., 2014Xu et al., 2014) scale. In the present study, we document for the first time the invasion history of H. halys on a regional scale in a time period of 4 years. We describe a remarkably high haplotype diversity across the whole range in South Tyrol with a high number of haplotypes that were increasing across the different years. The comparison of the haplotypes found in South Tyrol with previously published studies highlight a complex invasion history in this region, challenging the previously believed introduction hypothesis of H. halys in Europe. The most striking result of our study is the high number of haplotypes present in South Tyrol. General invasive populations express a lower genetic diversity than populations in their native range after experiencing a strong population bottleneck (e.g. Dlugosch and Parker, 2008). The founder effect was particularly strong in H. halys populations in North America where initially just one to two haplotypes were described in the USA Xu et al., 2014) and just three haplotypes in Canada Valentin et al., 2017). Although the number of haplotypes in North America increased likely because of additional introductions (Morrison et al., 2017;Valentin et al., 2017;Lee et al., 2018;Kapantaidaki et al., 2019), genetic diversity of H. halys in North America is still limited with a maximum of five haplotypes described (Valentin et al., 2017). Similarly, the number of haplotypes in various European countries is low, especially in Romania (one haplotype), Hungary (two haplotypes), France (three haplotypes) and Switzerland (four haplotypes) (Gariepy et al., 2015;Morrison et al., 2017;Valentin et al., 2017;Cesari et al., 2018;Lee et al., 2018). In contrast, populations in Greece show considerably Table 3. Overview of the haplotypes found in South Tyrol compared with descriptions in their native range (Xu et al., 2014, Valentin et al., 2017, invasive range in Europe and North America (Gariepy et al., , 2015Morrison et al., 2017, Kapantaidaki et al., 2019 and different regions in Italy
high genetic diversity with up to ten haplotypes described (Gariepy et al., 2015;Valentin et al., 2017;Lee et al., 2018;Kapantaidaki et al., 2019). Moreover, the use of multiple molecular markers revealed up to 22 haplotypes in Italy . On a regional scale up to ten different haplotypes were described in Piedmont and Lombardy based on a combined analysis of the mitochondrial COI and COII regions . Considering only the COI locus used in our study a maximum of eight haplotypes where described in Piedmont, seven in Lombardy and six in Veneto . Thus, with 15 haplotypes described in South Tyrol, this region represents the highest number and diversity of haplotypes outside the native range of H. halys.
H. halys was recorded in Italy for the first time in 2012, in the province of Reggio Emilia . The first individuals in South Tyrol were found four years later in March 2016, with some individuals detected in open fields and remarkably high numbers found in packing materials in shipments from Bergamo, Northern Italy (Unterthurner et al., 2017). By haplotyping 27 individuals (including six individuals studied by Cesari et al. (2018)) from 2016 we found six haplotypes present in South Tyrol: H1, H3, H8, H78, H153 and H164. By that time just two haplotypes, H1 and H3, were described in Bergamo, whereas H78 and H8 have been described in Brescia and in Como, respectively, approximately 50 km next to Bergamo. However, neither H153 nor H164 have been found in this geographic area . Although these findings are generally in concordance with a primary introduction of H. halys from Lombardy, the presence of haplotype H153 that has been described exclusively in Veneto (Morrison et al., 2017;Cesari et al., 2018) and H164, a haplotype that has not been described outside of South Tyrol indicates that additional introductions, likely secondary invasions from other areas have occurred.
Generally, we could not find a significant impact of geography on the haplotype composition, with H1 present in the majority of individuals from all populations and H3 and H153 that were present in most populations in minor frequencies. The other haplotypes were scattered across the region with no spatial structure. The haplotype diversity was generally higher in cities (Meran, Bozen and Klausen) than in villages with a particularly low nucleotide diversity in remote areas in Völlan, Laimburg and Montan. This emphasizes the important role of active and passive human-mediated transportation for the dispersal of H. halys but might also highlight a general preference for urban areas where a higher diversity of host plants and the presence of suitable overwintering shelter locations might benefit the insect (Hoebeke and Carter, 2003;Leskey et al., 2012;Wallner et al., 2014;Leskey and Nielsen, 2018).
Essentially all haplotypes described in different years were still present in the following years suggesting that introduced H. halys individuals were able to establish into the new territory. However, the number of haplotypes and the general haplotype diversity increased across the different years. While the number of haplotypes in 2016 and 2017 with six and four haplotypes was relatively low, it increased to nine and 14 haplotypes in the following two years. Although the higher sample size in 2019 might have affected this outcome, both haplotype and nucleotide difference increased across the years (table 2). This is in line with Cesari et al. (2018) who described an increasing haplotype diversity in Italy within a few years. We interpret this result as a consequence of repeated introduction of H. halys from different populations across Italy and other countries. Especially the detection of H33, a haplotype that outside its native range was found exclusively in Greece indicates that secondary invasions likely due to long-range human-mediated transportation might have occurred.
It is unclear, why the populations in a restricted area in South Tyrol show a higher haplotype diversity than other areas previously studied. Most of the studies describing the haplotype diversity of H. halys in invasive areas in Europe and North America describe populations sampled between 2011 and 2016 Morrison et al., 2017;Valentin et al., 2017;Cesari et al., 2018;Lee et al., 2018;Kapantaidaki et al., 2019) or earlier (Xu et al., 2014). It might be possible that additional introductions from the native range into Europe occurred recently and increased the haplotype diversity in other areas as well. Valentin et al. (2017) analyzed the invasion history of H. halys and concluded that most invasive populations in North America and Europe established from a direct introduction of H. halys populations from China with separate introductions into Eastern and Western USA and Canada, as well as into Switzerland and Greece. In contrast, the population in Northern Italy was established by a secondary introduction from Eastern USA into Emilia Romagna and the spread of the founder population in Switzerland to the neighboring countries. The authors found little evidence that populations outside China such as Japan and Korea established outside their native range (Valentin et al., 2017). Since half of the haplotypes found in South Tyrol have not been described in their native range, we performed a phylogenetic analysis to reconstruct their potential origin. Surprisingly, half of the haplotypes described in South Tyrol clustered with haplotypes that were previously described in Japan (Xu et al., 2014). This includes all four newly described haplotypes as well as four haplotypes that have been described previously in other regions in Italy . As limited commercial relations between South Tyrol and Japan exist, we consider a direct introduction of H. halys into South Tyrol unlikely. Nevertheless, our results suggest that an introduction of H. halys individuals from Japan might have occurred in Italy or neighboring countries recently, a hypothesis that needs to be tested in future studies.
Supplementary material. The supplementary material for this article can be found at https://doi.org/10.1017/S0007485320000553.