Halyomorpha halys adults were field-collected (by sweep net, hand-picking from vegetation or in pheromone traps) between 2016 and 2019 from locations with recently-established populations in Georgia (including Abkhazia; 2016–2018; n = 293), Kazakhstan (2017; n = 11), Russia (2016, 2018, 2019; n = 202), Serbia (2018; n = 129) and Ukraine (2017; n = 11) (Fig. 1). Note that Abkhazia is a disputed territory within the Caucasus; however, for the purpose of this study, we considered Abkhazia as located within Georgia, based solely on the geographic continuity of the agricultural landscape in this area. As such, throughout this manuscript, samples from Georgia and the disputed territory of Abkhazia will collectively be referred to as samples from Georgia. Complete specimen and collection data are publicly available at www.boldsystems.org (Project Halyomorpha halys in eastern Europe and Eurasia, EEUR) and are summarized in Suppl. material 1: Table S1. Individual insects were stored in 95% ethanol for subsequent molecular analysis.

Figure 1.

Map of Halyomorpha halys collection locations in Serbia, Ukraine, Russia, Georgia and Kazakhstan.

As described by Gariepy et al. (2014), a single leg was carefully removed from each insect using flame-sterilised forceps and placed in an individual 200 µl well of a 96-well microplate, along with 2 µl of proteinase K (20 mg/ml) and 100 µl of 5% Chelex 100 Molecular Grade Resin (Bio-Rad Laboratories, Hercules, CA, USA). A negative extraction control containing the Chelex and Proteinase K solutions, but no insect tissue, was included in each microplate. Sealed microplates were incubated overnight at 55 °C, followed by 10 min at 99 °C. Samples were centrifuged at 5800 g for 5 min to pellet the Chelex solution and 50 µl of supernatant (containing DNA) was transferred to wells in a new plate, taking care not to transfer the Chelex residue along with the sample. Microplates containing the extracted DNA were stored at – 20 °C until further analysis.

PCRs were performed in a 25 µl volume containing 0.125 µl of Taq Platinum, 2.5 µl of 10× PCR buffer, 1.25 µl of 50 mM MgCl₂, 0.125 µl of 10 µM dNTPs (Invitrogen, Carlsbad, CA, USA), 0.25 µl of 10 µM forward and reverse primer (respectively), 19.5 µl ddH₂0 and 1 µl of template DNA. A 658-bp sequence of the mitochondrial gene Cytochrome C oxidase subunit I (COI) was amplified by PCR using primers LCO1490 and HCO2198 (Folmer et al. 1994). Thermocycling conditions included initial denaturation at 94 °C for 1 min, followed by five cycles of 94 °C for 30 s, annealing at 45 °C for 40 s, extension at 72 °C for 1 min, followed by another 35 cycles of 94 °C for 30 s, 51 °C for 40 s and 72 °C for 1 min and a final extension period of 5 min at 72 °C.

PCR products were visualised with a QIAxcel Advanced automated capillary electrophoresis system (Qiagen, Hilden, Germany) using the DNA screening cartridge and method AL320. Results were scored with QIAXCEL SCREENGEL Software (version 1.2.0) and only those samples of the expected fragment size with a signal strength exceeding 0.1 relative fluorescent units were scored as positive.

Samples, scored as positive, were purified using ExoSAP-IT (Affymetrix, Santa Clara, CA, USA), following the manufacturer’s instructions. Purified PCR products were bidirectionally sequenced on an ABI 3730 DNA Analyser at the Robarts Research Institute (London Regional Genomics Centre, ON, Canada). Forward and reverse sequences were assembled and edited using CODONCODE ALIGNER, version 9.0.1 (Codon-Code Corporation, Centreville, MA, USA). Sequence data and trace files were uploaded to the Barcode of Life Datasystems (BOLD; www.boldsystems.org) in the Project Halyomorpha halys in eastern Europe and Eurasia (EEUR).

Samples were grouped, based on their country of collection (Georgia, Kazakhstan, Russia, Serbia and Ukraine) and standard measures of diversity were calculated for each group using DnaSP v.5.10.01 (Librado and Rozas 2009), including number of haplotypes, haplotype diversity (h, the probability that two randomly-selected haplotypes are different; Nei 1987) and nucleotide diversity (π, the average number of nucleotide differences per site between two randomly-selected DNA sequences; Nei and Li 1979).

Samples were grouped, based on their country of collection (Serbia, Ukraine, Russia, Georgia and Kazakhstan) and the proportion of each haplotype within each group (i.e. country) was calculated in order to obtain a representation of the haplotype composition.

Additionally, based on current data and previous publications (e.g. Gariepy et al. 2014, 2015; Cesari et al. 2015, 2018; Šapina and Jelaska et al. 2018; Schuler et al. 2020; Yan et al. 2021), the number of COI haplotypes and the identity of dominant haplotypes from invaded European and central Asian countries were tallied and used to generate an overview of the trends and dominant haplotypes in these areas.

Genetic diversity measures for H. halys collected in Georgia, Kazakhstan, Serbia, Russia and Ukraine are shown in Table 2. For all samples combined, there were five haplotypes with a total of nine polymorphic sites. Overall haplotype and nucleotide diversity was 0.223 ± 0.021 (mean ± SD) and 0.00052 ± 0.00007 (mean ± SD), respectively. Haplotype and nucleotide diversity was zero in samples collected from Russia, Georgia and Ukraine, as only a single haplotype was observed. Serbia had the most diverse population observed in the present study, with five haplotypes recorded and showed the highest haplotype and nucleotide diversity (Table 2).

A total of 646 samples were analysed and a 658-bp fragment of the DNA barcoding region of the COI gene was generated (Genbank Accession numbers MZ871818 - MZ872463). Collectively, five COI haplotypes (H1, H3, H8, H33 and H80) were identified. The majority of the samples were identified as haplotype H1 (87.6%), followed by H3 (9.6%), H8 (2.1%), H33 (0.5%) and H80 (0.2%).

The proportion of haplotypes from each country is shown Fig. 2 and Table 3. Haplotype H1 was recorded in all five countries and was either the dominant haplotype or the only haplotype in four of the five countries: Georgia (100%), Kazakhstan (82%), Russia (100%) and Ukraine (100%) (Table 3). In Serbia, H1 was the second most common haplotype, but nonetheless, represented 39.5% of the haplotype composition, making it a major contributor to the haplotype composition. H3 was the dominant haplotype from Serbia, representing 46.5% of the haplotype composition and was also recorded from Kazakhstan (18% of the haplotype composition). Three additional haplotypes were recorded from Serbia: H8 (10.9%), H33 (2.3%) and H80 (0.8%) (Fig. 2 and Table 3). The known global distribution of these haplotypes in their native and invasive ranges is presented in Table 4.

Figure 2.

Map of Halyomorpha halys collection locations, with the COI haplotype frequency shown in pie charts sized proportionally to the sample size from each country.

An overview of the number of haplotypes and the dominant haplotypes in non-native countries in Eurasia is presented in Table 5. The haplotype composition in the majority of countries in this invaded region consists of a single haplotype or relatively few haplotypes (≤ 4) (Fig. 3). Only four countries demonstrate a moderate (5–7 haplotypes: Switzerland and Serbia) to high (> 8 haplotypes: Italy and Greece) number of reported haplotypes (Fig. 3). All countries with H. halys COI haplotype data in Eurasia are dominated or co-dominated by H1, with the exception of Switzerland and France, which are dominated by H3 (Table 5, Fig. 4).

Figure 3.

Trends in the number of cytochrome oxidase I (COI) haplotypes in the invasive range of Halyomorpha halys in Eurasia.

Figure 4.

Trends in the distribution of dominant cytochrome oxidase I (COI) haplotypes in the invasive range of Halyomorpha halys in Eurasia.

Invasive species typically have reduced genetic variation due to the occurrence of genetic bottlenecks upon colonisation of new locations (Fauvergue et al. 2012). Nonetheless, some species have become very successful colonisers despite strong founder effects (Sax and Brown 2000), particularly when multiple introductions from different locations contribute to enhanced diversity (Miller et al. 2005; Dlugosch and Parker 2008; Lawson Handley et al. 2011). Previous studies in Europe have suggested establishment of H. halys from multiple sources, including the direct establishment of Asian populations, as well secondary invasions via previously-established populations through the bridgehead effect (Gariepy et al. 2015; Valentin et al. 2017; Schuler et al. 2020). A low genetic diversity in invasive H. halys has been observed in European countries where only a single haplotype is present (e.g. Romania; Cesari et al. 2018; Yan et al. 2021) or where one haplotype is dominant amongst a small number of haplotypes (e.g. Switzerland, Hungary and France; Gariepy et al. 2015). In these countries, H. halys populations have established either directly from Asia (e.g. from China to Switzerland) or by secondary invasion from neighbouring European countries (e.g. from Switzerland to France), with haplotype diversity (h) ranging from 0 to 0.27 and nucleotide diversity (π) ranging from 0 to 0.0008 (Gariepy et al. 2014, 2015; Cesari et al. 2018; Yan et al. 2021). In contrast to other invaded countries in Europe, H. halys populations in Italy and Greece are substantially more diverse (h = 0.702–0.724; π = 0.0036–0.0054; Gariepy et al. 2015; Cesari et al. 2018), with establishments that have originated from multiple source locations, including directly from China, Japan and Korea and/or via the bridgehead effect from established populations in the USA (Cesari et al. 2015, 2018; Gariepy et al. 2015; Morrison et al. 2017; Valentin et al. 2017; Schuler et al. 2020).

In the present study, the overall haplotype (h = 0.223) and nucleotide diversity (π = 0.00052) of the COI barcode region was relatively low and is consistent with the values mentioned above from previous studies in most European countries (excluding Italy and Greece). However, diversity spanned a broad range, with no diversity in samples from Russia, Ukraine and Georgia (where a single haplotype was found) to haplotype and nucleotide diversity values of 0.62 and 0.00163 (respectively) from samples collected in Serbia (where a total of five haplotypes were found). Yan et al. (2021) found similar results in terms of a lack of haplotype diversity in their samples collected in Georgia, where a single COI haplotype was found. Our larger number of samples from several additional areas in Georgia provides a more thorough assessment of the populations in this region and confirms the observations by Yan et al. (2021) in that a single haplotype is (currently) present in this country. The population from Kazakhstan yielded diversity values similar to those reported in European countries with < 5 haplotypes (Fig. 3); however, it is based on a small sample size (n = 11) collected from a single site and, therefore, may not accurately represent the diversity in that country. Goodall-Copestake et al. (2012) recommend sample sizes ≥ 25 for accurate comparisons of population-level COI diversity; nonetheless, our limited samples from Kazakhstan provide a baseline dataset from a previously-unrepresented country that can be built upon with further sampling. This may also be the case with our samples from Ukraine, as the sample size was small and originated from a single region (Odessa). However, given the fact that current and previous studies in neighbouring countries (Romania, Hungary and Russia) have also shown little to no haplotype diversity, it is not surprising that similar results were found in Ukraine in the present study, despite the small sample size. Nonetheless, more thorough sampling in additional, geographically-diverse locations would help to confirm this observation and/or record changes in diversity over time as the invasion progresses.

The most diverse population in the present study was recorded from Serbia. Although Yan et al. (2021) found lower values for haplotype and nucleotide diversity from their Serbian samples, their values were based on nine specimens collected from a single site, which likely underestimated the actual diversity in Serbia (as per Goodall-Copestake et al. 2012). This is supported by the fact that Yan et al. (2021) recorded only two COI haplotypes, whereas the present study recorded five haplotypes and consisted of a much larger sample size from several locations (Fig. 1; Table 5). The diversity measures that we recorded from Serbia (h = 0.62; π = 0.00163) are just slightly lower than those reported in Italy and Greece (h = 0.702–0.724; π = 0.0036–0.0054; Gariepy et al. 2015; Cesari et al. 2018), where several haplotypes from multiple source populations (including China, Japan, Korea and USA) are known to occur (Gariepy et al. 2015; Valentin et al. 2017; Cesari et al. 2018). As will be discussed more thoroughly below, we suspect that H. halys populations in Serbia are derived from multiple source populations from more than one of the surrounding European countries, resulting in diversity levels more similar to those found in countries with multiple sources of invasion (e.g. Italy and Greece).

Based on the present study and the collective dataset available in literature for H. halys (see Table 5), the occurrence of COI haplotypes is fairly uniform across most of eastern Europe and central Asia (Figs 3 and 4). In most countries, only a single haplotype (H1) is present (e.g. Croatia, Romania, Ukraine, Turkey, Russia and Georgia) or is dominant alongside a minor contributing haplotype (typically H3; Kazakhstan and Hungary) (Figs 3 and 4). The exception in the present study is Serbia, where a total of five haplotypes were detected, including two dominant haplotypes (H1 and H3) and three additional haplotypes (H8, H33 and H80), two of which were minor contributors (< 5%; H33 and H80). A similar exception in previous studies was observed in Greece, where two dominant haplotypes are also known (H1 and H33; Fig. 4 and Table 5), along with several minor contributing haplotypes (H3, H13, H22, H30, H31, H32, H158, H159 and H160; Gariepy et al. 2015; Morrison et al. 2017).

Greece and Italy are known hotspots of invasive H. halys haplotypes, with 11 (Gariepy et al. 2015; Morrison et al. 2017) and 20 COI haplotypes (Morrison et al. 2017; Cesari et al. 2018; Schuler et al. 2020), respectively. Further, within the invaded range, many of these haplotypes are unique to these two countries (i.e. are not found elsewhere in Europe). For example, in the native and invasive ranges of H. halys, H80 was previously only known from Shandong Province in China (Zhu et al. 2016) and from northern Italy (Cesari et al. 2018; Schuler et al. 2020). However, in the present study we found this haplotype at very low levels in Serbia. Given the relatively low occurrence of this haplotype in Italy (Cesari et al. 2018), it is difficult to speculate whether the source of H80 is the result of movement and spread of H80 from Italy or whether it is a separate establishment originating from China. Similarly, H33 was previously only known from Greece (Gariepy et al. 2015; Morrison et al. 2017), but was recently detected for the first time in Italy by Schuler et al. (2020). Although H33 is known from China (Shanxi, Shaanxi and Anhui Provinces; Zhu et al. 2016; Valentin et al. 2017; Cesari et al. 2018), Schuler et al. (2020) suggest that the movement and spread of this haplotype from the already-established population in Greece is more likely the source of H33 in Italy, especially given the prevalence and persistence of this haplotype in Greece (Gariepy et al. 2015; Morrison et al. 2017). The detection of H33 in Serbia in the present study also suggests movement and spread of H33 through secondary invasion from Greece; however, it is unclear whether it is due to passive dispersal or associated with commercial trade and travel (Konjević 2020). As the occurrence of H. halys haplotypes in most of the countries that share a border with Serbia and Greece has not yet been investigated, it would be important to determine the occurrence of this (and other) haplotype(s) in such locations where H. halys is also known to occur (Table 1), in particular Bulgaria, North Macedonia, Albania, Montenegro and Bosnia and Herzegovina. This may provide insight as to whether the distribution of H33 is widespread or continuous across the Balkan countries or whether it is primarily in urban centres associated with commercial trade and travel. All of our samples were collected in north and north-eastern Serbia (Fig. 1), which is bordered by countries where a single haplotype (H1) is known or dominant (e.g. Croatia, Hungary and Romania; Fig. 4) and where H33 and H80 have not been reported. The movement of H33 and H80 into Serbia is, therefore, unlikely from this direction – neither through natural dispersal nor from commercial trade or travel. As such, sampling in the southern portion of Serbia would provide a more thorough account of the haplotype distribution across the entire country and provide insight on the movement and spread of H33 and H80, particularly from countries along the southern border of Serbia. In contrast, the movement of haplotype H1 likely occurred from the spread of this haplotype from neighbouring countries to the west (Italy), north (Hungary) and/or northeast (Romania), where H1 is dominant. Halyomorpha halys populations in Serbia were first observed in areas near or along the Serbian-Romanian border (Šeat 2015) and in close proximity to the railway line that connects Bucharest, Romania with Belgrade, Serbia (Musolin et al. 2018). Thus, the trapping and interception records support the secondary invasion of H1 from neighbouring eastern European countries (as opposed to separate establishment events from China), possibly associated with railway travel or movement of commodities on railway cars, as suggested by Musolin et al. (2018). In terms of H3, which was dominant in Serbia, the prevalence of this haplotype in western Europe (in particular Switzerland, France and Austria), as well as the known presence (albeit at low levels) in Hungary, could indicate movement and spread from this direction. Similarly, H8 is the second most common haplotype in Switzerland and is also present (at lower levels) in France and Italy (Gariepy et al. 2014, 2015; Cesari et al. 2018; Schuler et al. 2020), suggesting the direction of movement of H8 is likely from western Europe to Serbia. In neighbouring Croatia, only H1 has been reported by Šapina and Jelaska (2018), but this is only based on two specimens; further analysis of samples in Croatia would help clarify the occurrence and diversity of additional haplotypes, in particular H3 and H8, which would be interesting in terms of evaluating the spread of haplotypes from this direction. Serbia is located directly at the centre of the invasive range of H. halys, surrounded by countries with different haplotype compositions. Although we cannot exclude the possibility of separate invasion(s) from China, the likely scenario (based on location, haplotype and trapping data) is that H. halys is entering Serbia from more than one direction simultaneously (through natural dispersal, via commercial/horticultural trade and/or travel). Although largely speculative, H1 may have initially arrived in Serbia from neighbouring countries to the east (possibly with additional invasions from other directions, given that surrounding countries all have a high proportion of H1), with H3 and H8 arriving from the western European countries and H33 and H80 from Greece and Italy.

The first established population of H. halys in the area of eastern Europe and central Asia occurred in Sochi City, Russia in 2013–2014 and at the time the pest was not present in neighbouring countries (see Table 1). The fact that this establishment was geographically disconnected from the rest of the invaded range in Europe suggests that secondary invasion from natural dispersal of the pest is unlikely (Musolin et al. 2018). Although the source of H. halys could be the result of a separate introduction from Asia, a secondary invasion via the bridgehead effect is more likely, based on the timing and location of the arrival of H. halys and the events surrounding its establishment. Musolin et al. (2018) suggest that the pest was accidentally introduced from Italy or Greece with infested plant material that was used in massive landscaping efforts associated with the 2014 Winter Olympic Games hosted in Sochi; plants from northern Italy were regularly imported in 2012–2013 due to a similarity between the climates of the two regions and low availability of local stock to meet the landscaping demands leading up to the Olympic Games. This coincides with the occurrence of high population levels of H. halys in northern Italy and reports of economic damage to crops (Pansa et al. 2013; Maistrello et al. 2014). The H. halys population in Italy is predominantly H1 in most regions, in particular Emilia Romagna and Lombardy (Cesari et al. 2015; Morrison et al. 2017; Cesari et al. 2018; Schuler et al. 2020). This information, combined with the fact that the present study demonstrated that H1 is the only haplotype in Russia (based on > 200 specimens), lends support to the theory proposed by Musolin et al. (2018) that H. halys in Sochi may have originated from locations in Italy where H1 is dominant. From the focal point of Sochi City in Russia, H. halys may have dispersed to other regions in the Caucasus. However, populations of H1 in some of these locations (e.g. Sevastopol), as well as neighbouring locations in Ukraine (Odessa) and Georgia (Abkhazia, Adjara and Samegrelo), could be due to the spread of H. halys from other locations associated with the movement of commercial goods, as all of these cities are important seaports within the region. A similar concern in Australia was also flagged as a threat, when a significant number of live H. halys were intercepted in shipments arriving from Italy; haplotype H1 and H23 were both identified from shipments originating from Italy (Horwood et al. 2019), demonstrating how easily secondary invasion could occur, even over long distances.

In Kazakhstan, H. halys was first reported in 2016 in Almaty and establishment was confirmed when populations continued to expand in the area in 2017 and 2018 (Esenbekova 2017; Temreshev et al. 2018). Although the present study consisted of a limited number of samples from Almaty, Kazakhstan (n = 11), two haplotypes (H1 and H3) were detected, with H1 being dominant. Almaty shares a border with Xinjiang Province in China; however, H. halys is not known to occur in Xinjiang or Qinghai (Yu and Zhang 2007), indicating that the source of this pest in Kazakhstan is not from natural dispersal near or along the Kazakhstan-China border. Similarly, Almaty is far removed from the distribution of H. halys in the Caucasus region, suggesting that natural dispersal from this region is not responsible for the occurrence of the pest in Kazakhstan. However, Almaty is the major commercial centre of Kazakhstan and, as such, we speculate that commercial trade (with China, Russia and/or other European countries) would be the source of H1 and H3 in this country. Kazakhstan has a relatively low number of invasive alien species currently recorded; however, the recent rise in international trade and oil, gas and mining development in Kazakhstan, Turkmenistan, and Uzbekistan will also likely result in an increase of invasive species in these countries (Turbelin et al. 2017). Given Kazakhstan’s geographic location (i.e. directly between the invasive range of H. halys in Europe / Russia and the native range in China), its importance as the hub of international trade in central Asia and its position as a major transportation hub linking China to Russia and western Europe by air, rail, road and sea (Selmier 2020), we speculate that H. halys invasion from Europe, Russia and China are all very likely. Another invasive hemipteran, Leptoglossus occidentalis Heidemann (Hemiptera: Coreidae), was recently discovered in Kazakhstan (Barclay and Nikolaeva 2018) and is similarly far-removed from other known established populations. Its arrival is likely due to passive transportation of adults as stowaways in cargo or through nursery trade (Barclay and Nikolaeva 2018), which is likely the same pathway of entry for H. halys in this region. The collection and COI haplotype analysis of additional H. halys samples from this area would provide a more thorough documentation of the haplotype diversity in Kazakhstan. However, given its separation from the centre of other H. halys invasions in Eurasia and without interception records to corroborate potential pathways of entry, additional COI haplotype analysis will only tell us which countries have a similar haplotype composition and is unlikely to clarify whether the pathway of entry is directly from the native range in Asia or whether it is the result of a secondary invasion via Europe, Russia or some combination thereof (but see future directions below).

The present study focused solely on the COI gene, as this gene has shown reliability in terms of revealing geographic patterning (O’Loughlin et al. 2008; Valade et al. 2009) and has been widely utilised in haplotype studies on H. halys and has the most comprehensive, publicly-available global haplotype network available for this species. However, more in-depth multilocus analysis, based on microsatellite DNA or high-resolution genomic data (e.g. single nucleotide polymorphisms, SNPs; restriction site associated DNA sequencing, RADseq), may reveal additional patterns regarding invasion pathways (Garnas et al. 2016; Sunde et al. 2020). As techniques for generating high-resolution genomic data become more mainstream in ecological studies (see Andrews et al. 2016), future research that makes use of these newer techniques to investigate the global diversity of H. halys may provide a more fine-tuned interpretation of patterns of dispersal, in particular in terms of dissecting pathways of entry from the area of pest origin versus movement and spread of already-established populations to new locations.