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

doi:10.21203/rs.3.rs-59349/v3

Background: 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 on 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.

Result: 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 exist between groups of populations: China, Japan, and the invasive populations.

Conclusion: 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 and a new stable genetic cluster has established in Europe and the USA. The study also revealed that 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, especially biocontrol.

Medical Genetics

Molecular Genetics

Population Genetics

Invasion

Genetic cluster

mitochondrial DNA

nucleotide diversity

haplotype diversity

Brown marmorated stink bug (BMSB), Halyomorpha halys (Stål, 1855) (Hemiptera: Pentatomidae) is a highly polyphagous insect, which feeds on more than 300 hosts plant [1, 2]. BMSB has been causing extensive damage to a wide variety of agricultural crops and poses a global economic threat [3, 4]. As per 2010 reports, the economic lose 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]. BMSB is considered as a nuisance pest as adults search for man-made structures to overwinter and discharge an unpleasant and long-lasting odour once disturbed [8].

BMSB is native to China (including Taiwan), Japan and the Korean peninsula [9-11] and has invaded to 30 countries [40] including most of the states in the USA [4, 5], Canada [12], Europe [12-19] and Chile [20]. Climate modelling and the known global distribution of BMSB categorizes New Zealand at high risk of establishment should a successful invasion occur [21, 22]. The climatic suitability for this organism is wide, ranging from 30 to 50 degrees of latitude, with an annual mean temperature range between 10 ℃ and 30 ℃ (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 [21, 22]. BMSB has increasingly been intercepted at the border and post-border scenarios due to the rise of international travel and trade [23] since 2005 when BMSB was first intercepted at border of New Zealand [41]. There were 1,620 recorded interceptions of BMSB at border of New Zealand since 2005 [41]. As agricultural exports comprise a significant proportion of the New Zealand Gross Domestic Product (GDP), such an event would be highly detrimental to the country, and BMSB is considered a “high risk” organism. The New Zealand Institute of Economic Research (NZIER) study predicts a worst-case scenario of 3.6 billion New Zealand dollars in agricultural losses by 2038 if BMSB establishes in New Zealand successfully [24]. Therefore, it is imperative to study the genetic diversity of BMSB to provide a method to trace the origin of captured individuals to better maintain border security and could also play a vital role in developing an appropriate pest management strategy, especially the biocontrol.

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 [25], and has been applied in a number of studies for BMSB [12, 19, 26-28]. Previous studies on the 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 invaded BMSB populations in the USA [27]. Additional studies, focusing on BMSB populations in Europe, used Cytochrome Oxidase I (COI) gene alone [12], or in combination with the COII gene [19, 28], and revealed additional information on the genetic diversity of BMSB [12, 38, 39]. To expand the knowledge on the genetic diversity of BMSB populations around the world, we collected BMSB 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 countries of Europe. In the current study, two mitochondrial barcoding genes, COI and COII, were sequenced from 463 BMSB specimens collected from 43 regions 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 shared 100% identity in the same region (657 bp) with those previous reported, the same haplotype names were given while new names were given accordingly for the rest sequences obtained in this study. All the new haplotypes identified were confirmed by BLAST. The result showed all the new haplotypes were unique in the region (657 bp). Further analysis showed that 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 is not accurate to assign the same haplotype name if the sequences are not the same length.

Further comparison with the COI sequences from those deposited in Barcode of Life Data System (BOLD) showed that all the sequence obtained in this study are belong to BIN AAM9563, with over 98% identity. In contrast, the sequences shared over 94% and 82% to the sequences of the 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), 29 haplotypes were identified including 20 novel and 9 previously reported ones (Additional file 2). BLAST search showed that the new haplotypes did not share any identical sequences in the 518 bp overlap region with that of the previous reported sequences.

The geographical distribution of the identified COI and COII haplotypes are shown in Figure 1 and Figure 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 (China and Japan) except N47 (Slovenia) (Additional file 2). All the novel COII haplotypes identified were detected from the two native countries, China and Japan (Additional file 2).

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 H45, N22 , H23 , N40 for COI and h11 for COII (Table 1). Outside of the native areas, 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), thus used for COI-COII combined haplotype analysis. The combined COI-COII haplotype analysis produced 59 haplotypes, in which only 5 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 Figure 3. In the native countries of BMSB, high haplotype diversity was observed with 24 haplotypes in China and 32 in Japan, without haplotypes shared between the two countries (Additional file 2 and Figure 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 the BMSB- 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).

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 for further 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, 5 sampling sites were studied, of which at two sites no haplotype diversity was observed, while 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. It indicates that the invasion of BMSB in Hungary may have come from genetically distinct populations.

In 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 4) 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 haplotypes diversity was slightly higher with an average of 0.63 for group B (Austria, Serbia, and Slovenia), but a Fu's Fs values of 1.453 (p > 10).

The Principle Coordinates analysis (PCoA) using the F_ST 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) (Figure 4, Additional file 3). The recent invasion in Slovenia showed genetic similarities to those from Hebei and Beijing provinces of China. 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.

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

The haplotype network of the BMSB individuals further revealed the widespread occurrence of H1h1 and H3h3, except the population from Japan excluded, whereas all the other haplotypes were mainly detected in the native countries (Figure 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 the combined haplotypes, the combined COI and COII dataset from the present and the previous studies [19, 28] were analysed together and resulted in a total of 81 haplotypes. The haplotype network analysis (Figure 6) indicated similar genetic relationships as previously reported except that few BMSB specimens from Italy had close relationship with Japanese populations (Figure 6).

This study revealed 51 COI haplotypes (36 novel) and 29 COII haplotypes (20 novel) from 463 BMSB individual 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 population 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 (2014) [27] 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 (2018) [26] identified COI haplotype H1 as the main COI haplotype for China (68%), Hungary (98.8%), Italy (80%) and the USA (92.5%). Similarly, the current study identified haplotype H1 as the predominant one, accounting for 27.6% of Chinese samples, 93.4% of Hungarian samples, 95.8% of 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 [27] but was not detected in this study (Table 2). Variations in the haplotype numbers and percentages between the current study and that of Cesari et al [28] were also observed. Cesari et al. (2018) [28] reported a total of twenty-six COI and COII combined haplotypes for BMSB specimens, mainly from Italy [28]. 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 81 known so far. Cesari et al. (2018) [28] 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 [28]. 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 [28] were also detected in China, Austria, and Japan, respectively, in this study. After combining the haplotype data from the two studies, haplotype network (Figure 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 [19, 28]. Therefore, it can be hypothesised that the BMSB populations in Italy possibly have originated in Asia, which can be supported by the extensive, ongoing cross-border travel and trade between Asia and Italy. However, the invasion from North America cannot be ruled out as H1h1 was found also in the USA. Furthermore, haplotype H8 was detected in Switzerland [12] in 2012, in France [29] and Northern Italy in 2013 [19], suggesting the possible invasion of BMSB to Italy was from Switzerland [28] based on the geographical relatedness. 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 the borders. However, the widespread distribution of H1 and H3 haplotypes in Austria (Table 1), opens the likelihood 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 study by Xu et al (2014) [27]. 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 in congruent with previous studies [19, 27]. The haplotype network highlighted that China and Japan had their own haplotype clusters, suggesting that there is limited or no interbreeding due to geographical barriers between these BMSB populations and that Chinese and could be evolving independently. The Principal Coordinates Analysis also supported this conclusion. According to the Principal Coordinates analysis, most of the Japanese and Chinese populations were clustered by themselves except populations from regions of Akita (JP) and Hainan (CN). This phenomenon also was found in Zhu’s study [39] that it believed the discontinuous distribution pattern between Japan and Hainan could be due to long-distance dispersal itself or human-mediated transportation.

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 [30]. The relatively low differences in F_ST value between the populations within the genetic group A indicates 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 the BMSB, where the 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 [4], then has spread in Europe and Chile recently. This study further confirmed 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 [30]. The haplotype diversity (Hd) and nucleotide diversity (π) of these three populations in group B was 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 exist of the common haplotype H1h1 and H3h3 in China and the invaded countries. The reasons for that need further investigation, one possibility is 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 also not clear whether there is a possibility that these BMSB haplotypes can adapt more easily to new environments. Since the COI and COII haplotype analysis is 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 the 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 invasion into Chile, Hungary, Georgia, Turkey, Romania and Italy potentially from the USA without ignoring the likely chances of possible invasion from China due to the presence of the ancestral predominant haplotype H1h1. Moreover, the BMSB populations from Austria, Serbia and Slovenia were possibly recent invasions from China. In conclusion, we believe that the novel haplotype information and enhanced knowledge of genetic diversity that have been observed among the global BMSB populations will lay down the foundation for future population genomic studies and could help in formulating an effective BMSB management strategies, especially in the biocontrol arena largely help future population genomic studies. This study will also help in tracing the origin of BMSB intercepted at the border in those countries, such as New Zealand, where they never established earlier.

Sample collection and DNA extraction

BMSB specimens were collected from 43 regions 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 ℃ 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, Table 4).

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 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 template. The LCO1490 (5’-GGTCAACAAATCATAAAGATATTGG-3’) and HCO2198 (5’-TAAACTTCAGGGTGACCAAAAAATCA-3’) [31] 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 [27]. 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 50mM 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 seconds, annealing at 50°C for 30 seconds, and extension at 72 °C for 45 seconds, followed by a final extension phase at 72 °C for 7 minutes. 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 (Biomatters, Auckland, New Zealand) [32]. 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 process 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 following steps. Firstly, previously reported COI and COII barcode sequences of BMSB were downloaded from the GenBank. The COI (UID: 1674561291, 1591437641, 1334761755, 1304534304, 1240496350, 1201369261, 1024298892, 985693878, 443298673, 537366792, 552099040) and COII (UID: 1334762135, 552098974) sequences were aligned separately using Geneious 10.2.5 [32]. Secondly, for COI, a 657 bp region was selected for further analysis while 518 bp selected for COII. The sequences were not from the same region for COI or COII were discarded. Finally, the aligned sequences were trimmed to 657 bp for COI and 518 for COII, then removed the duplicated sequences. The reminding unique sequences were the reference sequences for haplotype assignment. All the sequences obtained in this study were aligned with the reference sequences. All sequences obtained in this study were checked and edited in Geneious. The poor quality sequences were discarded, thus all sequences used in this paper have no missing data and ambiguous base-calls such as Y, W, N. A. A database of the haplotypes that have been reported was created as 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. The sequences are identical with those previously reported were assigned the same haplotype name (with a prefix ‘H for COI or h for COII’). The rest sequences which not shared 100% identity with the previous reported was given a new haplotype (with a prefix ‘N for COI and n for COII’). Finally, a excel file was created as additional file 1 and the number of haplotypes per population was calculated using excel. To further confirm those new haplotypes, all the sequences detected in this study were also BLAST searched in GenBank database [42]. 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 examined by calculating the number of haplotypes detected at the region/country, divided by the total number of individuals sampled using Excel. Population genetic diversity, as indexed by the number of haplotype (Hn), haplotype diversity (Hd) and nucleotide diversity (π) was performed using DnaSP v6 [33] 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 performed among the BMSB populations from China, Japan, group A and group B.

Genetic differentiation among populations was estimated by the fixation index (F_ST), and the overall genetic variance was estimated by AMOVA (Analysis of molecular variance). Both calculations were fulfilled using Arlequin 3.5 [34]. To show the relationships among the populations, a PCoA (Principal Coordinates Analysis) was conducted based on F_ST data (Additional file 3) using GenAlEx 6.5 [35]. The relationships among haplotypes were estimated using a parsimony network by applying the method described by TCS analysis [36] based on the COI-COII combined haplotypes and visualized using PopART [37]. To better illustrate the genetic relationships among a wide region of BMSB samples, the data obtained in this study were also combined with the reference sequences from the GenBank for those specimens containing both COI and COII sequences.

Analysis of molecular variance (AMOVA)

Brown marmorated stink bug (BMSB)

Cytochrome c oxidase I (COI)

Cytochrome c oxidase II (COII)

Mitochondrial DNA (mtDNA)

Ministry for Primary Industries (MPI)

New Zealand Institute of Economic Research (NZIER)

Principle Coordinates Analysis (PCoA)

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and materials

The in-house script can be found using the link below

https://1drv.ms/u/s!AvRFcQuxR5sIgpZzg_Ch02Lm2u4tJA?e=u9CoLx

Competing interests

The authors declared no conflicts of interest to this work

Funding

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’s contributions

DL, JY, CP, SG and RKB conceived and designed the study. QF, DA, DG, JZ, RG 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.

Acknowledgments

We would like to thank Drs David Waite and Luciano Rigano from 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 for 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 Fond 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 researches 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.

MacLellan R: Plants and environment brown marmorated stink bug: a potential risk to New Zealand. Surveillance (Wellington). 2013;40(1).
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-616.
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 on Pest Management. 2012;23(5):218-226.
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-237.
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.
Wiman NG, Walton VM, Shearer PW, Rondon SI, Lee JC: Factors affecting flight capacity of brown marmorated stink bug, Halyomorpha halys (Hemiptera: Pentatomidae). Journal of pest science. 2015;88(1):37-47.
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-651.
Hamilton GC: Brown marmorated stink bug. Am Entomol. 2009;55(1):19-20.
JOSIFOV M, IM K: Heteroptera aus korea. Ii. Aradidae, berytidae, lygaeidae, pyrrhocoridae, rhopalidae, alydidae, coreidae, urostylidae, acanthosomatidae, sautelleridae, pentatomidae, cydnidae, plataspidae. 1978.
Rider DA, Zheng LY, Kerzhner IM: Checklist and nomenclatural notes on the Chinese Pentatomidae (Heteroptera). II. Pentatominae. Zoosystematica Rossica. 2002;11(1):135-153.
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-641.
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. Journal of pest science. 2014;87(1):17-28.
Macavei LI, Baetan 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. Lucrări Ştiinţifice, Universitatea de Stiinte Agricole Şi Medicină Veterinară” Ion Ionescu de la Brad” Iaşi, Seria Agronomie. 2015;58:105-108.
Heckmann R: First evidence of Halyomorpha halys (Stal, 1855) (Heteroptera: Pentatomidae) in Germany. Heteropteron. 2012;36:17-18.
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.
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.
Wermelinger B, Wyniger D, Forster B: First records of an invasive bug in Europe: Halyomorpha halys Stal (Heteroptera: Pentatomidae), a new pest on woody ornamentals and fruit trees? Mitteilungen-Schweizerische Entomologische Gesellschaft. 2008;81(1/2):1.
Milonas PG, Partsinevelos GK: First report of brown marmorated stink bug Halyomorpha halys Stål (Hemiptera: Pentatomidae) in Greece. EPPO Bulletin. 2014;44(2):183-186.
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. Journal of Pest Science. 2015;88(1):1-7.
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-307.
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.
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. Journal of Pest Science. 2017;90(4):1033-1043.
Vandervoet TF, Bellamy DE, Anderson D, MacLellan R: Trapping for early detection of the brown marmorated stink bug, Halyomorpha halys, in New Zealand. New Zealand Plant Protection. 2019;72:36-43.
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.
Ficetola GF, Bonin A, Miaud C: Population genetics reveals origin and number of founders in a biological invasion. Mol Ecol. 2008;17(3):773-782.
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.
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-166.
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-1092.
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. Journal of Pest Science. 2015;88(4):665-673.
Fu Y: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147(2):915-925.
Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R: Phylogenetic uncertainty. Mol Marine Biol Biotechnol. 1994;3:294-299.
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-1649.
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-3302.
Excoffier L, Lischer HE: Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular ecology resources. 2010;10(3):564-567.
PE PRS: GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research–an update. Bioinformatics. 2012;28(19):2537.
Anonymous Proceedings of the Parallel and Distributed Processing Symposium, International: 2002.
Leigh JW, Bryant D: popart: full‐feature software for haplotype network construction. Methods in Ecology and Evolution. 2015;6(9):1110-1116.
Josifov M, Im K: Heteroptera aus Korea. II. Teil (Aradidae, Berytidae, Lygaeidae, Pyrrhocoridae, Rhopalidae, Alydidae, Coreidae, Urostylidae, Acanthosomatidae, Scutellieridae, Pentatomidae, Cydnidae, Plataspidae). Fragmenta Faunistica. 1978;23(9): 137–196.
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). Scientific reports. 2016;6(1):1-10.
EPPO: Halyomorpha halys (HALYHA). EPPO Global database, Paris, France: EPPO 2020. https://gd.eppo.int/taxon/HALYHA/distribution (accessed 3 August 2020)
Labware LIMS. Ministry for Primary Industries. 2020. (Accessed June 2020)
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ: Basic local alignment search tool. J Mol Biol. 1990;215(3):403-410.
Deardorff A: Tableau (version. 9.1). Journal of the Medical Library Association. 2016;104(2):182-183.

Table 1. The percentage of the dominant mtDNA COI and COII haplotypes of H. halys detected in 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 haplotypes in the country are listed in the table.

Percentage (%)	COI								COII
Percentage (%)	H1	H3	H22	H33	H23	H45	N22	N40	h1	h3	h11
China	28	10	14	22	-	-	-	-	26	46	-
Japan	-	-	-	-	16	16	14	9	-	-	81
Austria	25	50	-	-	-	-	-	-	20	67	-
Serbia	50	50	-	-	-	-	-	-	44	56	-
Slovenia	69	25	-	-	-	-	-	-	67	33	-
Chile	97	3	-	-	-	-	-	-	97	3	-
Georgia	100	-	-	-	-	-	-	-	100	-	-
Hungary	93	5	-	-	-	-	-	-	94	4	-
Italy	96	4	-	-	-	-	-	-	93	7	-
Romania	100	-	-	-	-	-	-	-	100	-	-
Turkey	100	-	-	-	-	-	-	-	100	-	-
the USA	100	-	-	-	-	-	-	-	100	-	-

Table 2. The sample information and mtDNA diversity. The total sample size (N) for each province and country are 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 and each genetic group. Asterisk (*) represents statistically significant difference (p<0.02).

Country	Province	N	Hn	Hd	π	Fu's Fs
China	Total	90	24	0.858	0.00327	-7.852*
	Anhui	5	3	0.7	0.00102
	Beijing	27	7	0.553	0.00113
	Hainan	7	5	0.857	0.00284
	Hebei	8	4	0.75	0.00225
	Jiling	3	2	0.667	0.00227
	Shannxi	31	9	0.751	0.00205
	Shanxi	8	2	0.429	0.00182
	Guizhou	1	1	NA	NA
Japan	Total	65	32	0.942	0.00238	-29.707*
	Akita	2	2	1	0.0017
	Chiba	2	2	1	0.0017
	Gifu	1	1	NA	NA
	Ibaraki	2	1	0	0
	Ishikawa	2	2	1	0.00255
	Iwate	4	4	1	0.00426
	Kagoshima	5	3	0.7	0.00102
	Kanagawa	3	3	1	0.00454
	Kyoto	10	5	0.756	0.00127
	Mie	5	4	0.9	0.00136
	Miyagi	3	3	1	0.0034
	Nagasaki	1	1	NA	NA
	Saga	7	4	0.857	0.00105
	Shizuoka	5	5	1	0.00221
	Tokushima	1	1	NA	NA
	Yamanashi	12	10	0.955	0.00253
Serbia	Senta	9	2	0.556	0.00095	2.302
Slovenia	Ljubljana	15	3	0.514	0.00115	1.626
Austria	Vienna	15	4	0.686	0.00206	1.84
Turkey	Arhavi	11	1	0	0	0
the USA	Total	24	1	0	0	0
	Maryland	15	1	0	0
	West Virginia	9	1	0	0
Georgia	Total	31	1	0	0	0
	Eki	28	1	0	0
	Samegrelo	3	1	0	0
Hungary	Total	90	3	0.107	0.00028	-0.195*
	Budapest	61	3	0.038	0
	Debrecen	10	1	0	0
	Pécs	6	1	0	0.00085
	Szeged	2	2	0.5	0.00136
	Szombathely	11	3	0.345	0
Italy	Total	24	2	0.083	0.00014	-0.192
	Codroipo (UD)	4	1	0	0
	Mantova	2	1	0	NA
	Pozzuolo del Friuli (UD)	17	1	0	0
	Trentino Alto	1	1	NA	0
Romania	Bucharest	23	1	0	0	0
Chile	Santiago	31	2	0.065		-0.426
group A		234	3	0.059	0.00045	-1.174*
group B		39	5	0.63	0.00158	1.453

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

Status:

Journal Publication

Version 3

Abstract

Figures

Background

Results

COI and COII haplotypes in the BMSB populations

COI-COII combined haplotypes of the BMSB populations

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

Discussion

Conclusions

Methods