Invasion History of Sirex noctilio Based on COI Sequence: The First Six Years in China

Sirex noctilio F. (Hymenoptera: Siricidae: Siricinae), a new invasive species in China, is a significant international forestry pest which, transported via logs and related wood packing materials, has led to environmental damage and substantial economic loss in many countries around the world. It was first detected in China in 2013, and since then infestations have been found in 18 additional sites. Using a 322 bp fragment of the mitochondrial barcode gene COI, we studied the genetic diversity and structure of S. noctilio populations in both native and invaded ranges, with a specific focus in China. Twelve haplotypes were found across the native and invaded distribution of the pest, of which three were dominant; among these there were only one or two mutational steps between each pair of haplotypes. No obvious genetic structure was found other than in Chinese populations. China has a unique and dominant haplotype not found elsewhere, and compared with the rest of the world, the genetic structure of Chinese populations suggested a multiple invasion scenario.


Introduction
Biological invasions are increasing exponentially as globalization progresses, without any sign of saturation, especially in insects [1,2]. At present, phytophagous species associated with woody plants constitute the majority of these insect invaders [3]. The problem warrants increased attention, as the invasion of alien species is likely to impact native biodiversity, and can also cause huge economic losses [4][5][6][7]. There is a statistically significant correlation between the number of major invasive alien species in China and Chinese trade imports from the year 1980 to 2015 (Pearson correlation r = 0.915, p < 0.01, N = 35). The current primary route for alien forest species invading China is through human activity. This is mainly unintentional introduction; transported seedlings and their propagation materials carry 46% of human activity-generated invaders, with packaging materials or wood accounting for 21% (National Bureau of Statistics, http://www.stats.gov.cn/tjsj/tjgb/ ndtjgb/) [8]. Some major invasive species, such as the North American bark beetle Dendroctonus valens, fall webworm Hyphantria cunea and loblolly pine mealybug Oracella acuta, have already caused the destruction of forestry resources in China, negatively influencing species diversity and ecological security, and resulting in major economic loss [9][10][11]. In order to monitor and manage the invasion of a pest, including possible biocontrol by natural enemies as well as the prevention of further introductions, it is essential to figure out its invasion history. However, observational records are usually difficult to get, and often biased, particularly when considering a recently-arrived invasive species. From this perspective, genetic analyses comparing populations in the native range with those in the invaded areas constitute an important tool for reconstructing the invasion history. for the first time 778 individuals from 6 continents. In particular, samples from original European populations and recent Chinese infestations were included. Overall, this study aims to elucidate the invasive routes of this pest as well as provide data enabling the prevention of the harmful consequences of S. noctilio spread.

Sample Collection and Selection
S. noctilio samples were collected from native European, non-Chinese invaders, and Chinese populations, comprising 44, 127, and 607 specimens, respectively. A summary of sampling sites and common parameters of the genetic diversity of S. noctilio in this study is provided in Table 1. GenBank/DRYAD provided 23 of the 44 native European specimens; others were dry or anhydrous ethanol-preserved individuals borrowed from museums or donated by local institutions. Non-Chinese invader populations were represented by 22, 39, 50 and 16 GenBank/DRYAD samples from Africa, North America, South America, and Oceania, respectively. All Genbank/DRYAD samples are retrieved from DRYAD: https://datadryad.org/stash/dataset/doi:10.5061/dryad.37mm8; from Genebank: https: //www.ncbi.nlm.nih.gov/nuccore/?term=sirex+noctilio. Finally, we obtained a total of 607 specimens from 19 Chinese populations (Table 1 and Figure 1), donated by various local forestry bureaus or obtained from the standard quarantine facility of Beijing Forestry University, which contains a large number of locally collected infested stems; these specimens were identified using taxonomic literature on Sirex woodwasps before DNA extraction [26]. As the Chinese source areas were infested at different times and to different degrees, the number of specimens and sampling year vary. Dried museum samples were handled outside our laboratory, as described below; specimens already in 99.7% pure anhydrous alcohol remained so; all remaining specimens were frozen at −80 • C.
Insects 2020, 11, x FOR PEER REVIEW 3 of 14 the first time 778 individuals from 6 continents. In particular, samples from original European populations and recent Chinese infestations were included. Overall, this study aims to elucidate the invasive routes of this pest as well as provide data enabling the prevention of the harmful consequences of S. noctilio spread.

Sample Collection and Selection
S. noctilio samples were collected from native European, non-Chinese invaders, and Chinese populations, comprising 44, 127, and 607 specimens, respectively. A summary of sampling sites and common parameters of the genetic diversity of S. noctilio in this study is provided in Table 1. GenBank/DRYAD provided 23 of the 44 native European specimens; others were dry or anhydrous ethanol-preserved individuals borrowed from museums or donated by local institutions. Non-Chinese invader populations were represented by 22, 39, 50 and 16 GenBank/DRYAD samples from Africa, North America, South America, and Oceania, respectively. All Genbank/DRYAD samples are retrieved from DRYAD: https://datadryad.org/stash/dataset/doi:10.5061/dryad.37mm8; from Genebank: https://www.ncbi.nlm.nih.gov/nuccore/?term=sirex+noctilio. Finally, we obtained a total of 607 specimens from 19 Chinese populations (Table 1 and Figure 1), donated by various local forestry bureaus or obtained from the standard quarantine facility of Beijing Forestry University, which contains a large number of locally collected infested stems; these specimens were identified using taxonomic literature on Sirex woodwasps before DNA extraction [26]. As the Chinese source areas were infested at different times and to different degrees, the number of specimens and sampling year vary. Dried museum samples were handled outside our laboratory, as described below; specimens already in 99.7% pure anhydrous alcohol remained so; all remaining specimens were frozen at −80 °C.      As there is an obvious sample size bias towards the six Chinese populations with high sample numbers (population codes DM, QQ, HG, FJ, JBT, CT), we used only 20 of the sequenced individuals, which included all the haplotypes found in each population. We selected individuals from each haplotype based on the selection ratio and the number of samples in each haplotype. The final total of 364 samples for analysis consisted of 44, 127, and 193 individuals from Europe, non-Chinese invaded areas, and China, respectively.

DNA Sequence Analysis
For the native European samples not supplied as sequences directly from GenBank, DNA extraction, PCR amplification and sequencing of the barcode fragment were performed for 15 specimens at the Laboratory of Forest Zoology URZF, INRAE (Orléans, Paris, France). DNA extracts were prepared from one adult hind leg or thorax muscle using NucleoSpin tissue XS kits (Macherey Nagel, Dylan, Germany) according to the manufacturer's protocol. The COI barcoding fragment was amplified via PCR using the forward primer LCO1490 (5 -GGT CAA CAA ATC ATA AAG ATA TTG G-3') and the reverse primer HC02198, (5'-TAA ACT TCA GGG TGA CCA AAA AAT CA-3') (Vrijenhoek, Rijswijk, Netherlands, 1994). PCR reagents were used and reactions performed following Sun [24]. PCR products were purified using the QIAquick PCR Purification Kit (Qiagen, Hilden, Germany). Sequencing was performed using the Sanger method with an ABI Prism BigDye Terminator v3.1 cycle sequencing kit (Thermo-Fisher, Waltham, MA, USA). The reagents and procedure used were: 1 µL 5 × bigdye buffer, 1 µL bigdye, 1 µL LCO/HCO, diluted to 3.3 µM, 2 µL DNA; 30 cycles of 10 s at 96 • C, 5 s at 50 • C, 3 min at 60 • C. The remaining 6 specimens were borrowed from museums and had to remain intact; for these the same COI fragment was amplified at the Canadian Centre for DNA Barcoding (CCDB-Biodiversity Institute of Ontario, University of Guelph) using a slightly different sequencing set (C-LepFolF/C-LepFolR), following the standard high-throughput protocol [27].
DNA was extracted from the leg or thorax muscle of each Chinese specimen using the Multisource Genomic DNA Miniprep Kit (Axygen, San Francisco, CA, USA), following the manufacturer's protocol to prevent cross-contamination, all tools were sterilized by flame or 75% pure anhydrous ethanol. DNA was eluted using 70 µL of elution buffer and stored at −20 • C until use. Species-specific COI amplification primers were synthesized by a commercial company (SinoGenoMax, Beijing, China). PCR primers LCO1490 and HC02198 were used as above. Reactions were performed in a total volume of 25 µL, containing 12.5 µL of 2 × GoTaq ® Green Master Mix (Promega, Madison, WI, USA), 1 µL of each primer, 1.5 µL of DNA template, and 9 µL of ultrapure water. The reaction procedure was as follows: initial denaturation for 2 min at 94 • C, followed by 35 cycles of denaturation at 94 • C for 30 s, annealing at 45 • C for 30 s, and extension at 68 • C for 2 min, with a final extension at 68 • C for 10 min. PCR products (3 µL) were analyzed by electrophoresis on a 1.5% (w/v) agarose gel (1 × TAE), alongside a DNA marker (D2000, Takara, Kyoto, Japan). After electrophoresis at 130 V for 20 min, the PCR products were visualized using ethidium bromide and ultraviolet light. Finally, products clearly visible after electrophoresis were sent to a commercial company (SinoGenoMax, China) for sequencing in both directions.
All electropherograms were checked manually in CodonCode Aligner v3.7.1. (CodonCode Corporation, Centerville, MA, USA) to assess their quality. The alignment was done independently using Clustal W implemented in MEGA v6.06 with default parameters, and no pseudogenes or stop codons were detected in the sequences [28,29]. All COI sequences were trimmed to the same length (322 bp) for final alignment.
Details on the collecting data for each specimen, as well as photographs, sequence records, trace files, and primer sequences used for PCR amplification, together with GenBank accession numbers, are available through the dataset in BOLD (www.boldsystems.org).

Statistical Analysis
Common parameters of genetic diversity, including the number of haplotypes, haplotype, and nucleotide diversity were calculated using DNAsp v6 [30]. The haplotype network was built in PopART [31][32][33]. Haplotype distribution and frequency were projected on a map using ArcGIS 9.3 (ESRI, Redlands, CA, USA).
Analyses of molecular variance (AMOVA) were performed on native European populations and Chinese populations to measure the partitioning of genetic variation between populations and groups of populations in Arlequin 3.5.2.2 [34,35]. The first AMOVA was performed on four clusters of Chinese populations in the northeast of China defined by geographic proximity: western (mountain area), middle (city area), which was divided into two parts, and eastern (city and mountain areas). A second AMOVA was performed on two clusters of Chinese populations grouped by latitude, north and south, in order to test the trend of spread.

Results
There were 11 variable sites among all insects sampled. Three dominant haplotypes, H1, H2, and H9, were present among the total of 12 haplotypes; the haplotype diversity was 0.4207. Only one or two mutation sites existed across the different haplotypes ( Figure 2). There were two mutation sites between H1 and H12. H1 and H9 were haplotypes that were made up of populations from different countries: H1 was widely located on all continents, whereas H9 was not found in China and Africa (Figures 2 and 3). H2, a major Chinese haplotype along with H1, was found only in Chinese samples ( Figure 2). H3, H4, H5, H6, H7, and H8 were unique to Chinese populations. H8, the haplotype found near the national border with Russia ( Figure 4) is worth noting. H10 was unique to South America, as were H11 and H12 to Europe.
Insects 2020, 11, x FOR PEER REVIEW 7 of 14 [31][32][33]. Haplotype distribution and frequency were projected on a map using ArcGIS 9.3 (ESRI, Redlands, CA, USA). Analyses of molecular variance (AMOVA) were performed on native European populations and Chinese populations to measure the partitioning of genetic variation between populations and groups of populations in Arlequin 3.5.2.2 [34,35]. The first AMOVA was performed on four clusters of Chinese populations in the northeast of China defined by geographic proximity: western (mountain area), middle (city area), which was divided into two parts, and eastern (city and mountain areas). A second AMOVA was performed on two clusters of Chinese populations grouped by latitude, north and south, in order to test the trend of spread.

Results
There were 11 variable sites among all insects sampled. Three dominant haplotypes, H1, H2, and H9, were present among the total of 12 haplotypes; the haplotype diversity was 0.4207. Only one or two mutation sites existed across the different haplotypes ( Figure 2). There were two mutation sites between H1 and H12. H1 and H9 were haplotypes that were made up of populations from different countries: H1 was widely located on all continents, whereas H9 was not found in China and Africa (Figures 2 and 3). H2, a major Chinese haplotype along with H1, was found only in Chinese samples ( Figure 2). H3, H4, H5, H6, H7, and H8 were unique to Chinese populations. H8, the haplotype found near the national border with Russia ( Figure 4) is worth noting. H10 was unique to South America, as were H11 and H12 to Europe.      AMOVA analysis of the genetic structure of Chinese populations (Table 2) showed that when grouped by geographic proximity (Figure 5), most of the variation (38.64%, p < 0.01) was between geographic clusters, with a relatively small part of the genetic variance (9.15%, p < 0.01) between populations within geographic clusters. When sites were grouped by latitude (Figure 6), the variation (25.58%, p < 0.01) between clusters was slightly higher than that (22.46%, p < 0.01) between populations within clusters. Within this study, no clear genetic structure emerged within native European populations.  AMOVA analysis of the genetic structure of Chinese populations (Table 2) showed that when grouped by geographic proximity (Figure 5), most of the variation (38.64%, p < 0.01) was between geographic clusters, with a relatively small part of the genetic variance (9.15%, p < 0.01) between populations within geographic clusters. When sites were grouped by latitude (Figure 6), the variation (25.58%, p < 0.01) between clusters was slightly higher than that (22.46%, p < 0.01) between populations within clusters. Within this study, no clear genetic structure emerged within native European populations.   The Chinese sampling sites are divided into four groups based on geographical proximity: west for one group, middle for two groups and east for one group. The middle and eastern groups, where there are cities, contained 4 haplotypes, whereas the western group, which is mainly mountainous, had 3 haplotypes. The first site discovered in China, DM, contained 40% H1 and a similar proportion, 55%, of H2. There is no H1 in the western group, although the proportion of H1 in the two middle area groups is quite high as is also found in the eastern group. There is a decreasing trend from west to east in the proportion of H2. The representative populations on this axis were HDQ, DM, and HG with 100%, 55%, and 20% H2, respectively. Among the other six unique haplotypes in China, H8 is located in MZL, on the border between China and Russia. Four of the remaining five haplotypes are distributed in the populations from which large samples were taken and are located in city areas. The sixth, HHEJ, is located in the area from which P. sylvestris var. mongolica seedlings are taken for distribution and planting elsewhere. H3 is present in 5% of DM individuals, and so does H4, H5, H6, H7, these haplotypes occurred in 10%, 5%, 7%, and 25% of sampled insects of HG, FJ, YS and HHEJ, respectively. The latitude 46.204 • N was used to group samples into the categories north and south. Haplotype diversity in the northern area is obviously higher than that in the south. The northern cluster has seven haplotypes, H1, H2, H3, H4, H5, H7, and H8, from nine populations, compared to two haplotypes (H1 and H6) from ten populations in southern China.

Host Unity of S. Noctilio in China
The host range of S. noctilio is rather large, although it is known that the pest is harmful to various kinds of conifers around the world, especially Pinus, Picea, Abies, Larix, and Pseudotsuga [19]. Both natural forests and pine plantations can suffer damage, especially when they are overstocked and stressed. In China, the origin of P. sylvestris var. mongolica is in the area of Honghuaerji in Inner Mongolia and north of the Huma River in Heilongjiang [36]. Field investigation found no S. noctilio in the wild distribution areas of P. sylvestris var. mongolica in either mountain or sandy ecological zones (Figure 1), other than two larvae found in the western, sandy area. This occurrence, and the unique haplotype, probably came from the considerable imports of Russian timber required for the construction of Honghuaerji Scenic Area.
Currently, the main infestation zones in China are those where P. sylvestris var. mongolica has been introduced. Two typical examples are Hegang (Heilongjiang) and Jinbaotun (Inner Mongolia). The former is a mixed forest of Pinus koraiensis, L. olgensis, Picea koraiensis and P. sylvestris var. mongolica. The latter contains Pinus tabuliformis and P. sylvestris var. mongolica. However, fieldwork has indicated that the only infested conifer species in China to date is P. sylvestris var. mongolica, even when other potential hosts are present in the same forest [37]. As is characteristic of invasive pests, S. noctilio appears to be harmless in the natural distribution area of its favorite host tree species.

Patterns of Multiple Invasion and Spread in China
There are currently three more frequent haplotypes, H1, H2, and H9, which have a worldwide distribution. It is thus unsurprising that they are in the center of the haplotype network and closely connected with other well-established haplotypes [38]. Interestingly, H1, despite being globally the most prevalent haplotype, is absent in western China, although strongly represented in the city regions of central and eastern China. This constitutes evidence that genetic communication is human activity-derived in both invasion and spread scenarios. Given that the western component is the primary seedling source of P. sylvestris var. mongolica, trees from this area will be planted elsewhere. We speculate that H1 is an invasion haplotype carried by humans that have entered the central and eastern urban areas of China. As it existed on all continents, its source is unclear. Based on our data, H2, H7 and H8 were found only in China, and given that they are unlikely to have originated there, we are unable to identify their origin. H2 is the main haplotype of the western component, with a distribution ratio that gradually decreases from west to east; this is in line with the characteristics of spread from a point of origin to other regions. This, as noted above, may represent an invasion originating with Russian wood imported for construction at Honghuaerji and spread eastward by planting. Severely endangered woodlands in Dumeng, Hegang, Fujin, and other eastern areas may have been invaded in this way. H7 and H8 are also unique to the western region; H8 occurs on the border with Russia.
Grouped by latitude, significantly more haplotypes exist in northern Chinese S. noctilio populations than in those to the south. Damage to woodlands in the northeastern region occurs discretely, in multiple forest farms, parks, and other wooded land. This damage was concentrated, as is characteristic of the middle to late stages of an invasion. In contrast, in newly-invaded southern areas, such as Jinbaotun, damaged forests were localized and individual affected trees were distributed discretely in the affected area, indicating the initial stage of introduction. Although the 1130 km between Honghuaerji and Jinbaotun is shorter than the 1400 km between Honghuaerji and the affected area in the northeast, as the haplotype H1 found in the southern region is not present in Honghuaerji, the south regional populations may originate from northern China.
S. noctilio has a relatively strong natural flight ability. In the southern hemisphere, the maximum spread rate of an invading population of S. noctilio is 82 km per year [39]. In China, the diffusion rate of the pest is 74.57 km/year [40]. The straight-line distances between Dumeng, the first site detected, and the other two severely endangered areas, Hegang and Jinbaotun, are 440 km and 390 km, respectively, which would take about 6 years to achieve at 74.57 km/year. Yet the woodwasp was already present in 2014 and 2015. The first recorded occurrence sites in all the affected areas are distributed along with the main traffic networks (Figure 1), which thus appear to have increased the impact caused by human activities, such as wood packaging and log transportation.

Focus for Further Research
When analyzing biological invasions, knowledge of the genetic structure of populations from the area of origin is vital in determining the source. Our target species is native to Eurasia and North Africa, but it does not cause serious harm there and has consequently received less attention, leading to difficulties in collecting samples. In particular, the lack of samples from the Russian region caused sampling breakpoints between the European and northeast Chinese populations. The observation that H2 occurred only in China may be due to the large sample size from some previously unsampled sites from this area. More reliable results may be obtained from sampling efforts focused on this region and on western China. Finally, microsatellite marker use involving cytochrome C oxidase subunit I (COI) could prove effective, showing its utility as a method for determining invasion pathways and spread.

Conclusions
In this study, a 322 bp fragment of the mitochondrial barcode gene COI was used to analysis the genetic diversity and structure of S. noctilio populations in both native and invaded ranges, with a specific focus in China. Twelve haplotypes were found across the native and invaded distribution of the pest, of which three were dominant; among these there were only one or two mutational steps between each pair of haplotypes. No obvious genetic structure was found other than in Chinese populations. China has a unique and dominant haplotype not found elsewhere, and compared with the rest of the world, the genetic structure of Chinese populations suggested a multiple invasion scenario.
Funding: This research was funded by Beijing's Science and Technology Planning Project, grant number Z191100008519004.