Identifying Morphs of the Yellow-Legged Hornet (Vespa velutina) and Other Pests of Quarantine Importance with Geometric Morphometrics

Abstract

We assess the accuracy of geometric morphometrics (GMM) for determining the origin of insects of quarantine importance using the Asian hornet (Vespa velutina Lep.1836) as a case study. This species is highly variable, has an extensive natural distribution, and has been transported to many regions of the world. Forewing landmarks were applied to a large sample of regionally specific color morphs (previously considered “subspecies”) from across the species’ native Asian range. We reconfirm that GMM can statistically distinguish geographic variants independent of the color patterns that have heretofore been used for provenance, but which have been suspected of being unreliable. Almost all morphs in our analyses were statistically different except the centrally located V. v. variana, whose range lies between the continental V. v. auraria Smith, 1852, and V. v. nigrithorax du Buysson, 1905 morphs, and the Malaysian and Indonesian morphs. Even with moderate-sized training samples, discriminant function analysis (DFA) was able to classify geographic morphos with about 90% accuracy (ranging from 60% to 100%). We apply these results to determine the origin of a dead wasp recently intercepted in a mail parcel in Utah. Both DFA and continuous-trait maximum-likelihood clustering suggest that the Utah specimen belongs to the nigrithorax morph, which is native to southern China but now invasive in Europe, Japan, and Korea. These results are also supported by DNA barcode analysis, which groups the Utah individual with nigrithorax populations in South Korea and Japan. The relationship between variation in wing shape and genetic differentiation deserves further study, but molecular data are consistent with the GMM results suggesting that morphometric comparisons may be able to identify and provenance intercepted specimens quickly and inexpensively when molecular sequences and taxonomic specialists are unavailable.

1. Introduction

The Asian hornet (Vespa velutina), also known as the Yellow-legged hornet, is an invasive species of social wasp that is now spreading rapidly around the world but was originally native to southeast Asia, where is it distributed from Afghanistan to East Timor in Indonesia [1,2] (Figure 1). Vespa velutina is known to prey on domestic honeybees, Apis mellifera L., and their stings can pose serious health risks to people who are allergic [3,4,5]. Individuals live in colonies of several hundred to a few thousand individuals. Their nests are made of pulp (cellulose) and are built above ground on the branches of trees. Their diet consists of arthropods, particularly insects, which are used as protein sources in addition to nectar from flowers, extra-floral nectaries, honeydew, ripe fruits, sap flows, and occasionally man-made sweet substances [6].

This species is one of the most invasive hornets known. In Asia, it was introduced into South Korea (around 2003) where it has become invasive [8,9], and it was reported as naturalized on Tsushima Island between 2010 and 2012 [10] and around 2015 on Kyushu Island [11], both of which are islands in Japan. In Europe, the Asian hornet was first introduced into France on commodities from China and quickly spread to Spain (2010), Portugal (2011), Italy (2012), Great Britain (2016), Belgium (2017), and Luxembourg (2020) at a speed of several hundred miles per year [3,12,13,14,15,16]. Vespa velutina has also been introduced into some Mediterranean Islands, including the Balearic Islands of Spain [17]. It has now reached as far north as Frankfurt, Germany, a distribution record that also represents the northernmost distribution for the species and which also shows its invasive potential [18].

While there were more than 50 interceptions of hornets representing at least six species at US ports of entry from 2010 to 2018 (see review [2]), the Asian hornet had not been intercepted in North America until recently. In early 2020, a dead hornet was found by a concerned citizen in Salt Lake City, Utah, inside an e-commerce package. The specimen was later identified as V. velutina, probably nigrithorax, by one of us (AHS-P) and was later confirmed as V. velutina sensu stricto by entomologists at the Systematic Entomology Laboratory (USDA-SEL) in Washington D.C. The package that contained the dead specimen was shipped from Germany, where Asian hornets are now established, but its ultimate origin is not clear, as the merchandise in the package seems to have been manufactured somewhere else. The interception of this Asian hornet, even in the form of a dead specimen, was significant for three reasons: (1) it was the first confirmed interception of V. velutina in the United States; (2) the specimen was first detected in Utah, not at one of the U.S. ports of entry, and (3) the species is one of the most invasive of the true hornets due to its behavioral plasticity and broad diet breadth.

Determining the geographic origin of intercepted Asian hornets is not a simple matter. Vespa velutina is variable in its coloration patterns to the extent that at least 17 subspecies or other named taxa have been described for the regional variants found across its southeastern Asian distribution [7,19,20,21]. Its coloration ranges from almost entirely yellow or orange to extensively black [20,22] and can be variable even within populations [7]. A recent investigation using mitochondrial DNA and microsatellite loci suggested that color patterns do not map cleanly onto genetic structuring within the species [7]. A genetic analysis identified two divergent groups: the first consisting of mainland wasps from China, Southeast Asia, Nepal, Kashmir, and northern Pakistan; the other consisting of wasps from the islands of Indonesia. Unfortunately, color patterns are shared in parallel or convergently in these groups, suggesting that color by itself may not be a good guide to geographic provenance. But because this species is rapidly expanding its range by hitchhiking in packages and containers that are shipped around the world, an inexpensive and accurate method for determining the taxonomic identity and provenance of intercepted Asian hornets is desirable. Officials in agricultural inspection agencies often have few tools, limited taxonomic expertise, and large budgets for identifying intercepted insects, including the Asian hornet. Geometric morphometric analysis of wing venation patterns has been long-established as an aid to the taxonomic identification of flying insects, e.g., [23,24,25,26,27] and recent work has shown that it has considerable potential for discriminating among the color morphs of Asian hornets [1,28].

Asian hornets have long been classified by color pattern, but the taxonomic status of color morphs is questioned [29] and coloration tends to grade from one region to another, calling into question its efficacy for identifying the ultimate natural source populations of invasive wasps [2,30]. Many competing taxonomic criteria like features of the genitalia or molecular sequencing require specialized knowledge or expensive analytical tests. However, forewing shape analysis, which can be carried out with unspecialized photographic and computational equipment, has promise for rapid identification of invasive individuals by local customs or wildlife inspection teams. Perrard and colleagues [1,7] used geometric morphometric analysis to show that forewing shape varies quantitatively among Vespa species and geographic sub-populations of Vespa velutina. They found that a naive clustering of wing shape from taxonomically and geographically heterogeneous samples of workers (females) returned groups congruent with species and confirmed that individuals can be reliably attributed to their genus, species, and populations based on their wing shape. They also confirmed that variation in coloration between populations was unrelated to molecular, geographic, or climatic differences and that it may be the product of convergent evolution that is as much influenced by constraints on aposematism and Müllerian mimicry as by abiotic pressures on melanism, a conclusion further supported by Do and colleagues [31], who found that genetic diversity and forewing shape in invasive populations in a mountainous area of South Korea were invariant with altitude or temperature. This work suggests that forewing shape may be a useful and inexpensive method for determining the provenance of intercepted invasive wasps.

In this paper, we evaluate the effectiveness of geometric morphometric methods (GMM) for distinguishing between geographic variants of Vespa velutina and classifying intercepted specimens using wing venation. Our objectives are as follows: (1) to assemble a database of forewing landmarks for all the morphological variants, or morphs (subspecies or races), of V. velutina from its entire geographic distribution (native and expanded by introduction); (2) to further evaluate the accuracy of GMM analysis of forewing venation patterns for classifying specimens from different source locations using an expanded data set; (3) based on this, to determine the race and likely origin of the specimen intercepted in the e-commerce package in Utah; and (4) to apply DNA barcoding to assess the accuracy of the GMM identification.

2. Materials and Methods

2.1. Geometric Morphometrics

Images of forewings specimens of V. velutina were taken by the senior author from specimens housed at the California State Collection of Arthropods at the California Department of Food and Agriculture in Sacramento (CDFA), the American Museum of Natural History, New York (AMNH, James Carpenter and Christine LeBeau), the Bohart Museum of Entomology at University of California, Davis (BME, Lynn Kimsey, and Steve Heydon). Other photographs were obtained from Dr. Adrien Perrard at the Institute of Ecology and Environmental Sciences IEES—Paris (PERR) from some of his work on hornets, including the material used in Perrard’s previous studies of genetic and color variation of this species [1,7]. Those photographs are available upon request from Dr. Perrard. The wings of museum specimens were imaged after being relaxed with water vapor for several days and the forewings were extended over a gel foam surface with 4 mil/0.1 mm clear mylar strips. Images for the borrowed specimens were taken with a Leica S9i dissecting microscope with a digital camera and the LAX software (Leica 1.4.6). All images were rotated and trimmed using Adobe Photoshop (Adobe CS4).

We selected specimens to broadly sample the native distributions of the morphs in the region of their type localities. Because of the uncertainty of how color morphs relate to genetic differentiation and geographic location, we used provenanced specimens and distinguished between morphs that have completely endemic distributions (such as in Taiwan, Malaysia, and some islands of Indonesia) from ones that coexist in broad ranges of sympatry, particularly in continental Vietnam, China, Laos, Myanmar, India, and Nepal (Figure 1). This strategy avoids ambiguity related to debates about hybridization, sympatry, and polymorphism in those areas where more than one morph has been identified and focuses on patterns of phenotypic and biogeographic similarity among unambiguous morphs. Because the photographs were taken from museum specimens collected by others in the past, we have no direct information about the location of the local hornet colonies but because the specimens were taken by different collectors on different dates, often separated by decades, we assume that each individual is from a different colony.

Our data set consisted of 233 individuals belonging to 11 morphs (see the list of synonyms and type localities in [2]): V. v. ardens: Lombox and Sumbawa: Indonesia (N = 11); V. v. auraria: northern India (N = 18); V. v. celebensis Pérez, 1910: Sulawesi: Indonesia (N = 6); V. v. divergens Pérez, 1910: Malaysia (N = 34); V. v. floresiana van der Vecht, 1957: Flores: Indonesia (N = 17); V. v. karnyi van der Vecht, 1957: Sumatra: Indonesia (N = 24); V. v. nigrithorax: Guangdong Province: China (N = 32); V. v. sumbana van der Vecht, 1957: Sumba: Indonesia (N = 25); V. v. timorensis van der Vecht, 1957: Timor: Indonesia (N = 8); V. v. variana van der Vecht, 1957: Thailand (N = 8); V. v. velutina Lepeletier, 1836: Java: Indonesia (N = 49), plus the intercepted individual from Utah (N = 1). This data set allowed us to focus on whether wing venation patterns can be linked to a distinct geographic origin. Finally, some analyses were conducted based on the mean shapes of these latter samples (N = 12). A list of all specimens is provided in Table S1.

Nineteen landmarks for each forewing were marked on each image using the program tpsDig2w64 [32] following the scheme developed by Perrard and colleagues [1] (Figure 2). All landmark coordinates were digitized by us (AHS-P). In the TPS files (S2 Appendix), specimen metadata is encoded in an ID tag of each specimen: morph name, geographic location, collection, and specimen ID. For example, “ID = ardens-Indo-Perr-16” indicates the specimen belongs to the morph ardens from Indonesia with a photograph by Perrard and colleagues, and it is the 16th individual from that source. We checked for outliers and other digitization errors using visual inspection and by comparing the cumulative distribution of the squared Mahalanobis distance between each specimen and the sample mean using MorphoJ’s outlier inspection tool [33,34,35]. Based on the observed shape variation (Figure S1) and the multivariate normal model (Figure S2), we excluded specimens with flaws in the wings and re-digitized those with landmarking errors to arrive at our N = 233.

We used three software packages to perform our analyses: MorphoJ v. 1.07a for outlier analysis [35], PAST v. 4.04 for DFA analysis [36], and Geometric Morphometrics for Mathematica v. 12.4 for the rest of the analyses [37].

Each analysis began with Procrustes superimposition to standardize the orientation and size of the landmark coordinates [38]. To circumvent analytical complications caused by covariances among the landmarks and the loss of degrees of freedom from Procrustes superimposition, landmarks were transformed into orthogonal shape variables by projecting them onto their principal components using the covariance method and retaining the axes with non-zero variance [34,39] For the first two data sets, this transformation produced 34 shape variables (2k − 4, where k is the number of landmarks), and 13 were produced for the third data set, which is sample limited (n − 1, where n is the number of shapes in the data set) [38,40].

Principal components analysis [41] was used to explore the variation of wing venation in V. velutina and to determine whether the Utah specimen fell within the phenotypic range of any of these pure-morph samples from their native ranges.

To determine and confirm whether the morphs themselves were statistically different from each other, we used multivariate analysis of variance (MANOVA) [42] with a non-parametric permutation test to test differences in shape among the groups for statistical significance. These tests used all the shape variables (PC axes) with non-zero variance, because when shape differences are small and many groups are being compared, small but meaningful differences can be found on any axis, even those that account for only a small proportion of the total variance [39,43,44]. Pairwise post hoc tests for significant differences between each pair of morphs with Bonferroni-style correction for multiple comparisons were also carried out.

To assign the intercepted Utah specimen to a morph based on its wing venation, we used discriminant function analysis (DFA), also known as linear discriminant analysis (LDA) or canonical variates analysis (CVA), to find the components of shape that best separate the morphs and to classify the intercepted specimen [45]. DFA finds the set of axes that best separate groups. It has a similar goal to between-groups principal components analysis, but the variances are standardized, and the axes are those that best separate the groups. Both methods are prone to over-fitting problems when the number of variables exceeds the number of groups [46,47,48] so only the first twelve principal components (one fewer than the number of groups) were used. Leave-one-out cross-validation was used to measure the accuracy of classification.

DFA always classifies an unknown to one of the training groups regardless of whether it is a good fit, so we also used continuous trait maximum-likelihood [49,50] to build a tree to show the hierarchical pattern of similarity in wing venation among the morphs and to determine which morph the Utah specimen is closest to. To construct the tree, the sample was Procrustes superimposed, a consensus shape was estimated for each of the 13 morphs plus the intercepted specimen from Utah, their means were superimposed again, and shape variables were obtained through PCA. The shape variables were used as traits in the continuous trait ML procedure. To retain undistorted shape information, the variables were not standardized. This method of constructing a morphometric tree is preferred over other commonly used methods like UPGMA or minimum spanning trees because it does not distort pairwise distances by imposing an ultrametric topology, it treats each shape variable as an independent trait, and it adopts an explicit evolutionary model, Brownian motion in this case [51]. The ML tree was projected into principal component space [52] to illustrate the multivariate pattern of similarity in reduced-dimension PC plots. The CONTML module of PHYLIP 3.695 was used to construct the tree [53].

2.2. Molecular Sequencing and Analysis

To test whether our morphometric identifications were correct, we used DNA barcoding to match the intercepted specimen from Utah to existing barcode databases for V. velutina. DNA was extracted using a Lucigen MasterPure DNA extraction kit (Lucigen Corp., Middleton, WI, USA). The whole dry specimen was removed from the pin, placed in a 50 mL conical centrifuge tube, immersed in approximately 500 mL of Tissue and Cell Lysis Solution and 1 μL of Proteinase K, and heated to 65 °C overnight in a water bath. After overnight incubation, the specimen was removed from the lysis buffer, rinsed with absolute ethanol, and returned to the pin. The remaining extraction was carried out according to the manufacturer’s instructions with an additional incubation at −20 °C to increase DNA yield, which took place after the addition of isopropanol, and with final elution in Buffer EB (QIAGEN Sciences, Germantown, MD, USA) to ensure downstream compatibility for use in PCR.

PCR reactions were performed with TaKaRa Ex Taq HS polymerase (Takara Bio, Shiga, Japan) in total volumes of 50 μL using the manufacturer’s recommended volumes of 10X Ex Taq buffer and dNTP mixture. The primers HCO/LCO [54,55] were used to amplify a 658 bp segment of cytochrome c oxidase I (COI) on a Bio-Rad C1000 Touch (Bio-Rad Laboratories, Inc., Hercules, CA, USA). PCR conditions included an initial denaturation step of 94 °C (3 min), 39 cycles of 94 °C (20 s)/50 °C (20 s)/72 °C (30 s), and an extension step of 72 °C (5 min). Amplicons were purified using a Qiaquick PCR Purification Kit and eluted into 35 μL of EB buffer. Sanger sequencing was performed by GENEWIZ (Azenta Life Sciences, South Plainfield, NJ, USA). Individual forward and reverse contigs were assembled using Geneious Prime 2021 (Biomatters Ltd., Auckland, New Zealand), manually trimmed, and examined for errors.

An additional 97 publicly available V. velutina sequences were downloaded from the Barcode of Life Data System website (BOLD) [56]. These were combined with the sequence of the Utah specimen and three additional Vespa/Vespula sequences that were used as outgroups. All sequences (101 total) were aligned with MAFFT ver. 6 using the G-INS-i algorithm [57]. To examine sequence similarity, a neighbor-joining (distance) tree was constructed under the Kimura 2 parameter model (K2P) for nucleotide substitutions using PAUP* [58]. Provenance data for the BOLD database specimens are reported in Table S3.

3. Results

The mean forewing shapes for each of our regional samples are illustrated in Figure 3. The most pronounced differences are between samples whose wings are broader anteroposteriorly (e.g., ardens, sumbana, velutina) and those that are narrower or more compressed in that direction (e.g., nigrithorax, variana). There are also visible differences in the shape and relative proportions of the submarginal cells. For example, the average nigrithorax wing has second and third submarginal cells that are subequal in area, with the third marginal having parallel medial and lateral bars of approximately equal length, while the average celebensis has a proportionally larger third marginal cell with bars that are not parallel. The specimen intercepted in Utah has an anteroposteriorly compressed wing with second and third submarginal cells like nigrithorax.

A principal component morphospace of all specimens is shown in Figure 4. The shape differences represented by each of the first four principal components (PC) axes are illustrated in Figure 5. PC 1 separates the anteroposteriorly narrow wings (negative end) from the broad-winged morphs. PCs 2–4 describe differences in specific wing cells. PC 2 primarily relates to differences in the marginal and second submarginal cells, PC 3 to the mediolateral proportion of the first submarginal cell, and PC 4 to the anteroposterior breadth of the cell lying between the median and subcubital veins. The morphospace demonstrates that some color morphs are very distinct from one another, but there is considerable overlap between morphs on any given axis of the shape space. For example, timorensis, a broad-winged morph, is completely non-overlap** divergens, a narrow-winged morph, on PC 1 (Figure 4A,D), but divergens overlap with nigrithorax to at least some degree on each of the first four PC axes. The morphospace has 34 dimensions with non-zero variance (19 landmarks × 2 dimensions − 4 dimensions lost to translation in each dimension, scaling, and rotation in the Procrustes alignment) and its first four axes account for 62.6% of the total variance with PC1 = 28.8%, PC2 = 14.7%, PC3 = 12.2%, and PC4 = 7.0% (Table S2, Figure S6).

Despite the overlap in shape variation, the MANOVA test confirmed that the majority of the regional samples are statistically distinct from each other in the mean shape of their wings (

Table 1. MANOVA p-values from pairwise post hoc comparisons between morphs. Gray values are not statistically significant.

								SS	df	MS
auraria	0.00						Between	0.023	10	0.0023
celebensis	0.09	0.00					Within	0.035	221	0.0002
divergens	0.00	0.00	0.00				Total	0.058	231	0.0003
floresiana	0.00	0.00	0.01	0.00			F-ratio = 14.59
karnyi	0.00	0.00	0.00	0.00	0.00		P(non-parametric) < 0.001
nigrithorax	0.00	0.01	0.00	0.00	0.00	0.00		R2 = 0.40
sumbana	0.00	0.00	0.01	0.00	0.00	0.00	0.00
timorensis	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
variana	0.00	0.13	0.02	0.00	0.02	0.00	0.13	0.01	0.00
velutina	0.00	0.00	0.01	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	ardens	auraria	celebensis	divergens	floresiana	karnyi	nigrithorax	sumbana	timorensis	variana

). The other pairs that were not significantly different from each other were ardens and celebensis, variana and auraria, and variana and nigrithorax, all of which are each other’s geographic nearest neighbors (Figure 1).

To determine which aspects of wing shape best separate the regional samples and thus would best classify the intercepted specimen from Utah, we performed a canonical variates analysis (CVA; Figure 6). Each canonical axis separates two or more groups from one another. The loadings of landmark variables on those axes can be combined into a discriminant functions analysis (DFA) that can then be used to find the best classification of an unknown specimen to one of the groups in the analysis. Because there are 13 morphs, there are as many as twelve meaningful canonical axes. The distribution of groups in the first two is shown in Figure 6. Axis 1 separates the narrow-winged morphs (e.g., nigrithorax) from the broad-winged ones (e.g., ardens, velutina, timorensis, celebensis). Axis 2 separates divergens from timorensis, which further differ in the proportional size of the marginal cell and the shapes of the second and third submarginal cells.

The CVA plot strongly affiliates the Utah specimen with the narrow-winged morphs and the associated DFA classified it as nigrithorax.

We also assessed the morphometric affinity of the Utah specimen by constructing a continuous-trait maximum likelihood tree using the wing shape variables for it and the mean shapes of the 13 groups. The intercepted specimen is closest to the nigrithorax, as seen in both the tree itself (Figure 7A) and the connections among the taxa in the morphospace plot (Figure 7B).

We confirmed the affinities of the Utah specimen using DNA barcoding using the cytochrome c oxidase subunit I (COI) region (Figure 8). Although data for more than 370 individuals of V. velutina are now present in the BOLD database, sequences from many locations are 100% identical, so we did not include dozens of identical sequences (primarily from South Korea) for this analysis. Many specimens in BOLD were identified as subspecies (color morph) by the original authors, eight of which are represented (color-coded in Figure 8 to match Figure 1; precise locations shown in Figure S7; Table S3), but some were not (indicated with black type). Like the findings from Perrard and colleagues [1,7], the molecular tree does not unambiguously cluster the specimens by subspecies (or color morph), indicating that COI barcodes do not fully differentiate between morphs. As such, putative subspecies labels were removed, and sequences are listed by country in Figure 8. The specimen from Utah clusters with a large group from South Korea and Japan (and one individual from the vetulina morph from Java). This group is separated from another large group containing specimens from China, India, and Southeast Asia, and from another large group consisting of invasive nigrithorax specimens from Europe and unprovenanced individuals from China. Together, these three groups contain the majority of specimens putatively identified as nigrithorax. The DNA barcode data are therefore consistent with the morphometric estimation of the origin of the Utah specimen but are more precise in aligning it specifically with the invasive populations of Korea and Japan (for which no morphometric data were available).

4. Discussion

Our results reconfirm the conclusion of Perrard and colleagues [7] that statistically significant geographic variation in forewing shape exists within the natural range of V. vetulina and extend that work by demonstrating that accurate determinations of the geographic origin of an individual wasp can be made using forewing shape. These results suggest that forewing shape analysis can be applied to real-world interception problems quickly and inexpensively without relying on assumptions of the geographic significance of color patterns, specialist knowledge needed to classify individuals using genitalia characters, or the expense and specialized training required to genotype them with molecular sequencing. By using a forewing shape to immediately identify an intercepted wasp, rapid quarantine decisions can be made at ports of entry.

4.1. Accuracy of GMM Identifications

Using DFA, GMM was able to classify specimens to geographically specific Asian hornet color morphs based on forewing shape with a fairly high degree of accuracy. The accuracy of the DFA classification was assessed with leave-one-out cross-validation, a procedure that iteratively recalculates the CVA by drop** a different specimen each time. Each dropped specimen is treated as an unknown and classified based on the shape data.

In our analysis, 207 out of 232 specimens were correctly classified, which is an overall accuracy of 89%, indicating that considerable confidence can be placed in the morphometric classification. The ability of DFA to accurately classify some morphs is greater than others, as shown in the “confusion matrix” from the cross-validation (

Table 2. Discriminant function confusion matrix with tallies of classifications. Numbers in rows are correct (gray) and incorrect (white) classifications; numbers in columns are true identities of specimens classified as a morph. The far-right column reports how frequently a morph was classified correctly; the bottom row reports how often a DFA classification was correct.

	ardens	auraria	celebensis	divergens	floresiana	karnyi	nigrithorax	sumbana	timorensis	varians	velutina	Total	% Correctly Classified
ardens	11	0	0	0	0	0	0	0	0	0	0	11	100%
auraria	0	13	0	0	0	1	1	1	0	2	0	18	72%
celebensis	0	0	5	0	0	1	0	0	0	0	0	6	83%
divergens	0	0	0	33	0	0	0	1	0	0	0	34	97%
floresiana	0	2	0	0	15	0	0	0	0	0	0	17	88%
karnyi	0	1	0	0	1	22	0	0	0	0	0	24	92%
nigrithorax	0	1	0	0	0	2	27	0	0	2	0	32	84%
sumbana	0	0	0	0	0	0	0	25	0	0	0	25	100%
timorensis	0	0	0	0	1	0	0	0	7	0	0	8	88%
varians	0	1	0	0	1	0	1	0	0	5	0	8	63%
velutina	2	0	1	0	2	0	0	0	0	0	44	49	90%
Total	13	18	6	33	20	26	29	27	7	9	44	232
% classifications correct	85%	72%	83%	100%	75%	85%	93%	93%	100%	56%	100%

). For each morph (rows), the table reports the number of individuals classified as one morph or another (columns). For example, all 11 ardens specimens were correctly classified as ardens, indicating that this morph is very likely to be classified and identified. However, there were also two velutina specimens classified as ardens, which means that there is an 11 out of 13 chance (85%) that a DFA classification of ardens is correct. Thus, the percentages in the right column represent the number of times a morph was classified correctly, and the bottom row represents the number of times any given classification is correct.

We found that celebensis, divergens, and timorensis could be classified with virtually complete confidence. Most other classifications had confidence of 80% or higher.

While the DNA barcode data confirm the GMM identification of the Utah specimen as nigrithorax, it is interesting to note that the topology of the GMM tree is only partially congruent with the barcode tree, suggesting that the former may not be good at reconstructing phylogeny. For example, in the barcode tree, divergens from southern peninsular Malaysia and velutina from Java lie close to one another relative to other morphs (Figure 8), but in the morphometric M-L tree, these morphs are distant from each other as well as in the PCA morphospace and (Figure 4). On the other hand, similarities between some morphs, as seen in the barcode tree, are recovered in the morphometric tree. Assessing the ability of wing venation patterns to recover phylogenetic topology is beyond the scope of this paper, and readers should exercise caution in drawing evolutionary conclusions from the tree in Figure 7.

4.2. Strengths and Weaknesses of the DNA Barcode Analysis

Not surprisingly, our DNA analysis provided similar results to past studies using much of the same data [7,59]. Although DNA barcodes are frequently used to separate similar taxa and “cryptic species”, e.g., [60], traditional DNA barcoding in Hymenoptera can be thwarted by issues that limit its usefulness, such as mtDNA introgression between taxa facilitated by Wolbachia infection [61]. Indeed, bee-specific primers have been designed to prevent the amplification of Wolbachia genes [61]. It is not known if similar issues are affecting Vespa barcodes, but success varies depending on the specific taxon. In general, DNA barcodes can separate species of Vespa, e.g., [59], but their utility at the subspecific level is less certain. While morphs or subspecies of V. mandarinia can be reliably separated with barcode data, other methods are necessary for species like V. velutina. Expanding beyond basic COI DNA barcodes is one solution, but geographic sampling is difficult and the overall number of taxa and sequences in databases such as BOLD is hard to match. This study demonstrates that augmenting DNA sequence data with GMM could be a way to avoid the shortfalls of both methods and that DNA results may be helpfully augmented with morphological considerations.

4.3. The Provenance of the Specimen Intercepted in Utah

DFA classified the Utah specimen as nigrithorax based on its forewing shape, a group for which 93% of classifications are correct. One auraria and one variana specimen were also incorrectly classified as nigrithorax, but no other morphs were inappropriately classified this way. The DFA results suggest that the Utah specimen is most likely to be nigrithorax but has a low probability of being an auraria or variana. This identification is consistent with P-values associated with Mahalanobis squared distances of the Utah specimen from each of the morphs that form the basis of the DFA classification (

Table 3. Mahalanobis distances squared (D2) and probabilities (P) that the Utah specimen belongs to each morph based on the first five dimensions of the shape space. A significant chance of group membership is indicated in boldface.

Taxon	D2	P (same)
ardens	62.8	<0.0001
auraria	35.2	<0.0001
celebensis	892.1	<0.0001
divergens	66.4	<0.0001
floresiana	79.7	<0.0001
karnyi	28.1	<0.0001
nigrithorax	12.7	0.01
sumbana	65.26	<0.0001
timorensis	696.9	<0.0001
variana	44.4	<0.0001
velutina	24.7	<0.0001

). The Utah specimen lies outside all groups except nigrithorax and pruthii, with p < 0.0001, and has a 0.01 probability of belonging to the latter two. As discussed above, pruthii has N = 2, so the covariance matrix needed for the Mahalanobis distance is unreliably estimated, but with N = 29 the nigrithorax P value is more reliable. Note, however, that all of the P values in Table 3 assume that the group is multinormally distributed, which is not necessarily the case and is not assumed for the other tests presented in this paper. That nigrithorax is the best morphometric identification of the Utah specimen is also supported by the M-L tree of forewing shape similarity, which also links the Utah specimen to nigrithorax. Note that this tree does not provide a phylogenetic analysis but rather an analysis of similarity based on the geometry of the wings. Although it is not possible to definitively separate velutina subspecies or color morphs using COI barcode data [7], our barcode data confirm this identification by clustering the Utah haplotype closely with the invasive nigrithorax specimens from South Korea and Japan, as well as more distantly with the invasive nigrithorax in Europe and the native individuals in China (Figure 8).

The Utah specimen is somewhat of an outlier on the GMM tree, as indicated by its long branch length and distance from the other morphs in the morphospace plot (Figure 7). Its divergence is at least partly due to it being represented only by a singular specimen, since the morphs are represented by the means of their respective samples (compare Figure 7B with Figure 4A, where the full range of variation for each taxon is shown), but some of its divergence may also be due to it being derived from invasive populations in Korea or Japan, which themselves might be expected to be diverged from their source populations in China due to founder genetic effects and subsequent population size expansion, c.f., [62,63]. The long branch length leading to the Utah forewing shape in Figure 7 might, in principle, be filled in by the forewing shape from the invasives in Korea and Japan, a proposition that can be tested by future studies.

The GMM and DNA results collectively suggest that the Utah hornet most likely originated from nigrithorax populations that are now established in Korea and Japan. By themselves, the GMM results cannot distinguish sources in the native range in China and Southeast Asia from those in Europe, Korea, Japan, or elsewhere where nigrithorax has become established because we do not have forewing shape data from those populations. Furthermore, the DFA results also suggest a small probability that the Utah specimen belongs to the auraria or variana morph, which would indicate an origin somewhere on the Southeast Asian mainland. Regardless, the GMM results rule out many other possible origins. DNA barcode analysis confirms the GMM identification of the Utah specimen as nigrithorax by placing it in a cluster of invasive nigrithorax specimens from South Korea and Japan identified (along with one individual from Indonesia identified as subspecies velutina). Similar to other recent studies, e.g., [13,18,64,65,66], we found that haplotype diversity within areas invaded by velutina is relatively low due to a single introduction of only one or a few individuals. As such, it is possible to determine that the Utah specimen likely did not originate in Europe from where its package was shipped.

The invasive populations in Europe are themselves derived from at least one founder queen belonging to the nigrithorax group that arrived in France in 2004 in trade products from China and then spread into Germany [13,67,68]. Whether the intercepted hornet originated in Germany or at a manufacturing point in SE Asia remains ambiguous because our morphometric results can be deemed consistent with either option. Mitochondrial DNA suggests that the invasive populations in Europe originated in Zhejiang, Jiangsu, or Shanghai provinces in China [13], which are natively inhabited by the nigrithorax morph as it is recognized by us.

Our nigrithorax GMM sample is from Guandong province (Hong Kong) a few hundred kilometers to the south. In morphometric terms, the specimen from Utah best matches this Guandong sample, while the DNA sequence results indicate that it originates from a cluster that includes Korea and Japan, which may explain why it appears to be a partial outlier from the main nigrithorax cluster in our morphometric analysis (Figure 4C and Figure 7). Whether the difference is due to the Utah specimen having originated from a European invasive population, which is no doubt distinguished by founder effects from its parent population in China, or whether it is sampled from a Korean (based on the molecular data results) or Chinese population that differs from Guandong by regional variation is currently unknown.

4.4. Using GMM to Determine the Provenance of Invasive Insects

Our study adds to a growing set of tools that will help agriculture inspection agents identify invasive species. GMM has been used to classify specimens of unknown taxonomic affinity for at least 20 years, e.g., [25,69,70,71,72]. The ability of these methods to accurately identify specimens varies with the taxonomic group, body part, degree of divergence among the candidate taxa, and sampling, e.g., [51,73,74]. Despite its potential, GMM has only been applied to the identification of intercepted quarantine species or assessed for accuracy in that context in a few cases (e.g., [75] on snails, [76] on lady beetles). Our results suggest that the standard use of wing landmarks to assess shape differences of different populations has the potential to be used to make quick and inexpensive preliminary identifications to exclude potentially invasive species at ports of entry where molecular sequencing is impractical and where taxonomic experts are unavailable. If a more extensive database for forewing landmarks is combined with automated landmarking techniques, e.g., [26,77,78,79,80], machine learning—particularly convolutional neural networks—based on the use of extensive datasets of images of the different species, e.g., [81,82,83] a project we are currently working on, and phone apps, e.g., [84], it is feasible to provide instant identification of Asian hornets at the point of interception.

5. Conclusions

These results are promising for controlling the continued spread of the Asian hornet because they suggest that it may be possible to easily and inexpensively determine the ultimate origin of individuals that are transported with goods if we know their population differences (in terms of the shape of the forewings) and their distributions. Knowing the origin of an invasive pest is important to establish control or safeguarding plans at the origin and prevent their dissemination to other areas; in addition, in cases like V. velutina, knowing the origin of the invading population (race or subspecies) can provide valuable information about the invasive pest’s behavior, phenology, and potential expansion range. The Asian hornet has one of the largest distribution ranges of any species of Vespa, and has the most morphological diversity within the genus, making it difficult or impossible for local inspection officials to make the fine-scale taxonomic determinations that would indicate where an intercepted individual originated. Because V. velutina has been introduced in more countries outside its natural distribution than any other species of the genus, this GMM tool is needed. Such tools are also valuable when other data sources, such as commonly used DNA barcodes, are insufficient for separating populations or even taxa due to overall invariability in the barcode region or the founder effect for new introductions.

Our results from geometric morphometrics suggest that the hornet intercepted in Utah belongs to the nigrithorax group, but they also indicate that the individual is an outlier relative to the training sample of nigrithorax from its home range in China. The DNA barcode analysis confirmed that the Utah individual has a nearly identical sequence to nigrithorax individuals, but it groups it with invasive populations in Korea, a location from which we have no morphometric data. The barcode thus serves to confirm that the morphometric data are both capable of assigning a V. vetulina species to its progenitor morph as well as identifying potential discrepancies that indicate that the precise source is not represented in the training data. Because our training data were intentionally limited to the native ranges of morphs, where they live in isolation from other morphs, the GMM data by themselves are incapable of resolving whether the Utah specimen originated in the native range of nigrithorax in China, from invasive populations in Europe (including Germany, from where the intercepted package was shipped), or from invasive populations in Korea and Japan. Even though the DNA barcode data seem to indicate that the Utah specimen originated from Japan or Korea and not from Germany or China, to determine this morphometrically would require training data from established invasive populations around the world. Ultimately, expanding the Vespa wing shape database to include both native and invasive populations of all taxa and extending geographic sampling over the native range of these species will help further resolve these ambiguities and improve our ability to accurately and inexpensively provenance intercepted specimens in situations where more sophisticated solutions like molecular sequencing are not feasible in the time an agent has to make a quarantine decision.