Skip to main content
  • Research article
  • Open access
  • Published:

Genetic diversity analysis of brown marmorated stink bug, Halyomorpha halys based on mitochondrial COI and COII haplotypes



In the past decade, the brown marmorated stink bug (BMSB), Halyomorpha halys (Hemiptera: Pentatomidae) has caused extensive damage to global agriculture. As a high-risk pest for many countries, including New Zealand, it is important to explore its genetic diversity to enhance our knowledge and devise management strategies for BMSB populations. In this study, two mitochondrial genes, Cytochrome c oxidase I (COI) and Cytochrome c oxidase II (COII) were used to explore the genetic diversity among 463 BMSB individuals collected from 12 countries.


In total, 51 COI and 29 COII haplotypes of BMSB were found, which formed 59 combined haplotypes (5 reported and 54 novel). Of these, H1h1 was the predominant haplotype. The haplotype diversity (Hd) and nucleotide diversity (π) were high while the neutrality (Fu’s Fs) values were negative for the BMSB populations in the native countries, China, and Japan. For the BMSB populations from the invaded countries, the Fu’s Fs values were negative for populations from Chile, Georgia, Hungary, Italy, Romania, Turkey, and USA, indicating that those populations are under demographic expansion. In comparison, the Fu’s Fs values were positive for the populations from Austria, Serbia, and Slovenia, revealing a potential population bottleneck. Analysis of molecular variance (AMOVA) suggested that significant genetic difference exists among the BMSB populations from China, Japan, and the invasive countries.


This study revealed that the haplotype diversity of the BMSB populations was high in those two studied countries where BMSB is native to (China and Japan) but low in those countries which have been invaded by the species. The analysis indicated that multiple invasions of BMSB occurred in Europe and the USA. The study also revealed three ancestral lines and most of the novel haplotypes were evolved from them. Moreover, we observed two genetic clusters in the invasive populations that are formed during different invasion events. Our study provided a comprehensive overview on the global haplotypes distribution thus expanding the existing knowledge on BMSB genetic diversity that potentially could play an important role in formulating feasible pest management strategies.


Brown marmorated stink bug (BMSB), Halyomorpha halys (Stål, 1855) (Hemiptera: Pentatomidae) is a highly polyphagous insect, which feeds on more than 300 host plants [1, 2]. BMSB has been causing extensive damage to a wide variety of agricultural crops and poses a global economic threat for agricultural and horticultural industry [3, 4]. As per 2010 reports, the economic losses caused by BMSB valued more than 37 million USD in North America [5]. BMSB is capable of long-distance flight as well as local walking dispersal during the growing season [6, 7]. They are also considered nuisance pests as adults search for human-made structures to overwinter and discharge an unpleasant and long-lasting odor once disturbed [8].

BMSB is native to China (including Taiwan), Japan and the Korean peninsula [9,10,11] and has invaded 30 countries [12] including most of the states in the USA [4, 5], Canada [13], many countries in Europe [13,14,15,16,17,18,19,20] and Chile [21]. Climate modelling predicted that New Zealand at high risk for BMSB establishment should a successful invasion occur [22, 23]. The climatic suitability of this species is wide, ranging from 30 to 50 degrees of latitude, with an annual mean temperature range between 10 °C and 30 °C (NOAA-CIRES Climate Diagnostics Center, Boulder, Colorado). It covers wide geographical regions including northern Europe, north eastern part of North America, southern Australia, and the North Island of New Zealand [22, 23]. BMSB has increasingly been intercepted at the border and post-border scenarios due to the rise of international travel and trade [24] since 2005 when BMSB was first intercepted at the border of New Zealand [25]. BMSB is considered as “high risk” organism, and there were 1620 recorded interceptions of BMSB at the border of New Zealand since 2005 (accessed in June 2020) [25]. As agricultural exports comprise a significant proportion of the New Zealand Gross Domestic Product (GDP), the establishment of the pest would be highly detrimental to the country. The study by New Zealand Institute of Economic Research (NZIER) predicts a worst-case scenario of 3.6 billion New Zealand dollars in agricultural losses by 2038 if BMSB successfully establishes in New Zealand [26]. Therefore, it is imperative to study the genetic diversity of BMSB to trace the origin of captured individuals to better manage the border biosecurity risks and predict the pathways for trade.

Mitochondrial DNA (mtDNA) sequence analysis is one of the most widely used methods to examine genetic diversity and determine the locality of origin of an invasive species [27], and has been applied in a number of studies for BMSB [13, 20, 28,29,30]. Previous studies on BMSB populations from Asia and the USA have sequenced the mitochondrial Cytochrome Oxidase II (COII) and the ribosomal 12S genes to trace the origin of the invading BMSB populations in the USA [27]. Additional studies, focusing on BMSB populations in Europe, used Cytochrome Oxidase I (COI) gene alone [13], or in combination with the COII gene [20, 30], and revealed additional information on the genetic diversity of BMSB [13, 31]. To expand the knowledge on the genetic diversity of BMSB populations around the world, we collected BMSB adult specimens from different geographical locations, including the countries where BMSB is native to, China and Japan and the countries which have been invaded by the pest, such as the USA, Chile, and several European countries. In the current study, two mitochondrial barcoding genes, COI and COII, were sequenced from 463 BMSB specimens collected from 43 regions/provinces in 12 countries (2 native and 10 invaded) across four different continents to develop new insights into BMSB genetic diversity and their potential pathways of invasion.


COI and COII haplotypes in the BMSB populations

A total of 441 COI sequences (657 bp each) were obtained from 463 BMSB individuals collected from the 12 countries (Additional file 1). We identified 51 haplotypes using COI, consisting of 36 newly identified and 15 previously reported haplotypes (Additional file 2). For the sequences sharing 100% identity in the same region (657 bp) with those previously reported, the same haplotype names were given while new names were given accordingly for the rest of the sequences obtained in this study. All the new haplotypes identified were confirmed by BLASTn search [32]. The result showed that all the new haplotypes were unique in the COI region (657 bp). Further analysis showed that the haplotype N22 shared the same sequences with two shorter reference sequences, KY710432 (651 bp) and KY710450 (648 bp). However, it is not clear whether the missing bases from the two reference sequences are the same or different from the sequence we obtained, thus the sequence was considered as a new haplotype. The analysis also indicated that it might not be accurate to assign the same haplotype name for the sequence with different length.

There are three Barcode Index Numbers (BINs) for the COI sequences of BMSB specimens submitted in Barcode of Life Data System (BOLD) [33]: BIN AAM9563 (containing 682 sequences), ADT6053 (one sequence) and AAK5312 (six sequences). Further comparison with those COI sequences in BOLD showed that all the sequences obtained in this study belong to BIN AAM9563, with over 98% identity. In contrast, the sequences shared over 94 and 82% identity to the sequences of two other BINs, ADT6053 and AAK5312, respectively.

A total of 450 COII sequences (518 bp) were obtained from the 463 BMSB individuals (Additional file 1), which formed 29 haplotypes including 20 novel and 9 previously reported (Additional file 2). BLASTn search [32] showed that the new haplotypes did not share any identical sequences in the 518 bp overlap COII region with that of the previous reported sequences.

The geographical distribution of the identified COI and COII haplotypes are shown in Figs. 1 and 2. Of the identified haplotypes, H1 (61.9% of the total individuals) and h1 (61.7% of the total individuals) were predominant for COI and COII, respectively, and were detected in all the countries studied except Japan (Table 1 and Additional file 2). Haplotypes H3 (7% of the total individuals) and h3 (16% of the total individuals) were the second most predominant haplotypes detected in China, Austria, Chile, Hungary, Italy, Serbia, and Slovenia. In addition, haplotypes H8 and H48 for COI were only detected in Austria. The newly identified haplotypes were mainly observed in the native countries (Fig. 1a for China and Fig. 1b for Japan) except N47 in Slovenia (Additional file 2). All the novel COI and COII haplotypes identified were detected from the two native countries, China (Figs. 1 and 2a) and Japan (Figs. 1 and 2b).

Fig. 1
figure 1

Geographical distribution of mtDNA COI haplotypes (657 bp fragment of the COI) of BMSB in native countries - (a) China, and (b) Japan. Provinces where BMSB was collected from are shaded in grey. The size of pie is proportional to the frequency of haplotypes. Each colour represents a different haplotype. The figure was generated using Tableau 2019 ( [34]

Fig. 2
figure 2

Geographical distribution of mtDNA COII haplotypes (518 bp fragment of the COII gene) of BMSB in native countries in the current study – (a) China, and (b) Japan. Regions where BMSB was collected from are shaded in grey. The size of pie is proportional to the frequency of haplotypes. Each colour represents a different haplotype. The figure was generated using Tableau 2019 ( [34]

Table 1 Percentage (%) of the dominant mtDNA COI and COII haplotypes of H. halys detected in the studied countries. The COI haplotypes are named in uppercase letters. The haplotype names with a prefix ‘H’ represent the previously reported haplotypes while those with ‘N’ are the newly detected haplotypes identified in this study. The percentages of the individuals for each dominant haplotype in the country are listed in the table

Overall, high haplotype diversity was observed in China. The main haplotypes from China were H1, H33, H22, H3 for COI and h3 and h1 for COII (Table 1). The predominant haplotypes from Japan were H23, H45, N22, N40 for COI and h11 for COII (Table 1). Outside of the native regions, low haplotype diversity was observed, and H1, H3 for COI, h1, h3 for COII were the main haplotypes detected in those countries. Only one haplotype of each (H1 and h1) was detected in Georgia, Romania, Turkey and the USA (Table 1).

COI-COII combined haplotypes of the BMSB populations

In total, 428 individuals were identified with both COI and COII sequences (Additional file 1), and thus used for COI-COII combined haplotype analysis. The combined COI-COII haplotype analysis produced 59 haplotypes, in which only five were previously reported and 54 were novel (Additional file 2). All these newly identified haplotypes were detected in China and Japan except a single haplotype in Slovenia (N47h3). The predominant haplotype H1h1 (62.6%) was observed in all the countries except Japan (Additional file 2). The geographical distribution of the identified COI-COII combined haplotypes is shown in Fig. 3. In the native countries of BMSB, high haplotype diversity was observed with 24 haplotypes in China (Fig. 3a) and 32 in Japan (Fig. 3b), with no haplotypes shared between the two countries (Additional file 2 and Fig. 3). In comparison, out of the 32 haplotypes identified in Japan, 31 were uniquely detected in Japan, and one haplotype, H41h15 was shared with an individual from Hungary (Additional file 2). Similarly, 22 out of 24 haplotypes detected in China were unique, and two haplotypes (H1h1 and H3h3) were also predominantly shared with the BMSB samples from those invaded countries (Additional file 2). In the invaded countries, H1h1 was the predominant haplotype, identified in more than 90% of the studied samples from most of the BMSB-invaded countries, including Chile, Georgia, Hungary, Italy, Romania, Turkey and the USA (Additional file 2).

Fig. 3
figure 3

Geographical distribution of mtDNA COICOII haplotypes (1175 bp fragment) of BMSB in native countries in the current study, (a) China, and (b) Japan. Provinces where BMSB was collected from are shaded in grey. COI-COII combined haplotypes were formed by combing COI and COII haplotypes. The size of pie is proportional to the frequency of haplotypes. Each colour represents a different haplotype. The figure was generated using Tableau 2019 ( [34]

Population genetic analysis based on the combined haplotypes of COI and COII

Japan and China had the highest haplotype diversity (Hd), with Hd values of 0.942 and 0.858, and nucleotide diversity (π) values of 0.00238 and 0.00327, respectively (Table 2). Outside of the native regions of BMSB, the highest haplotype diversity was observed in Austria (Hd = 0.686, π = 0.00206), Serbia (Hd = 0.556, π = 0.00095) and Slovenia (Hd = 0.514, π = 0.00115). In contrast, little to no haplotype diversity was observed in the BMSB samples collected from Chile, Georgia, Hungary, Italy, Romania, Turkey and the USA. Therefore, two genetic groups were defined based on the Hd values obtained from the haplotype analysis: group A (Chile, Georgia, Hungary, Italy, Romania, Turkey and the USA) and group B (Austria, Serbia and Slovenia). It is noteworthy that in Hungary five sampling sites were studied, of which at two sites no haplotype diversity was observed, while the other three sites showed variable diversity with an Hd value from 0.038 to 0.5 and a π value from 0.00085 to 0.0017036, with an overall Hd value of 0.107 and a π value of 0.00028. This indicates that the invasion of BMSB in Hungary may have come from genetically distinct populations.

Table 2 The sample information and mtDNA diversity. The total sample size (N) for each province and country is listed. The number of haplotypes (Hn), haplotype diversity (Hd) and nucleotide diversity (π) were calculated based on the combined haplotypes of COI and COII. Two genetic groups were identified based on the COI and COII diversity from invaded populations. The genetic group A comprises Chile, Georgia, Hungary, Italy, Romania, Turkey, and the USA. The genetic group B comprises Austria, Serbia, and Slovenia. The analysis was conducted for each population, the populations from each country and each genetic group. Asterisk (*) represents statistically significant difference (p < 0.02)

In the neutrality test, the Fu’s Fs statistic values were very low in the two native countries of BMSB, China and Japan, with − 7.852 (p < 0.02) and − 29.707 (p < 0.02) (Table 2) while for the BMSB-invaded countries, Fu’s Fs statistic value was − 1.174 (p < 0.02) for group A (Chile, Georgia, Hungary, Italy, Romania, Turkey and the USA) suggesting that group A was under population expansion. In comparison, the haplotype diversity was slightly higher with an average of 0.63 for group B (Austria, Serbia, and Slovenia), but with Fu’s Fs values of 1.453 (p > 10), indicating that group B was under bottleneck.

The Principle Coordinates analysis (PCoA) using the FST values showed that there were at least three population clusters, namely China, Japan and group A (Chile, Georgia, Hungary, Italy, Romania, Turkey and the USA) (Fig. 4, Additional file 3). The recent invasion in Slovenia showed genetic similarities to those from Hebei and Beijing provinces of China (Fig. 4). The BMSB populations from Austria and Serbia were also closely related to the Chinese populations of Shanxi and Anhui. The population from the Chinese province of Hainan also showed close relationship with a population from the Japanese province of Akita.

Fig. 4
figure 4

Principal Coordinates Analysis (PCoA) plot based on population pairwise genetic distances. The points from Austria, Chile, Georgia, Hungary, Italy, Romania, Serbia, Slovenia, Turkey and the USA represent the samples collected from one country while the points from China and Japan represent the samples collected from one province. The colour represents the countries where the samples were collected from. The provinces are labelled with name tags. The X axis is the value of Coord.1 while Y axis is the value of Coord.2. Percentage of variation explained by coord.1, Coord.2 was 43.39 and 12.99%, respectively. The figure was generated using Tableau 2019 ( [34]

The AMOVA (Analysis of molecular variance) showed that the genetic variation among the 12 populations contributed 71.26% while variation within population contributed 28.74%. The overall FST value was 0.71 (p < 0.05), indicating that the genetic variation among populations was high.

The haplotype network of the BMSB individuals further revealed the widespread occurrence of H1h1 and H3h3, except the populations from Japan, whereas all the other haplotypes were mainly detected in the native countries (Fig. 5). The analyses showed that there were three ancestral lines found in this study namely h1, h3 and h11. Most of the other haplotypes mutated from these three lines with differences of several base pairs. Moreover, an interesting phenomenon was observed that some haplotypes (N3n3, N5n3, N4n4, N5n5) detected only in the Hainan population (China) was highly isolated and closer to Japanese populations rather than to Chinese populations. To further explore the distribution of those haplotypes, the combined COI and COII dataset from the present and previous studies [20, 30] were analysed together and resulted in a total of 80 haplotypes. The haplotype network analysis (Fig. 6) indicated similar genetic relationships as previously reported except that a few BMSB specimens from Italy had close relationship with Japanese populations (Fig. 6).

Fig. 5
figure 5

Haplotype network derived from the TCS analysis using COI-COII combined haplotypes (1175 bp fragment) of BMSB around the world. Each pie represents one haplotype with the haplotype name next to it. The size of pie is proportional to the frequency of haplotypes. The colour represents the countries where the samples were collected from. The tick marker on each line represents a base pair difference. The figure was constructed based on the combined COI and COII haplotypes obtained in the current study

Fig. 6
figure 6

Haplotype network derived from the TCS analysis using COI-COII combined haplotypes (1175 bp fragment) of BMSB around the world. Each circle represents one kind of haplotype. The size of the circle represents the frequency of each haplotype. The colour represents the countries where the samples were collected. The tick marker on each line represents a base pair difference. The figure was constructed based on the haplotypes obtained in the current study and those from the study of Cesari et al. [30]. The haplotypes obtained from Cesari et al. [30] are labelled with TH


This study revealed 51 COI haplotypes (36 novel) and 29 COII haplotypes (20 novel) from 463 BMSB individuals of 12 countries. However, most of these haplotypes (80%) were detected only once (Table 1 and Additional file 2), indicating that these new haplotypes are less abundant in the populations we studied. The haplotype analysis of mtDNA sequences of the BMSB populations from 12 countries provided genetic information for the identification of the pathways of invasion and the possible sources of origin.

In terms of haplotype distribution, the predominant COII haplotypes for Beijing (China), Shaanxi (China), Japan and the USA were h1 (67.8%), h3 (75%), h11 (81.1%) and h1 (100%). Xu et al. [29] made a similar conclusion that the major COII haplotypes for Beijing (China) and Shaanxi (China), and Japan and the USA were h1 (50%), h3 (100%), h11 (38%) and h1 (100%). Lee et al. [28] identified COI haplotype H1 as the main COI haplotype for China (68%), Hungary (98.8%), Italy (80%) and the USA (92.5%). Similarly, this current study identified haplotype H1 as the predominant one, accounting for 27.6% of the Chinese samples, 93.4% of the Hungarian samples, 95.8% of the Italian samples, and all samples from the USA. The percentage of each haplotype in those countries varies slightly between the current and previous studies, which is likely due to the differences in the sample size. On the other hand, some divergences were also observed between the current and previous studies. For example, COII haplotype h14 was identified as the predominant haplotype (33%) in the Japanese population in Xu’s study [29] but was not detected in this study (Additional file 2). Variations in the haplotype numbers and percentages between the current study and that of Cesari et al. [30] were also observed. Cesari et al. [30] reported a total of 26 COI and COII combined haplotypes for BMSB specimens, mainly from Italy [30]. In contrast, only five haplotypes [H1h1 (TH1), H3h3 (TH4), H33h3 (TH25), H8h11 (TH11) and H45h11 (TH22)] were detected in our study, thus, the total number of known BMSB COI and COII combined haplotypes (both known and novel) has increased to 80 known so far. Cesari et al. [30] included relatively large number of samples from Italy (209 samples from 10 regions) and identified 22 unique COI and COII combined haplotypes. Our study was unable to identify any additional COI and COII combined haplotypes beside the common ones (H1h1, H3h3) in Italy due to the small sample size (24 samples from 3 locations) studied. However, there were five shared haplotypes [H1h1 (TH1), H3h3 (TH4), H33h3 (TH25), H8h11 (TH11) and H45h11 (TH22)] among the samples studied here and by Cesari et al. [30]. Of these, the first two haplotypes were the two most common haplotypes. Interestingly, the last 3 haplotypes [H33h3 (TH25), H8h11 (TH11) and H45h11 (TH22)] detected in Italy by Cesari et al. [30] were also detected in this study from China, Austria, and Japan. After combining the haplotype data from the two studies, haplotype network (Fig. 6) revealed that the BMSB populations in Italy had genetic relatedness to Japanese populations, sharing the same haplotype H45h11 (TH22). It has been shown that the predominant haplotypes, such as H1h1 and H3h3, found in Italy were also widespread in China, the USA, and other European countries [20, 30]. Therefore, it can be hypothesised that the BMSB populations in Italy possibly have originated from Asia, which can be supported by the extensive, ongoing cross-border travel and trade between Asia and Italy. However, invasion from North America cannot be ruled out as H1h1 was found also in the USA. Furthermore, haplotype H8 was detected in Switzerland [13] in 2012, in France [31] and Northern Italy in 2013 [20], suggesting the possible invasion of BMSB to Italy was from Switzerland [30] based on the geographical proximity. The detection of the same haplotype of H8 in Austria in the current study raised the possibility that the invasion in Austria might have originated from the neighbouring countries of Switzerland or Italy as they share borders. However, the widespread distribution of H1 and H3 haplotypes in Austria (Table 1) opens the likelihood of invasion from China as well.

The combined data for COI and COII led to the observation of several ancestral lineages, including h1 (H1h1), h3 (H3h3, H22h3, H33h3) and h11 lines (H45h11, N22h11, N40 h11). The haplotype networks support that most of the less abundant haplotypes were possibly evolved from these lines.

The haplotype diversity of BMSB from the native regions is much higher than that of the invaded populations. The haplotype diversity (Hd) of BMSB populations from China and Japan was 0.858 and 0.942, respectively, which clearly indicates that the genetic diversity of these two populations was much higher than that of the most of the BMSB populations in the invaded countries (Hd < 0.182). This conclusion agreed with the observations from the study by Xu et al. [29]. Another interesting result from our study was the absence of shared haplotypes from the two neighbouring native countries of BMSB, China and Japan, which was congruent with previous studies [20, 29]. The haplotype network highlighted that China and Japan had their own haplotype clusters, suggesting that there is limited or no interbreeding, probably the geographical barriers of the strait restricted gene flow among populations. It is also possible that China and Japan have put in place a strict quarantine inspection for BMSB to prevent human-mediated transportation with current intensive trade activities. As a result, these BMSB populations could have been evolving independently. The Principal Coordinates Analysis also supported this conclusion, in which, most of the Japanese and Chinese populations were clustered by themselves except populations from regions of Akita (JP) and Hainan (CN). This phenomenon was also found in Zhu’s study [35] that the BMSB from Hainan and Japan were in the same clade.

The PCoA also revealed that the genetic group A (Chile, Georgia, Hungary, Italy, Romania, Turkey, and the USA) could have become a relatively independent genetic group. The neutrality test supported this, where the Fu’s Fs value of the genetic group A (− 1.174, p < 0.02) was negative as that of the ancestral Chinese (− 7.852, p < 0.02) and Japanese populations (− 29.707, p < 0.02), indicating that the populations belong to genetic group A could have been under population expansion stage [35]. The relatively low differences in FST values between the populations within the genetic group A indicate that these populations could have originated from the same ancestral line. The close genetic relatedness among the BMSB populations from Chile, Georgia, Hungary, Italy, Romania and Turkey with the USA populations suggests that the late detection in those European countries might have originated from the USA. This also aligns well with the past invasion history of BMSB, where BMSB was first detected in the USA in mid-90s, and then spread throughout the country. Ten years after establishing in the USA, BMSB was detected in Switzerland in 2007 [36], then spread in Europe [14,15,16,17,18,19, 37] and Chile [21] recently. This study further suggested the secondary invasion from the USA to the European countries such as Georgia, Hungary, Italy, Romania and Turkey, and to Chile. In contrast, the populations of the genetic group B (Austria, Serbia and Slovenia) were clustered with the Chinese populations, but were genetically distant from the group A, indicating that these populations in group B originated from a different pedigree line or genetic group. The neutrality test’s Fu’s Fs value of genetic group B was positive (1.453, p > 10), which is consistent with the more recent detection of BMSB in these locations as the positive Fs indicated that these populations were under a population bottleneck [38]. The haplotype diversity (Hd) and nucleotide diversity (π) of these three populations in group B were higher than 0.5 and 0.001, respectively, implying that these countries have been invaded multiple times by BMSB from different origins.

This study further revealed that the predominantly common haplotypes H1h1 and H3h3 exist in China and the invaded countries. The reasons for this need further investigation, though one possibility could be due to the dominant distribution of these haplotypes in Asia, and thus they have higher chances to be transported passively to other regions. It is not clear whether there is a possibility that these BMSB haplotypes can adapt more easily to new environments than other haplotypes. Since the COI and COII haplotype analyses are only based on the information from the female lineages, further study on the genomic level using high-throughput sequencing techniques might be able to provide more information on the genetic diversity of the BMSB populations and may help clarify the past invasion scenarios in the future.


The present study has revealed genetic diversity among BMSB populations using combined COI and COII datasets and provided better understanding of their potential invasion pathways. The genetic diversity among the BMSB populations from the native regions was much higher than those from the BMSB-invaded countries. The haplotype analysis further indicated that the invasion of BMSB has occurred multiple times in the past, probably at least partially due to international trade and travel. BMSB populations from the invaded countries, such as Chile, Georgia, Hungary, Italy, Romania, Turkey, and the USA were genetically close, but well separated from the Chinese populations. However, the BMSB populations from Austria, Serbia and Slovenia were more closely related to the Chinese populations. The results indicated that some individuals of the recent invasions into Chile, Hungary, Georgia, Turkey, Romania and Italy potentially originated from the USA without ignoring the likely chances of possible invasions from China due to the presence of the ancestral predominant haplotype H1h1. Moreover, the BMSB populations from Austria, Serbia and Slovenia were possibly of recent invasions from China. In conclusion, we believe that the novel haplotype information and the genetic diversity among the global BMSB populations will lay down the foundation for future population genomic studies and could help in formulating effective BMSB management strategies. This study will also help in tracing the origin of BMSB intercepted at the border in those countries, such as New Zealand, where the species has not yet established.


Sample collection and DNA extraction

BMSB specimens were collected from 43 regions/provinces in 12 countries (Austria, Chile, China, Georgia, Hungary, Italy, Japan, Romania, Serbia, Slovenia, Turkey, and the USA). All BMSB specimens were stored in 95% ethanol at − 20 °C and the species identity was confirmed by morphological characteristics by MPI entomologists. Total genomic DNA was extracted from each individual specimen using QIAGEN DNeasy® Blood & Tissue Kit with QIAGEN RNase A (Qiagen, Valencia, CA, the USA). The DNA quality and purity were determined using NanoDrop™ (CA, the USA) and quantified using QuantiFluor™ (CA, the USA) (Additional file 1).

Polymerase Chain Reaction (PCR) and sequencing

The genetic diversity of the BMSB populations was determined by analysing the mitochondrial DNA (mtDNA) cytochrome c oxidase subunit I (COI) and cytochrome c oxidase subunit II (COII). The two markers were chosen due to their fairly high variability and large number of sequences previously reported. Partial sequence of the COI (657 bp) and COII (518 bp) genes were amplified using the genomic DNA as a template. The LCO1490 (5′-GGTCAACAAATCATAAAGATATTGG-3′) and HCO2198 (5′-TAAACTTCAGGGTGACCAAAAAATCA-3′) [39] primer pairs were used for amplification of the COI region. Similarly, the HhalysCO2F2 (5′-TAACCCAAGATGCAAATTCT-3) and HhalysCO2R2 (5′-CCATATATAATTCCTGGACGA-3′) primer pairs were used for amplification of the COII region [29]. The PCR reagents (reaction volume of 20 μl) for COI and COII were the same except the primers, comprising of 4.4 μl sterile deionized water, 10.0 μl 2X GoTaq® Green Master Mix (Promega), 1.0 μl of 5 μM Forward primer, 1.0 μl of 5 μM Reserve primer, 0.6 μl of 50 mM MgCl2, 1.0 μl of 10 μg/μl BSA. Thermocycling conditions for both reactions comprised of an initial denaturation step of 94 °C for 5 mins, followed by 40 cycles of denaturation at 94 °C for 15 s, annealing at 50 °C for 30 s, and extension at 72 °C for 45 s, followed by a final extension phase at 72 °C for 7 min. All reactions were performed using a Veriti 96-well thermal cycler (Life Technologies). In this study, 463 specimens were used for COI and COII sequencing. The amplified DNA fragments in the final PCR product was evaluated on 1% agarose gel against a 100 bp DNA ladder (Invitrogen™) in TAE buffer stained with SYBR safe (Life Technologies) and visualised using a Gel Doc Software system (BioRad, Hercules, CA, the USA). The resulting product was diluted 5-fold with sterile water, and sent for Sanger sequencing Macrogen (Soul, South Korea) using the amplification primers in both directions. The quality of the Sanger sequencing dataset was manually examined and analysed in Geneious software (Biomatters, Auckland, New Zealand) [40]. The resulting quality-filtered COI and COII gene sequences for each haplotype were submitted to GenBank under the accession numbers, MT517228 - MT517274 for COI and MT490838 - MT490860 for COII. All the metadata for each specimen and their sequences were submitted into BOLD under project code BMSB. The BOLD processed IDs for each specimen are listed in Additional file 1. All the sequences were used to create a dataset in BOLD, under the DOI:

Public COI and COII sequence data acquisition and haplotype calling

A comparative barcoding analysis of COI and COII genes from the BMSB cohort of the current study and the publicly available COI and COII sequences from BMSB were conducted. The reference COI and COII sequences were created, respectively, by the following steps. Firstly, previously reported COI and COII barcode sequences of BMSB were downloaded from the GenBank. The COI (UID: 1674561291, 1,591,437,641, 1,334,761,755, 1,304,534,304, 1,240,496,350, 1,201,369,261, 1,024,298,892, 985,693,878, 443,298,673, 537,366,792, 552,099,040) and COII (UID: 1334762135, 552,098,974) sequences were aligned separately using Geneious software v10.2.5 [40]. Secondly, for COI, a 657 bp region was selected for further analysis while a 518 bp region was selected for COII. Sequences not from the same region of COI or COII were discarded. Finally, the aligned sequences were trimmed to 657 bp for COI and 518 for COII and duplicated sequences removed. The remaining unique sequences were used as the reference dataset for haplotype assignment.

All the COI and COII sequences obtained in this study were checked and edited in Geneious software [40] to remove poor quality sequences. A 657 bp COI and a 518 bp COII regions were used for further analysis, thus no missing data and ambiguous base-calls such as Y, W, N. S in the sequences obtained. Using the previously reported BMSB haplotypes as a reference, an in-house python script was developed to allocate the haplotype name for each obtained sequence. The script automatically allocated haplotype name to each individual of this study by searching the database. Sequences identical with those previously reported were assigned the same haplotype name (with a prefix H for COI or h for COII). The remaining sequences which did not share 100% sequence identity with the previously reported were given a new haplotype name (with a prefix N for COI and n for COII). Finally, a Microsoft Excel file was created (Additional file 1) and the number of haplotypes per population was calculated in Microsoft Excel. To further confirm those new haplotypes, all the sequences detected in this study were also BLAST searched in GenBank database [32]. To better illustrate the population genetics, the COI and COII sequences were concatenated into one linked longer haplotype, namely COI-COII combined haplotype. All the subsequent genetic analyses were conducted based on the COI-COII combined haplotypes.

Population genetic analysis

The genetic diversity as the percentage of each haplotype present from different regions and/or countries was estimated by calculating the number of haplotypes detected at the country, divided by the total number of individuals sampled using Microsoft Excel. Population genetic diversity, as indexed by the number of haplotype (Hn), haplotype diversity (Hd) and nucleotide diversity (π) was estimated using DnaSP v6 [41] to quantify the degree of genetic diversity. Analyses of haplotype and nucleotide diversity were conducted separately for each population as well as for populations in one country and one genetic group. To examine the historical demographic expansion, a neutrality test was performed under DnaSP v6. Based on the Hd values, the BMSB populations in the invaded countries were divided into two groups: genetic group A (Chile, Georgia, Hungary, Italy, Romania, Turkey and the USA) and genetic group B (Austria, Serbia and Slovenia). Therefore, the neutrality test was also performed among the BMSB populations from China, Japan, group A, and group B.

Genetic differentiation among the BMSB populations was estimated by the fixation index (FST), and the overall genetic variance was calculated by AMOVA (Analysis of molecular variance). Both calculations were fulfilled using Arlequin 3.5 [42]. To show the genetic relationships among the populations, a Principal Coordinates Analysis (PCoA) was conducted based on the FST data (Additional file 3) using GenAlEx 6.5 [43]. The relationships among the COI and COII combined haplotypes were examined using a parsimony network by applying the method described by TCS analysis [44] based on the COI-COII combined haplotypes and visualized using PopART [45].

Availability of data and materials

The in-house script can be found using the link below.!AvRFcQuxR5sIgpZzg_Ch02Lm2u4tJA?e=u9CoLx

The dataset for the COI and COII sequences in BOLD can be found under the DOI:



Analysis of molecular variance


Brown marmorated stink bug


Cytochrome c oxidase I


Cytochrome c oxidase II


Mitochondrial DNA


Ministry for Primary Industries


New Zealand Institute of Economic Research


Principle Coordinates Analysis


  1. MacLellan R. Plants and environment brown marmorated stink bug: a potential risk to New Zealand. Surveillance (Wellington). 2013;40(1):34–6.

    Google Scholar 

  2. Nielsen AL, Hamilton GC. Life history of the invasive species Halyomorpha halys (Hemiptera: Pentatomidae) in northeastern United States. Ann Entomol Soc Am. 2009;102(4):608–16.

    Article  Google Scholar 

  3. Leskey TC, Hamilton GC, Nielsen AL, Polk DF, Rodriguez-Saona C, Bergh JC, Herbert DA, Kuhar TP, Pfeiffer D, Dively GP. Pest status of the brown marmorated stink bug, Halyomorpha halys in the USA. Outlooks Pest Manag. 2012;23(5):218–26.

    Article  Google Scholar 

  4. Hoebeke ER, Carter ME. Halyomorpha halys (Stål) (Heteroptera: Pentatomidae): a polyphagous plant pest from Asia newly detected in North America. Proc Entomol Soc Wash. 2003;105(1):225–37.

    Google Scholar 

  5. Leskey TC, Nielsen AL. Impact of the invasive brown marmorated stink bug in North America and Europe: history, biology, ecology, and management. Annu Rev Entomol. 2018;63:599–618.

    Article  CAS  Google Scholar 

  6. Wiman NG, Walton VM, Shearer PW, Rondon SI, Lee JC. Factors affecting flight capacity of brown marmorated stink bug, Halyomorpha halys (Hemiptera: Pentatomidae). J Pest Sci. 2015;88(1):37–47.

    Article  Google Scholar 

  7. Lee D, Nielsen AL, Leskey TC. Dispersal capacity and behavior of nymphal stages of Halyomorpha halys (Hemiptera: Pentatomidae) evaluated under laboratory and field conditions. J Insect Behav. 2014;27(5):639–51.

    Article  Google Scholar 

  8. Hamilton GC. Brown marmorated stink bug. Am Entomol. 2009;55(1):19–20.

    Article  Google Scholar 

  9. Josifov MV, Kerzhner IM. Heteroptera aus Korea. II. Teil (Aradidae, Berytidae, Lygaeidae, Pyrrhocoridae, Rhopalidae, Alydidae, Coreidae, Urostylidae, Acanthosomatidae, Scutelleridae, Pentatomidae, Cydnidae, Plataspidae). Fragmenta Faunistica. 1978;23(9):137–96.

    Google Scholar 

  10. Rider DA, Zheng LY, Kerzhner IM. Checklist and nomenclatural notes on the Chinese Pentatomidae (Heteroptera). II Pentatominae Zoosystematica Rossica. 2002;11(1):135–53.

    Google Scholar 

  11. Lee D, Short BD, Joseph SV, Bergh JC, Leskey TC. Review of the biology, ecology, and management of Halyomorpha halys (Hemiptera: Pentatomidae) in China, Japan, and the Republic of Korea. Environ Entomol. 2013;42(4):627–41.

    Article  Google Scholar 

  12. EPPO: Halyomorpha halys (HALYHA). EPPO Global database, Paris: EPPO 2020. (Accessed 3 August 2020).

  13. Gariepy TD, Haye T, Fraser H, Zhang J. Occurrence, genetic diversity, and potential pathways of entry of Halyomorpha halys in newly invaded areas of Canada and Switzerland. J Pest Sci. 2014;87(1):17–28.

    Article  Google Scholar 

  14. Macavei LI, Bâeţan R, Oltean I, Florian T, Varga M, Costi E, Maistrello L. First detection of Halyomorpha halys Stål, a new invasive species with a high potential of damage on agricultural crops in Romania. Lucrări Ştiinţifice, seria Agronomie. 2015;58(1):105–8.

    Google Scholar 

  15. Heckmann R. First evidence of Halyomorpha halys (Stal, 1855) (Heteroptera: Pentatomidae) in Germany. Heteropteron. 2012;36:17–8.

    Google Scholar 

  16. Vétek G, Papp V, Haltrich A, Rédei D. First record of the brown marmorated stink bug, Halyomorpha halys (Hemiptera: Heteroptera: Pentatomidae), in Hungary, with description of the genitalia of both sexes. Zootaxa. 2014;3780(1):194–200.

    Article  Google Scholar 

  17. Maistrello L, Dioli P, Vaccari G, Nannini R, Bortolotti P, Caruso S, Costi E, Montermini A, Casoli L, Bariselli M: First records in Italy of the Asian stinkbug Halyomorpha halys, a new threat for fruit crops. Atti, Giornate Fitopatologiche, Chianciano Terme (Siena), 18–21 marzo 2014, Volume primo. 2014;283–288.

  18. Wermelinger B, Wyniger D, Forster B. First records of an invasive bug in Europe: Halyomorpha halys Stål (Heteroptera: Pentatomidae), a new pest on woody ornamentals and fruit trees? Mitteil Schw Entomol Gesellschaft. 2008;81:1–8.

    Google Scholar 

  19. Milonas PG, Partsinevelos GK. First report of brown marmorated stink bug Halyomorpha halys Stål (Hemiptera: Pentatomidae) in Greece. EPPO Bul. 2014;44(2):183–6.

    Article  Google Scholar 

  20. Cesari M, Maistrello L, Ganzerli F, Dioli P, Rebecchi L, Guidetti R. A pest alien invasion in progress: potential pathways of origin of the brown marmorated stink bug Halyomorpha halys populations in Italy. J Pest Sci. 2015;88(1):1–7.

    Article  Google Scholar 

  21. Faúndez EI, Rider DA. The brown marmorated stink bug Halyomorpha halys (Stål, 1855) (Heteroptera: Pentatomidae) in Chile. Arquivos Entomolóxicos. 2017;17:305–7.

    Google Scholar 

  22. Zhu G, Bu W, Gao Y, Liu G. Potential geographic distribution of brown marmorated stink bug invasion (Halyomorpha halys). PLoS One. 2012;7(2):e31246.

    Article  CAS  Google Scholar 

  23. Kriticos DJ, Kean JM, Phillips CB, Senay SD, Acosta H, Haye T. The potential global distribution of the brown marmorated stink bug, Halyomorpha halys, a critical threat to plant biosecurity. J Pest Sci. 2017;90(4):1033–43.

    Article  Google Scholar 

  24. Vandervoet TF, Bellamy DE, Anderson D, MacLellan R. Trapping for early detection of the brown marmorated stink bug, Halyomorpha halys, in New Zealand. N Z Plant Protect. 2019;72:36–43.

    Google Scholar 

  25. Labware LIMS. Ministry for Primary Industries. 2020. (Accessed June 2020).

  26. Ballingall J, Pambudi D. Quantifying the economic impacts of a Brown Marmorated stink bug incursion in New Zealand: a dynamic computable general equilibrium modelling assessment; 2017.

    Google Scholar 

  27. Ficetola GF, Bonin A, Miaud C. Population genetics reveals origin and number of founders in a biological invasion. Mol Ecol. 2008;17(3):773–82.

    Article  CAS  Google Scholar 

  28. Lee W, Guidetti R, Cesari M, Gariepy TD, Park Y, Park C. Genetic diversity of Halyomorpha halys (Hemiptera, Pentatomidae) in Korea and comparison with COI sequence datasets from East Asia, Europe, and North America. Fla Entomol. 2018;101(1):49–54.

    Article  Google Scholar 

  29. Xu J, Fonseca DM, Hamilton GC, Hoelmer KA, Nielsen AL. Tracing the origin of US brown marmorated stink bugs, Halyomorpha halys. Biol Invasions. 2014;16(1):153–66.

    Article  CAS  Google Scholar 

  30. Cesari M, Maistrello L, Piemontese L, Bonini R, Dioli P, Lee W, Park C, Partsinevelos GK, Rebecchi L, Guidetti R. Genetic diversity of the brown marmorated stink bug Halyomorpha halys in the invaded territories of Europe and its patterns of diffusion in Italy. Biol Invasions. 2018;20(4):1073–92.

    Article  Google Scholar 

  31. Haye T, Gariepy T, Hoelmer K, Rossi J, Streito J, Tassus X, Desneux N. Range expansion of the invasive brown marmorated stinkbug, Halyomorpha halys: an increasing threat to field, fruit and vegetable crops worldwide. J Pest Sci. 2015;88(4):665–73.

    Article  Google Scholar 

  32. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10.

    Article  CAS  Google Scholar 

  33. Ratnasingham S, Hebert PDN. BOLD: The Barcode of Life Data System ( Mol Ecol Notes. 2007;7:355–64.

  34. Deardorff A. Tableau (version. 9.1). J Med Lib Assoc. 2016;104(2):182–3.

    Article  Google Scholar 

  35. Zhu G, Ye Z, Du J, Zhang D, Zhen Y, Zheng C, Zhao L, Li M, Bu W. Range wide molecular data and niche modeling revealed the Pleistocene history of a global invader (Halyomorpha halys). Sci Rep. 2016;6(1):1–10.

    Article  Google Scholar 

  36. Wermelinger BEAT, Wyniger D, Forster BEAT. First records of an invasive bug in Europe: Halyomorpha halys Stal (Heteroptera: Pentatomidae), a new pest on woody ornamentals and fruit trees? Schweiz Entomol Ges. 2008;81(1/2):1.

    Google Scholar 

  37. Güncan A, Gümüş E. Brown marmorated stink bug, Halyomorpha halys (Stål, 1855) (Hemiptera: Heteroptera, Pentatomidae), a new and important pest in Turkey. Entomol News. 2019;128(2):204–10.

    Article  Google Scholar 

  38. Fu Y. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147(2):915–25.

    Article  CAS  Google Scholar 

  39. Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R. Phylogenetic uncertainty. Mol Mar Biol Biotechnol. 1994;3:294–9.

    CAS  PubMed  Google Scholar 

  40. Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, Buxton S, Cooper A, Markowitz S, Duran C. Geneious basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28(12):1647–9.

    Article  Google Scholar 

  41. Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, Sánchez-Gracia A. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol. 2017;34(12):3299–302.

    Article  CAS  Google Scholar 

  42. Excoffier L, Lischer HE. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and windows. Mol. Ecol. Resour. 2010;10(3):564–7.

    Google Scholar 

  43. PE PRS: GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research – an update. Bioinformatics. 2012;28(19):2537.

    Article  Google Scholar 

  44. Clement M, Snell Q, Walker P, Posada D, Crandall K. TCS: estimating gene genealogies. In proceedings of the 16th international parallel and distributed processing symposium (IPDPS '02), vol. 311: IEEE computer society; 2002.

  45. Leigh JW. Bryant D: popart: full-feature software for haplotype network construction. Methods Ecol Evol. 2015;6(9):1110–6.

    Article  Google Scholar 

Download references


This work is dedicated to the memory of the late Dr Vétek Gábor, a co-author of this paper.

We would like to thank Drs David Waite and Luciano Rigano from the Plant Health and Environment Laboratory (PHEL), Ministry for Primary Industries (MPI) for critically reviewing the manuscript. We would like to thank Drs Zhidong Yu, Claire McDonald and Sue Escott-Brown for their support, suggestions, and feedback to the project. Our special thanks go to the Entomology team members of the PHEL of MPI for supplying the BMSB specimens intercepted at the New Zealand border or from the post-border surveillance and assisting in morphological identification of the specimens. Finally, we would like to thank the following researchers for providing the BMSB specimens from a range of different countries:

Dr. Catherine Duthie of Readiness Programmes, MPI, New Zealand; Ms. Xinyao Gu from Institute of Entomology, Guizhou University, China; Ms. Jianyun Wang from Environment and Plant Protection Institute, Chinese Academy of Tropical Agriculture sciences, China; Dr. Shu-Jun Wei from Institute of Plant and Environmental Protection, Beijing Academy of Agriculture and Forestry Sciences, China; Dr. Jun Ma from Institute of Plant Quarantine, Guangzhou Customs Technology Center, China; Dr. Yongliang Fan from Northwest A & F University, China; Dr. Dong Liu from Chines Academy of Sciences, China; Dr. Xiaofen Xue from Nanjing Agricultural University, China; Dr. Lixia Xie from Shandong Agricultural University, China; Dr. Gonzalo Avila from Plant & Food Research, New Zealand; Dr. Kyo Itoyama from Meiji University, Japan; Valerio Mazzoni from Edmund Mach Foundation, Iris Bernardinelli from Servizio fitosanitario e chimico, ricerca, sperimentazione e assistenza tecnica, Italy; Drs Laura Nixon and Kim Hoelmer from US Department of Agriculture, the USA; Dr. Roxana Ciceoi from University of Agricultural Sciences and Veterinary Medicine of Bucharest, Romania; Dr. Ilania Astorga Leiva from Subdepto Vigilancia y control de plagas agrícolas, División Protección Agrícola y Forestal, Servicio Agrícola y Ganadero, Gobierno de Chile; Dr. Attila Torma from University of Szeged, Dr. Ákos Varga from Government Office of Pest County, Hungary; Dr. Temel Göktürk from Artvin Coruh University, Turkey; Dr. Maka Murvanidze from Agricultural University of Georgia, Georgia; Dr. Wolfgang Rabitsch from University of Vienna, Austria.

This research was partly supported by the Hungarian Ministry for Innovation and Technology within the framework of the Higher Education Institutional Excellence Program (NKFIH-1159-6/2019) in the scope of plant breeding and plant protection research of Szent István University, and the NKTH 2017-2.3.3-TÉT-VN-2017-00006 (Biological control of invasive pest species in Vietnam and Hungary) research project. Finally, we would like to thank the editor and two anonymous reviewers for their suggestions in improving the manuscript.


The research was funded by the Operational Research programme under the project number 405731 from the Ministry for Primary Industries (MPI), New Zealand. The funding body provided the financial supports for conducting the experiments and generating the data, and the costs for publishing the findings. The funding body played no role in the design of the study and collection, analysis, and interpretation of data and in the preparation of the manuscript.

Author information

Authors and Affiliations



DL, JY, CP, SG and RKB conceived and designed the study. QF, DA, DG, JZ and GV collected the samples and conducted the morphological identification. JY conducted the laboratory experiments, bioinformatics, and statistical analyses. CP assisted the bioinformatics analysis. JY drafted the manuscript with input from DL, CP, GV, PF, AB. All authors have read, edited, and approved the final version of the manuscript.

Corresponding author

Correspondence to Dongmei Li.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declared no conflicts of interest to this work.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

Metadata information of the collected BMSB cohort.

Additional file 2.

Haplotype information of the BMSB populations.

Additional file 3.

FST values of all the BMSB populations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Yan, J., Pal, C., Anderson, D. et al. Genetic diversity analysis of brown marmorated stink bug, Halyomorpha halys based on mitochondrial COI and COII haplotypes. BMC Genom Data 22, 7 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: