Genetic and Morphological Variation of Vespa Velutina Nigrithorax, an Invasive Species in Mountainous Areas

The yellow-legged hornet (Vespa velutina nigrithorax) is an invasive species in South Korea with negative economic, ecological, and public health impacts. We investigated genetic and morphological variation in the species populations on Mt. Jiri, the tallest mountain in South Korea. We hypothesized that a high-altitude would be negatively correlated with the genetic diversity of the hornet population, and hornet wing morphology would change with an increase in altitude. Our results showed that the genetic diversity of yellow-legged hornets did not decrease as altitude increased. Regardless of the altitude, the inbreeding coecient was high at the newly colonized sites. A single genetic population occurred in the mountainous areas examined and gradually expanded its range. Wing morphology, especially shape, did not change with an increase in altitude or decrease in temperature. Although snow cover and cool temperatures at high altitudes could limit nest-building activities, they did not prevent the extension of the range of the species. Therefore, the yellow-legged hornet cannot be controlled naturally by climate or topography; combined approaches, including chemical control, nest removal, and bait-trapping techniques should be implemented.


Introduction
The yellow-legged hornet (Vespa velutina nigrithorax), native to China, was introduced in South Korea in 2003 [1] and it has rapidly colonized a large part of the country [2]. The spread of V. velutina nigrithorax has had negative economic, ecological, and public health effects [3][4][5]. However, no approaches to control the species have been identi ed or adopted [6].
High-altitude mountainous areas could serve as biogeographical barriers, as high altitudes may be negatively correlated with the occurrence of the yellow-legged hornet [7,8]. In particular, low-temperature conditions and snow cover in mountainous areas prevent the overwintering of mated queens. Nevertheless, the hornet is still found at high altitudes in mountain sites in South Korea. The country has a high proportion of mountainous areas when compared with other countries where the yellow-legged hornet successfully colonized, such as France, the UK, and Spain. Human-mediated transportation may explain the nest distribution at considerable distances from the invasion front; however, this dispersal may not be responsible for range expansion at high-altitude sites on mountains [9]. Nests of the yellow-legged hornet are more frequently observed at the edges of the mountains than at high-altitude sites; however, the hornet tends to gradually expand its range toward mountainous areas [4]. Investigating whether the yellow-legged hornet is genetically and morphologically adapted to new mountainous ranges could facilitate the establishment of effective control strategies for this invasive species.
In this study, we examined genetic and morphological variations in yellow-legged hornet populations that have successfully colonized in the mountainous area in South Korea. We hypothesized that hornets at high-altitude sites have less genetic variation than hornets that have colonized low-altitude sites because those populations would have colonized the areas relatively recently, and the low temperatures and snow cover may interrupt overwintering and primary nesting after overwintering. In addition, we hypothesized that hornets that have colonized high-altitude sites have distinct wing morphology. Cold temperatures lower the wing beating frequency and thus power output, resulting in reduced ight performance. If ight is vital for colony maintenance, and wing loading increases at low temperatures, hornets could reduce wing loading by increasing wing size [10].

Genetic distances and population structures
The Mantel test showed no signi cant genetic isolation (Mantel statistic r = -0.006, p = 0.538) based on geographical distances among yellow-legged hornets from the 13 sampling sites (Fig. 1). According to the constructed unweighted pair-group method with arithmetic (UPGMA) tree, yellow-legged hornets from G13 were clearly separated from the other groups (Fig. 2a). However, the STRUCTURE analysis revealed that the population genetic structures of the yellow-legged hornets were not different among the 13 sampling sites. The 39 yellow-legged hornets from all sites were grouped into the same K-cluster (Fig. 2b). In Discriminant analysis of the principal components (DAPC), discriminant function 1 (DF1) explained 61.30%, and discriminant function 2 (DF2) explained 12.38% of the total genetic variation in the hornets. DF 1 indicated that yellow-legged hornets in G8, G11, and G12 were separated from the other groups, and the rest were grouped into three clusters (G1, G9, G10/G6, G7, G13/G2, G3, G4, and G5). Conversely, DF2 indicated that the hornets associated with G2 were isolated, whereas all the other groups were grouped into successive clusters (Fig. 3a). Finer-scale structures detected through DAPC separated the yellow-legged hornets from the 13 sites into 13 distinct groups (Fig. 3b).

Morphological differences in wing shape
Centroid size, which is a proxy for wing size, was signi cantly different (one-way analysis of variance (ANOVA), F=2.304, p = 0.036) among the hornets from the 13 sampling sites (Fig. 4). Wing size in G4 was signi cantly larger (Tukey test, p < 0.05) than wing size in G10. There were no signi cant differences in wing size among the other groups.
The morphological distance is expressed based on Procrustes distance ( Table 2). Wing shape in G7 was signi cantly distant (p < 0.05) from wing shape in G8. In addition, wing shape in G10 was signi cantly distant (p < 0.05) from wing shape in G11, whereas wing shape in G13 was signi cantly distant (p < 0.05) from wing shape in G7, G10, and G12.

Discussion
The yellow-legged hornet is expanding its range in mountainous areas, regardless of altitude [2]. In the present study, adult hornets were collected at altitudes as high as 1100 m a.s.l. in mountainous areas. In Italy, adult hornets have been found at altitudes of 1200 m a.s.l, however, most nests of the yellow-legged hornet have been reported within 700 m a.s.l. [11]. The yellow-legged hornet can be observed at such high altitudes because its colony foraging radius is greater than 700 m [12]. Furthermore, in autumn, when foraging activity is the greatest, workers are found at sites relatively far from their nests [13]. Although snow cover and cool temperatures at high altitudes could limit hornet nest-building activities, Mt. Jiri, the highest mountain in South Korea, has not prevented the extension of the range of the yellow-legged hornet.
In this study, yellow-legged hornet genetic diversity did not decrease signi cantly with an increase in elevation. Yellow-legged hornet only exhibited a high inbreeding coe cient at a speci c site. The high level of inbreeding is consistent with the genetic bottleneck experienced by hornets following their introduction in South Korea [14]. The genetic bottlenecks were not drastic enough to limit the expansion of the ranges of the yellow-legged hornet into mountainous areas. In addition, the hornets are likely to experience temporary genetic bottlenecks at the local level when expanding their colonies annually. In South Korea, there is a debate on whether the yellowlegged hornet invaded from a single or several sites [15]. If the yellow-legged hornet was introduced from various sites and dispersed in many directions, populations with different genetic structures should be observed on Mt. Jiri, which is in the southern center of the Korean Peninsula. The yellow-legged hornet on Mt. Jiri is a single genetic population. The single genetic population is spreading at a slightly slower rate than in other countries, without being restricted by mountainous areas as geographical barriers.
Yellow-legged hornets that expand their ranges across or into mountainous areas may have no morphological adaptations. Although intraspeci c variability in wing morphology may exist in the population investigated in the present study, wing morphology was not correlated with altitude. Additionally, yellow-legged hornets prefer cooler highlands and mountains in their original habitats [16]. A decrease in temperature due to an increase in altitude does not seem to be su cient to alter wing shape.
In conclusion, South Korea does not have a climate or topography that can naturally prevent the spread of wasps. In the early stages of the invasion of the yellow-legged hornet, the spread and establishment of the hornet in mountainous habitats was hampered by competition with an ecologically similar hornet, Vespa simillima [4]. However, V. simillima, which is less aggressive than the yellowlegged hornet, could have been forced out of the mountainous area, and the hornet occupied its place mountainous areas in South Korea [4]. Typically, nest removal and bait-trapping are used to control the yellow-legged hornet. However, eradicating and controlling hornets that occupy large mountainous areas using such techniques is challenging and impractical [17]. Therefore, chemical control techniques, including the use of insect growth regulators that have not been considered so far, could be adopted in mountainous areas, which is a core habitat of the yellow-legged hornet in South Korea.

Sampling sites
We collected yellow-legged hornets using traps at Mt. Jiri National Park in September 2019 when the colony was at its maximum size and the production of new queens resulted in high demand for hornet workers. Sugared water with acetic acid and 95% ethanol was used as the trap solution in 2 L plastic containers. The traps were monitored every two weeks. The collected hornets were preserved in 95% ethanol for use in later molecular and morphological analyses. Traps were placed in 30 sites and hornets were captured from 13 sites (Fig. 6, Table 3). We selected three individuals from each site, totalizing 39 similar-sized workers.

DNA extraction and microsatellite genotyping
We collected muscle samples from the hind legs of each worker for genomic DNA extraction. DNeasy Blood and Tissue kits (Qiagen, Hilden, Germany) were used to extract genomic DNA according to the manufacturer's protocol. We measured the concentrations and quality of genomic DNA using a NanoDrop 2000 spectrophotometer (Thermo Scienti c, Wilmington, USA), and the concentrations of genomic DNA were diluted between 10 to 20 ng/µL range. Seven polymorphic microsatellite loci were ampli ed using primers developed by [18] for the yellow-legged hornet (D2-185, D3-15, R1-137, R1-36, R4-33, R1-75, and R1-169). We performed ampli cations using the procedures for each locus detailed in [18]. Seq-Studio Genetic Analyzer (Applied Biosystems, Foster City, USA) was used to visualize the amplicons, and GeneMapper version 6.0 was used to evaluate the dataset for genotype errors and the presence of null alleles.

Analysis of genetic diversity and distance, and population structure
The genotype datasets were used to analyze genetic diversity and population structure. GENEPOP version 4.7 [19] was used to con rm deviations in Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium of seven microsatellite loci. Null alleles, signi cant deviations from HWE, or evidence of linkage disequilibrium were not observed in any of the seven loci. In addition, seven loci could su ciently separate 39 samples of yellow-legged hornets (Fig. 7).
The MS Excel add-in, GenAlEx version 6.5 [20], and Arlequin version 3.5 [21], were used to calculate genetic diversity and diversity indices, namely mean number of alleles (N A ), effective number of alleles (N E ), observed heterozygosity (H O ), expected heterozygosity (H E ), Shannon's information index (I), molecular diversity (h), and inbreeding coe cient relative to the sub-population (F IS ) of each population. The UPGMA based on Euclidean distances was used to construct a dendrogram based on Nei's genetic distances obtained using GenAlEx version 6.5, using PAST 3 [22]. Paired population differentiation (F ST ) was performed using Arlequin version 3.5. The Bayesian clustering was performed using STRUCTURE version 2.3.4 [24]. An admixture model, which allows us to con rm whether ancestors in population k have passed a portion of their genetic material to individual i. Simulations (100,000) were performed in each analysis after an initial burn-in of 100,000 simulations. We used the ΔK method [25] in STRUCTURE Harvester [26], which supported the estimation of the best-identi ed k value. The range of one to 14 possible clusters with three independent runs each was employed in STRUCTURE Harvester. In the STRUCTURE harvest analysis results, the best-supported K value among hornets across 13 populations was identi ed as 7 (Fig. 8).
DAPC [27], which is a multivariate clustering algorithm, was used to analyze the population structures of 39 hornets from 13 sites using the "adegenet" package in R [28]. The analysis describes the greatest amount of variation through a discriminant function representing a linear combination of correlated alleles in a linear discriminant analysis based on principal components generated after reducing the dimension of the genetic variation using principal component analysis [29]. The method was used to determine how distinct each population among the 13 sites was.

Morphological analysis
Landmark-based geometric morphometrics were used to analyze the wing shapes of the hornets. TpsDig [30] was used to digitize the 19 landmark points in the wing vein (Fig. 9). Morpho J version 1.07a (Manchester, UK) was used to convert the digitized landmark coordinates into Procrustes coordinates. Centroid size, which used the size proxy in landmark morphometrics [31], was measured to compare wing sizes among the 13 groups. CVA was used to compare morphological differences among the 13 groups and to calculate the f contribution (%) of each canonical variate (CV) in CVA. Variation in wing shape among the hornets from 13 sites was visualized in a wireframe graph along the rst two (CV1 and CV2) axes. We also identi ed the values and signi cant distances of the Procrustes to explore the morphological similarity between groups. One-way ANOVA was used to test signi cant differences in centroid size, CV1, and CV2, among the hornets using GraphPad Prism version 7.0 for Windows (GraphPad Software, San Diego, USA). When signi cant differences were found in the ANOVA test, Tukey's post hoc tests were performed. All differences were considered signi cant at p < 0.05. Tables   Table 1. Genetic diversity and diversity indices of yellow-legged hornets (Vespa velutina nigrithorax) based on seven nuclear microsatellite loci estimated in 13 sampling sites: mean number of alleles (N A ), number of effective alleles (N E ), observed (H O ) and expected (H E ) heterozygosity, Shannon's information index (I), inbreeding coe cient relative to the population [29], and molecular diversity (h). Correlation between the linearized pairwise FST and geographical distance for the 13 yellow-legged hornet sampling sites in Mt. Jiri. The Mantel test revealed no signi cant genetic isolation based on geographical distances among 13 sampling sites.

Figure 2
Results of Nei's genetic distance and Bayesian clustering analysis in 39 yellow-legged hornets from 13 sites: (a) unweighted pair group method with arithmetic mean (Euclidean distance) tree among yellow-legged hornets from 13 sampling sites; (b) population clusters obtained by STRUCTURE analysis in yellow-legged hornets from 13 sampling sites.   Lowercase letter and box color indicate differences (p < 0.05) in centroid size based on Tukey's post-hoc tests.

Figure 5
Page 13/15 Canonical variate analysis (CVA) and wireframe graph with 19 landmarks illustrating morphological variance in wing shape among yellow-legged hornet from 13 sites: (a) wireframe graph illustrating morphological differences based on canonical variate 1 (CV1   Number of multilocus genotypes based on number of loci from 39 individual of yellow-legged hornet. When the number of loci was 6 or more, the 39 individuals could be separated about 100%.

Figure 8
Optimal k-means for STRUCTURE analysis obtained by the ΔK (delta K) method in STRUCTURE harvest. Here, 100,000 simulations after burn-in of 100,000 simulations were used to obtain the three independent runs for range of 2 to 14 possible clusters. The highest delta K value was 7 (3.102).

Figure 9
Locations of the landmark coordinates in the right forewings of a yellow-legged hornet. Nineteen landmarks were designated in the wing vein.