Images of forewings specimens of V. velutina representing 13 different color morphs were obtained from specimens housed at 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), and 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 (Perrard et al. 2014a,b).
We selected specimens to broadly sample the native distributions of the morphs and the type localities of named taxa. 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 (Fig. 1). Our full sample thus included the following named taxa: V. v. var. velutina Lepeletier, 1836; V. v. var. auraria Smith, 1852; V. v. var. ardens du Buysson, 1905 (1904); V. auraria var. nigrithorax du Buysson, 1905 (1904); V. v. var. celebensis Pérez, 1910; V. mongolica var. divergens Pérez, 1910; V. flavitarsus Sonan, 1929; V. auraria flavitarsis Ma, 1937; V. v. variana van der Vecht, 1957 v. karnyi van der Vecht, 1957 v. sumbana van der Vecht, 1957 v. timorensis van der Vecht, 1957 v. floresiana van der Vecht, 1957d v. pruthii van der Vecht, 1959.
We partitioned our data into three groupings. The first included all the specimens from the morphs listed above (N = 368). This data set allowed us to test whether sympatric morphs are distinct in their wing venation patterns. The second data partition was restricted to specimens from non-sympatric locations where only one morph is found (N = 238; see list of synonyms and type localities in Smith-Pardo et al. 2020). This reduced data set included 13 morphs: V.v. ardens: Lombox and Sumbawa: Indonesia; V.v. auraria: northern India; V.v. celebensis: Sulawesi: Indonesia; V.v. divergens: Malaysia; V.v. flavitarsus: Taiwan; V.v. floresiana: Flores: Indonesia; V. v. karnyi: Sumatra: Indonesia; V. v. nigrithorax: Guangdong Province: China; V. v. pruthii: northern India; V. v. sumbana: Sumba: Indonesia; V. v. timorensis: Timor: Indonesia; V. v. variana: Thailand; V. v. velutina: Java: Indonesia, plus the intercepted individual from Utah. This data set allowed us to focus on whether wing venation patterns can be linked to a distinct geographic origin. Finally, for some analyses, we analyzed the mean shapes of these latter samples (N = 14). This data set allowed us to assess patterns of phenotypic and biogeographic similarity among the morphs. A list of all specimens is provided in Supplemental File S1.
Wings of museum specimens were imaged after relaxing them with water vapor for several days and extending the forewings with translucid Millard strips over a gel foam surface. Images were taken with a Leica S9i dissecting microscope with a digital camera and the LAX software (Leica). Images were rotated and trimmed using Adobe Photoshop (Adobe CS) and later open in the program tpsDig2w64 for landmarking (Rohlf, 2006).
Nineteen landmarks for each forewing were marked on each image following the scheme developed by Perrard et al. 2014a (Fig. 2). The coordinates of the landmarks were recorded in TPS filed format for further analysis. In the TPS files (Supplemental File S2), specimen metadata is encoded in the 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 genus Vespa, species velutina, subspecies ardens, from Indonesia, is from the material described by Perrard et al (2014a) and is the 16th individual in that series.
We performed a series of geometric morphometric analyses (Bookstein 1991, Dryden and Mardia 1998) to clean data, explore the general pattern of wing shape variation, assess the distinctness of wing venation between color morphs, assess patterns of similarity among the morphs, and to classify the Utah specimen. We used three software packages to perform our analyses: MorphoJ v. 1.07a (Klingenberg 2011), PAST v. 4.04 (Hammer et al. 2001), and Geometric Morphometrics for Mathematica v. 12.4 (Polly 2022).
Each analysis began with Procrustes superimposition to standardize the orientation and size of the landmark coordinates (Rohlf and Slice 1990). To circumvent analytical complications caused by covariances among the landmarks and the loss of degrees of freedom from Procrustes superimposition, landmarks were transformed to orthogonal shape variables by projecting them onto their principal components using the covariance method and retaining the axes with non-zero variance (Dryden and Mardia 1998, Polly and Motz 2016). For the first two data sets this transformation produced 34 shape variables (2k-4, where k is the number of landmarks) and 13 for the third data set, which is sample limited (n-1, where n is the number of shapes in the data set) (Rohlf and Slice 1990, Zelditch et al. 2004).
Principal components analysis (Pearson 1901) was used to explore the variation of wing venation and to inspect data for outliers. Outlier inspection was also facilitated by comparing the cumulative distribution of squared Mahalanobis distance between each specimen and the sample mean using MorphoJ’s outlier inspection tool (Klingenberg 2011).
Multivariate analysis of variance (MANOVA) (Wishart 1928) with a non-parametric permutation test was used 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 could be found on any axis, even those that account for only a small proportion of the total variance (Uyeda et al. 2015, Polly and Motz 2016, Adams and Collyer 2018). Pairwise post-hoc tests for significant differences between each pair of morphs with Bonferroni-style correction for multiple comparisons were also carried out.
Discriminant function analysis (DFA), also known as linear discriminant analysis (LDA) or canonical variates analysis (CVA), was used to find the components of shape that best separate the morphs and to classify unidentified specimens (Rao 1948). This method finds the set of axes that best separate groups, similar in its goal to between-groups principal components analysis but standardized. Both methods are prone to over-fitting problems when the number of variables exceeds the number of groups (Mitteroecker and Bookstein 2011, Cardini et al. 2019, Cardini and Polly 2019) 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 based on the wing venation discriminant functions.
Continuous trait maximum-likelihood (Felsenstein 1973, 2004) was used to build a tree to show the hierarchical pattern of similarity in wing venation. 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 (Caumul and Polly 2005). The ML tree was projected into principal component space (Rohlf 2002) 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 (Felsenstein 1989).