Genetic and morphological variation of Vespa velutina nigrithorax which is an invasive species in a mountainous area

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 coefficient 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.

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, approaches to control the species are still lacking and need to be studied continuously 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 responded to elevation gradients in 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 highaltitude 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. This phenomenon can lead to isolated populations and low genetic diversity, which may destabilize the population and act as a mechanism to control the distribution or maintenance of the population [10][11][12][13] . 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 flight performance. If flight is vital for colony maintenance, and wing loading increases at

Genetic distances and population structures.
The Mantel test showed no significant 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 STRU CTU RE analysis revealed that the population genetic structures of the yellow-legged hornets www.nature.com/scientificreports/ 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 yellowlegged 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 significantly 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 significantly larger (Tukey test, p < 0.05) than wing size in G10. There were no significant differences in wing size among the other groups. Canonical variate analysis (CVA) results revealed two major morphological variations among the hornets from 13 sampling sites based on the two major axes CV1 and CV2 (Fig. 5a, b). CV1 explained 32.11% of wing morphological variance, and CV2 explained 20.20% of wing morphological variance. Changes in CV1 represent major variations in landmarks 3, 13, 8, and 9. As CV1 increases, the wing margin shape slightly changes. Changes in CV2 represent the major variations in landmarks 3, 4, 9, 10, 13, and 16. The shape of the wing slightly enlarges as CV2 increases. The wing shapes of hornets at each site were separated by CV1 and CV2 (Fig. 5c). CV1 (F = 45.22, p < 0.0001) and CV2 (F = 28.44, p < 0.0001) were significantly different among yellow-legged hornets from the 13 groups. Hornets in groups G3, G5, G8, and G9 were distinct (Tukey test, p < 0.05) from all the other groups based on CV1 (Fig. 5d). Based on CV2 (Tukey's test, p < 0.05), hornets of groups at seven sites (G2, G4, G6, G8, G9, G10, G13) were separated from those of other groups (Fig. 5e).
The morphological distance is expressed based on Procrustes distance from CVA (Table 2). Unlike direct comparison of CV1 axis and CV2 axis, most of group did not have a significant morphological distance overall. Significant morphological distances were found only in some groups. Wing shape in G7 was significantly distant (p < 0.05) from that in G8. Wing shape in G10 was significantly distant (p < 0.05) from that in G11 and wing shape in G13 was significantly distant (p < 0.05) from those in G7 and G8.

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. 15 . The yellow-legged hornet can be observed at such high altitudes because its colony foraging radius is greater than 700 m 16 . Furthermore, in autumn, when foraging activity is the greatest, workers are found at sites relatively far from their nests 17 . Although snow cover and cool temperatures at high altitudes could limit www.nature.com/scientificreports/ 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 significantly with an increase in elevation. Yellow-legged hornet only exhibited a high inbreeding coefficient at a specific site. The high level of inbreeding is consistent with the genetic bottleneck experienced by hornets following their introduction in South Korea 18 . 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 yellow-legged hornet invaded from a single or several sites 19 . In our study, although many hornets and nests were observed at relatively lower altitudes, their genetic structure did not show altitude-specific patterns. The yellow-legged hornet on Mt. Jiri appeared as a single 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 intraspecific minor 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 20 . A decrease in temperature due to an increase in altitude does not seem to be sufficient to alter wing shape. However, our study might be limited due to insufficient number of samples at high altitudes. We conducted permutation tests to resolve this limitation. Nevertheless, no altitude-specific pattern of morphological variation was observed. Our results demonstrated www.nature.com/scientificreports/ genetic and morphological adaptations specific to altitude in the hornets of Mt. Jiri; however, an extension of these results may be dangerous.
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 .  www.nature.com/scientificreports/ However, V. simillima, which is less aggressive than the yellow-legged 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 21 . 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.

Material and methods
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. At a lower altitude (approximately 200-800 m) and a higher altitude (approximately 800-1200 m), a total of 20 and 3 hornet nests were identified, respectively. 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. Same individuals were used in both genetic and morphological analysis, for tracking response of individual units in two types of variations. 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 Scientific, Wilmington, USA), and the concentrations of genomic DNA were diluted between 10 and 20 ng/µL range. Seven polymorphic microsatellite loci were amplified using primers developed by 22 for the yellow-legged hornet (D2-185, D3-15, R1-137, R1-36, R4-33, R1-75, and R1-169). We performed amplifications using the procedures for each locus detailed in 22 . 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 23 was used to confirm deviations in Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium of seven microsatellite loci. Null alleles, significant deviations from HWE, or evidence of linkage disequilibrium were not observed in any of the seven loci. In addition, seven loci could sufficiently separate 39 samples of yellow-legged hornets (Fig. 7). Table 2. Morphological distances of yellow-legged hornets from 13 sampling sites. Right value represents paired p values of Procrustes distances that were obtained using canonical variate analysis (CVA) after 10,000 permutation rounds among yellow-legged hornets from 13 sampling sites, and left value represents paired Procrustes distances from CVA. Significantly different distances are indicated by gray boxes. Bayesian clustering was performed using STRU CTU RE version 2.3.4 28 . STRU CTU RE is an admixture model, which allows us to confirm 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 29 in STRU CTU RE Harvester 30 , which supported the estimation of the best-identified k value. The range of one to 14 possible clusters with three independent runs each was employed in STRU CTU RE Harvester. In the STRU CTU RE harvest analysis results, the best-supported K value among hornets across 13 populations was identified as 7 (Fig. 8).
DAPC 31 , 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 32 . 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  Morphological analysis. Landmark-based geometric morphometrics were used to analyze the wing shapes of the hornets. TpsDig 34 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 35 , 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 first two (CV1 and CV2) axes. We also identified the values and significant distances of the Procrustes from CVA to explore the morphological similarity between groups. Significant morphological distances were obtained using CVA after 10,000 permutation rounds. Oneway ANOVA was used to test significant differences in centroid size, CV1, and CV2, among the hornets using GraphPad Prism version 7.0 for Windows (GraphPad Software, San Diego, USA). When significant differences were found in the ANOVA test, Tukey's post hoc tests were performed. All differences were considered significant at p < 0.05.   www.nature.com/scientificreports/