Phylogeographic Analysis to Understand the Origin, Speciation and Biogeographic Expansion of the Invasive Asian Hornet, Vespa velutina

1. Introduction

The genus Vespa, natively distributed in tropical, subtropical, and temperate Asia, features several invasive species such as V. velutina [1,2,3] and potentially invasive species such as V. mandarinia [4,5,6,7]. Several hornet species are mass murderers of honeybees [8,9], and the diversity and distribution of hornets closely mirror those of honeybees. This dependence of the hornets on honeybees for food led to an interesting evolutionary twist. An orchid species, Dendrobium sinense, has evolved a function of emitting a honeybee pheromone to attract hornets for pollination [10].

V. velutina is the most notorious invasive species among hornets. Its invasion of Europe [11] severely affected European apiculture, leading to tens of millions of dollars in management costs [12]. While the Asian honeybee (Apis cerana) has evolved a special thermal defense against local hornets [13,14], the European honeybee (Apis mellifera) remains largely defenseless against this new predator [9,15,16], although a defense by asphyxiation has been observed in Cyprian honeybees, A. mellifera cypria [17]. Probably because of the weaker defense in A. mellifera than in A. cerana against hornet predators, the bee-hawking success rate is three times greater when V. velutina prey on A. mellifera than on A. cerana [18].

There are 22 morphologically distinguishable hornet species in the genus Vespa living in Southeastern Asia [19,20,21], but why does V. velutina appear more invasive than others? Of the 11 to 13 recognized subspecies [21,22,23], why was only the subspecies V. v. nigrithorax successful in Europe, South Korea and Japan [24,25,26,27]? Is it an accidental introduction that could have happened to any other hornet species or subspecies, or is it because V. v. nigrithorax possesses features that are more conducive to colonization? What features are required for a successful invasive species?

A successful invasive hornet would need to overcome at least seven obstacles to start a new population in a new habitat. First, the mated queen (known as foundress) needs to have a means of long-distance dispersal, either by self-propelled flight or by hitchhiking. Second, the queen needs to find a nest site and build its first nest. Third, she and her offspring must find suitable food. Fourth, The colony has to grow rapidly to avoid accidental loss of the colony. Fifth, because the foundresses are likely few or even just one, at least some of her fertilized eggs should be relatively free of deleterious recessive alleles so that some of her descendants will not suffer from inbreeding depression. Sixth, the colony needs to neutralize new hymenopteran pathogens [28] and escape the new predators such as the European honey buzzard [29]. Seventh, the newly produced queens need to find mates that are relatively free of deleterious recessive alleles to propagate the colonies into future generations.

Although hornets are not as strong fliers as locusts or katydids that can fly across the Pacific to colonize Hawaiian Islands [30], they are known to hitchhike over long distances, which must be how V. mandarinia foundresses reached North America from Asia [5,6,16,31]. Foundresses of Vespa velutina tend to choose small and dark cavities to hibernate [32]. They can therefore be transported long distances while hibernating. Hitchhiking to Europe must have occurred multiple times in V. velutina based on the observation recorded in the GBIF database [33]. V. velutina was first observed in France and Belgium in 1915 and 1923, respectively (Figure 1), but then must have disappeared. It was again observed in the Netherlands in 1980 (Figure 1) but again disappeared. It was also observed twice in 2003 and four times in 2004 according to GBIF records [33], but only the observation in 2005 was described in detail and published [34,35,36] with the conclusion that V. velutina was already well established in several counties in France. One should not confuse these locally collected specimens with those specimens collected elsewhere but stored in local museums. For example, the Natural History Museum stores a V. velutina specimen (NHMUK010636232) collected in Sumatra in 1914 which should not be misinterpreted as an observation of V. velutina in the UK in 1914.

While food is often a key environmental factor limiting the geographic distribution of many species, it is not in the case of hornets. V. velutina prey mostly on honeybees [24] which contribute two-thirds of the diet of V. velutina in an urban environment and about one-third in farmland and forests [24]. The global distribution of honeybees (Figure 2) shows the availability of food for hornets. The highest species diversity of honeybees is observed in Southeastern Asia (Figure 2), which corresponds to the subspecies richness in V. velutina (Figure 1).

After V. velutina was observed on Nov 1, 2005 [34,35], the invasive species quickly expanded its range in France and from France to other European countries. Where is the source population from which the first founding queen initiated the successful colonization? When did the foundresses of the last successful colonization arrive in France? Is it much earlier than 2005? How rapidly can the population size increase over the years? These can be addressed by estimating the initial invasion time and the instantaneous rate of growth.

The source population is often identified by using DNA barcoding data [37] when two requirements are met. The first is that the populations need to be genetically unique. If all individuals of V. velutina in Southeastern Asia are genetically identical, then the source population will encompass the entire Southeastern region, which is not informative. The second requirement is the characterization of the DNA signature of all genetically unique populations. If a population X of V. velutina is genetically different from all other populations and if the European V. velutina individuals share identical genetic signatures of this population X, then we may infer that the source population of the European V. velutina is population X. The method has been applied successfully to the identification of the source population for V. mandarinia individuals found in Washington State in the US and British Columbia in Canada, with the source population in South Korea and Japan, respectively [16,31].

2. Material and Methods

Invasion time (${\mathit{t}}_0$) and instantaneous rate of growth (r): The population size of one or more founding individuals will increase in size after the first colonization. If the population size is monitored over time long enough, then this change in population size over time can be modeled by a logistic equation, so that $t_0$ and r can be estimated:

$${\mathit{N}}_{\mathit{t}}={\displaystyle \frac{{\mathit{N}}_0\mathit{K}}{{\mathit{N}}_0+\left(\mathit{K}-{\mathit{N}}_0\right){\mathit{e}}^{-\mathit{r}\left(\mathit{t}-{\mathit{t}}_0\right)}}}$$

(1)

The French government has sponsored insect-monitoring programs to create and update a complete insect inventory ever since the 19^th century, especially during the Nepoleonic period when the British imposed naval blockade so that France was not able to import insect-derived products such as cochineal from the Americas [38]. The number of reported observations (n) of V. velutina over the years might be proportional to population size, i.e., $n_t=c{N}_t$, where c is assumed to be constant, and can be used to obtain a reasonable estimation of $t_0$ and r. The resulting K (carrying capacity) and N₀ (initial population size) are not meaningful because we do not know the constant c. The number of observations of V. velutina in France over years was downloaded from the Global Biodiversity Information Facility [33].

The assumption of a constant c above could be problematic. If a country initiates a large-scale effort to monitor insects in year T, then the number of observations of insects after T would be greater than that before T simply because of the change in observation effort, not necessarily because of the change in population size. There are also many other factors that could affect observation efforts, e.g., wars. In the French case, there are no known significant changes in observation effort so the assumption of a constant c might be reasonable. However, the spread of V. velutina in France and in other European countries, and the detrimental effect on apiary, has resulted in a large-scale coordinated effort to eliminate the invasive hornet in 2012. This is legitimized by the French government’s classification of V. velutina as an invasive alien species harmful to beekeeping, in response to the repeated requests from beekeepers' associations for several years. This classification allowed the development and implementation of mandatory control programs at the national and departmental levels. Although the elimination effort failed, the population size of V. velutina decreased dramatically, so I included only observation data up to 2011.

DNA sequence data and analysis: I downloaded the mitochondrial sequences and geographic coordinates from GenBank and the DNA barcoding Bold System [37]. One of the sequences annotated as from V. velutina (JQ780459) must have resulted from a species misidentification because the sequence is nearly identical to V. affinis but very different from all known V. velutina. I wish to highlight this partly because the sequence has been used in other studies [39,40] as a V. velutina sequence and partly because of the many misannotations in GenBank sequences [41]. For specimens with complete mitochondrial genomes, the COX1 sequences were extracted from the GenBank file with DAMBE [42]. COX1 sequences from V. simillima, V. bicolor and V. vivax were used as outgroups. Previous studies [43,44,45,46] have shown V. simillima, V. bicolor and V. vivax are closely related to V. velutina. V. simillima appears to occupy an ecological niche to V. velutina because the invasion of V. velutina into South Korea displaces the native V. simillima [27].

The Bold System does not include proper geographic coordinates for some specimens, but instead lists the center of the country as the geographic coordinates. These are obviously unsatisfactory and were replaced by approximate geographic coordinates of the sampling location recorded in the original publication. The sequence identification for downloaded sequences, the latitude and longitude, and GC% were included in Table S1 in the supplemental file.

The resulting COX1 sequences were aligned using MAFFT [47] with the most accurate LINSI option (‘–localpair’ and ‘–maxiterate = 1000’). For phylogenetic reconstruction with PhyML [48], the GTR + Γ model was used with four discrete rate categories for approximating a continuous gamma distribution [49]. This model was chosen based on the information-theoretic index AIC and the likelihood ratio tests [50,51] among the nested HKY [52], TN93 [53], and GTR [54,55] models with or without the discrete gamma distribution to accommodate rate heterogeneity in substitution rate among sites. However, the TN93 + Γ model generated the same tree with negligible difference in branch lengths. The tree improvement option (‘-s’) was set to ‘BEST’ (best of NNI and SPR search). The ‘-o’ option was set to ‘tlr’, which optimizes the topology, the branch lengths, and rate parameters. MAFFT and PhyML are included in DAMBE and called to analyze sequences with a consistent user interface.

For sister lineages that have evolved over the same amount of time, a sister lineage with a wider geographic distribution than the other indicates that the former is likely more invasive than the latter. To measure the area of geographic distribution of a taxon, I used the convex hull algorithm [56] implemented in the PGT software version 1.0.0 [57]. Approaches to alleviate the confounding factors were taken and discussed in the result interpretation.

Geophylogenies were produced with PGT [57] which takes two types of information: a phylogeny and a list of specimens with associated latitude and longitude. When specimens from two subspecies share exactly the same latitude and longitude, the map marker of one specimen would completely cover the other. For this reason, their geographic coordinates are slightly shifted from each other to avoid complete overlap. The objective of a geophylogeny is to visualize phylogenetic relationships and biogeographic distribution of evolutionary lineages [30,31]. We used PGT software version 1.0.0 [57] to generate geophylogenies for visualization. PGT makes use of both Google Maps and Microsoft Bing Maps with regular map view and satellite terrain view.

3. Results

3.1. Estimating Invasion Time ($t_0$) and Instantaneous Rate of Growth (r)

A logistic equation was fitted to the number of observations of V. velutina in France over years downloaded from GBIF [33]. Only data from 2003 to 2011 were used because The French government officially classified V. velutina as an invasive pet to facilitate the implementation of mandatory control programs at the national and departmental levels [58]. An independent invasion event has happened in South Korea by the same subspecies, V. v. nigrithorax [27], with the first observation recorded in 2005 [33]. It has since increased in number and displaced native hornet species with similar ecological requirement and life history straits, such as V. simillima [27]. I estimated $t_0$ and r for both.

As I have explained in the Methods and Materials, if we take the risk of assuming that the number of observations (N_obs) of V. velutina recorded in GBIF [33] is proportional to population size (N_t), then we can estimate the four parameters in Eq. (1), i.e., $N_0,K,t_0$ and r, where $t_0$ is the time of the last successful invasion, and r is the instantaneous rate of population growth. The four parameters were estimated by the least-squares method (Figure 3). $t_0$ was in early 2002.40, and $r=1.3667$, in France (Figure 3A). The corresponding values for South Korea are 2002.77 and 0.2812, respectively (Figure 3B). Because we assume that $N_{obs}=c{N}_t$ but we do not know the constant c, the values of $N_0$ and $K$ in Figure 3 are not meaningful because the true $N_0$ and $K$ should be multiplied by the constant c that we do not know.

Two points are worth making. First, N_obs is much greater for France than for South Korea, so the observed data points are much closer to the fitted curve in the population in France than that in South Korea. Second, r is much smaller for the population in South Korea than in France, which is consistent with the following interpretation. There are at least seven Vespa species and three subspecies in South Korea as reported in 1994 [26], so the ecological niche for Vespa species is perhaps already crowded. This is consistent with the observation that the expansion of V. velutina in South Korea is associated with the displacement of the native V. simillima [27]. This interspecific competition would slow down the growth in population size for V. velutina as a new immigrant. In contrast, there is only one native hornet species, V. crabro, in France, and V. velutina and V. crabro do not affect each other negatively [60], suggesting niche divergence between the two species.

I should emphasize the potential problems with the estimated parameters in (Figure 3). First, the assumption of $N_{obs}=cN$, where N is the population size, may not be true as one could list multiple factors that could lead to violation of the assumption. Second, even if the assumption is roughly true, the estimation in Figure 3B could be unreliable because of the small N_obs values as well as 0 observations in seven consecutive years (Figure 3B).

3.2. Phylogeographic Analysis

A phylogeographic analysis is fundamental for understanding invasive species. When two sister species (X and Y) diverged from a common ancestor, they would spread to their suitable habitats. After the same amount of divergence time, if species X gains a much wider geographic distribution than species Y, then we may infer that species X has a greater invasiveness than species Y.

V. v. nigrithorax is the subspecies that invaded Europe [34,35] and South Korea [59], and it features the widest distribution among all V. velutina subspecies [22] even without counting the area it recently colonized in Europe. This is also clear from the distribution of the four most sampled subspecies (Figure 4). However, it would be wrong to therefore conclude that V. v. nigrithorax is the most invasive. There are two confounding factors. The first is the timing of the subspecies recognition which could be controversial [20]. If a subspecies were only recognized yesterday, then its samples would be few, and its distribution represented by its small samples would necessarily be much more limited than a subspecies that have been recognized and well-sampled for a hundred years. This problem can be overcome by using sampling points after the subspecies of interest have all been recognized. The other confounding factor is the divergence time. If subspecies X originated one million years ago whereas subspecies Y originated only 1000 years ago, then subspecies X is expected to have a wider distribution than species Y even if they are equally invasive. Thus, valid comparisons in their invasive potential can only be made between sister lineages that have diverged for the same amount of time. For these reasons, phylogenetic analyses delineating sister lineages are essential for understanding invasive species.

The sister subspecies of V. v. nigrithorax is V. v. variana, both being within the blue-shaded clade (Figure 5). Visual inspection of Figure 4 suggests that V. v. nigrithorax has a wider distribution than V. v. variana. To obtain a quantitative measure of geographic distribution for comparison, I have computed the area of geographic distribution by the convex hull algorithm [56] implemented in PGT [57]. In calculating the area for V. v. nigrithorax, I excluded samples from both Europe as well as those from Japan and South Korea because these also resulted from recent invasion events in South Korea [26,27] and Japan [61] [25,40,62]. The resulting area is 360.5129 square degrees (because the input is latitude and longitude in degrees) for V. v. nigrithorax and 47.1508 for V. v. variana, consistent with the interpretation that V. v. nigrithorax tends to be more invasive than V. v. variana.

One problem with the comparison above is the difference in sampling points between the two subspecies. V. v. nigrithorax has 170 sampling points whereas V. v. variana has only 20. If one sampled only 20 houses in Ottawa, the area enclosed by the convex hull polygon would underestimate the geographic distribution of Ottawa residents far more than an equivalent polygon from 2000 houses. This effect can be visualized if we subsample the 170 sample points for V. v. nigrithorax (Figure 6). When the number of sample points increases, the area enclosed by the resulting convex hull polygon also increases. For this reason, it is not fair to compare the area enclosed by 170 sampling points for V. v. nigrithorax with an area enclosed by only 20 points V. v. variana.

The phylogenetic tree (Figure 5) is consistent with previous studies [39,40]. All specimens from South Korea and Japan share the same haplotype with two specimens from Zhejiang, China (JQ780452 and JQ780449), all colored red in Figure 5. All specimens from Europe share the same haplotype as two specimens from Jiangsu and Zhejiang in China (JQ780451 and JQ780453), all colored in blue in Figure 5. The colonization of Europe and South Korea/Japan is by different lineages (Figure 5), i.e., by independent colonization events. The V. velutina found in South Korea and Japan are geographically close and located at two locations connected by high volumes of transportation, so the V. velutina in South Korea may be the secondary source population for those established in Japan [40].

Figure 6 points to an approach for a fair comparison. When taking a subsample of 20 points from the 170 points for V. v. nigrithorax, the area calculated from this subsample is comparable to that from the 20 points for V. v. variana. The average area of repeated subsamples is 63.0385, with the 95% confidence interval of (57.3055, 304.4059). Because the calculated area for V. v. variana is only 47.1508, lower than the 95% lower limit, we may conclude that the distribution area of V. v. nigrithorax is indeed wider than that of V. v. variana. I believe that this conceptual framework can be used to compare invasiveness between any pair of sister lineages.

The phylogenetic tree of V. velutina is generally consistent with other studies. It was previously reported that V. velutina first diverged into two clades: an Indonesian-Malaysian clade and an Asian Continent Clade [40,63]. The latter lineage and subsequent derivatives colonizing Pakistan and India are generally sedentary, but a new invasive lineage (shaded in blue, Figure 5) arose, moved east and colonized most of Southeast Asia. This is better visualized with a geophylogeny (Figure 7). The first split is between the Indonesia/Malaysian clade and the continental clade. The second split is between the sedentary India/Pakistan clade and an invasive clade that expanded not only to Southeast Asia including Japan and South Korea but also to Europe (Figure 7).

Contrasting the geographic points in Figure 4 and Figure 7, one would immediately notice the scarcity of points in Figure 7, i.e., many V. velutina populations are still not represented by molecular sequences. The scarcity of sequences would cause difficulties to trace the source population. Two conditions are required to trace an invasive individual to its source population, i.e., populations should be genetically distinct and should be genetically characterized. If the two sequences from Jiangsu and Zhejiang in China (JQ780451 and JQ780453, colored blue in Figure 5) were not sampled, we would have to infer that the source population is represented by the three sequences colored in green (Figure 5). These three sequences were all sampled from Yunnan in Southwest of China, thousands of kilometers from Jiangsu and Zhejiang.

When invasive individuals of V. mandarinia were caught in the Washington State and British Columbia, they were found to be genetically nearly identical to specimens sampled in South Korea and Japan, respectively, and the source populations were therefore inferred to be somewhere in South Korea and Japan [16,31]. However, the populations along the eastern coast of China were poorly sampled. It is possible that the invasive individuals are genetically more similar to those V. mandarinia populations in the eastern coast of China. I highlight this in the hope that more targeted samples will be sequenced in the future.

4. Discussion

Invasive species are interesting mainly because of three consequences [64]. The first is the ecological and economic consequence mediated by the damage caused by the invading species on the fauna and flora, as well as threats to human health. V. velutina is a mass murderer of honeybees causing significant damage to the honey and bee industry not only in Asia but more so in Europe [11,12,44,45,65,66]. A. mellifera has come in contact with hornets in two ways. The first is through the introduction A. mellifera into Japan and South Korea where hornets exist either natively or through accidental introduction [40]. The second is through the accidental introduction of Asian hornets into Europe. If different honeybee colonies in Europe exhibit anti-hornet behaviors with different efficiency, then one should be able to observe the effect of selection mediated by the hornet predation on honeybee evolution, i.e., those with more efficient anti-hornet behaviors would be favored by the selection.

The invasion of V. velutina may also exert selection pressure on local competitors, e.g., V. crabro in Europe or the multiple hornet species in South Korea. This may lead to three possible outcomes. First, when the population size is small for both species, then they may both increase in population size with increasing availability of prey species and other food items. When the population size of V. velutina (N_Vv) and that of V. crabro (N_Vc) are both high, and if they share the same ecological niche and the same carrying capacity of (N_Vv + N_Vc), then one would expect N_Vv to be negatively associated with N_Vc. Specifically, one would expect a successful invasive species to displace the native with a nearly identical ecological niche. This has been observed in South Korea in which the invasive V. velutina displaces the native V. simillima [27]. However, there is no clear indication that the invasive V. velutina in Europe displaces native European hornet, V. crabro [60], suggesting that V. velutina in Europe does not face a local competitor sharing a similar ecological niche. This interpretation is consistent with the instantaneous rate of growth (r, Figure 3). The invasive individuals in Europe and South Korea should have similar biological features. However, r is much greater for V. velutina in Europe (with little competition) than in South Korea (with strong competition) (Figure 3). The second possible outcome of competition between an invasive species and its native competitor is that the two may be expected to differ in seasonal activities or daily rhythms, and they may diverge in nest sites and diet. The third possible outcome is that they may exhibit differential habitat preferences. In addition, the invasive V. velutina also provides new food for local predators such as the European honey buzzard, Pernis apivorus [29]. All these provide research opportunities for evolutionary ecologists.

The direct effect of V. velutina on human health appears minor. While the related Asian giant hornet, V. mandarinia, was known to attack humans and can kill the victim with an average of only 59 stings [67], V. velutina has not been a threat to human health. In France, only one severe case of envenomation was attributed to V. velutina, when the victim was stung 12 times on the head [68].

The second consequence of an invasive species is [69,70,71]. As Wright (1942) envisioned, a small population made of one or a few foundresses would facilitate the fixation of non-adaptive alleles or allele combinations. In most cases, the few foundresses and their descendants may just go extinct, leading to failed invasions. As seen in Figure 1, the observations of V. velutina in 1915, 1923 and 1980 likely represent failed invasions. However, if the few foundresses managed to generate a viable population, then any existing genetic variation or new mutations would evolve in a unique genetic background independent of the source population in a new habitat potentially with a selection regime different from that of the source population. This may lead to a new evolutionary trajectory and facilitate speciation.

The source population may generate many invasive foundresses to explore new habitats. Some of these foundresses may carry many deleterious recessive alleles and their offspring would suffer from inbreeding depression and go extinct. Some other foundresses may harbor few deleterious recessives and succeed in colonizing new habitats. This inter-colony selection is depicted in Figure 8. A colony that has purged off most of the deleterious recessive mutations can spread and replace colonies with more deleterious recessive mutations. Such a successful colony would generate successful invasive foundresses because their descendants will have a better chance of escaping inbreeding depression when a few foundresses or just a single foundress moves into a new habitat with all genetic variation limited within her fertilized eggs. In this context, one may expect a successful invasive species should have low instead of high genetic variation, i.e., the descendants should distribute widely but with little genetic variation. V. v. nigrithorax fits these two criteria as one can see in Figure 4 and Figure 7.

The third consequence of invasive species is the alteration of established biogeographic patterns mediated by human-facilitated dispersal. For example, the native distribution of V. mandarinia is in Asian countries and the Russian Far East [20,72,73] and is unlikely to disperse from Asia to North America by natural means. However, globalization and modern cargo transportation would make this possible, and multiple incidents of V. mandarinia queens found in British Columbia in Canada and Washington State in the USA have been documented and traced to source populations in Asia [4,16,31]. Vespine hornets appear to be highly capable of survival and reproduction in new environments and include multiple invasive species [2,3].

There has been no explicit index of invasiveness, but my results suggest the conceptual framework in which such an index can be formulated. Designating the distribution area of a taxon as A, and the divergence time of the taxon as T, the index of invasiveness (I_inv) can be expressed as

$${\mathit{I}}_{\mathit{i}\mathit{n}\mathit{v}}={\displaystyle \frac{\mathit{A}}{\mathit{T}}}$$

(2)

The area A could be measured by the convex hull method that I have used in the paper and T can be from phylogenetic dating. For taxa that evolve with a roughly constant rate, then T can be operationally replaced by branch lengths for comparison of I_inv between sister lineages.

5. Conclusions

The time of the last successful invasion of Vespa velutina into France and South Korea is 2002.40 and 2002.77, respectively, estimated by fitting a logistic equation to the number of observations over time. The instantaneous rate of increase is 1.3667 for V. velutina in France and 0.2812 in South Korea, consistent with the interpretation of little competition in France and strong competition from local hornet species in South Korea. The invasive potential of two sister lineages can be compared by their distribution area when proper statistical adjustments are made to account for differences in sample size. The ancestor of V. velutina split into two lineages, one found in Indonesia/Malaysia and the other colonizing the Asian continent. The second lineage split into a sedentary clade inhabiting India and Pakistan and the other invasive one colonizing much of the Southeast Asia. This latter lineage gave rise to the subspecies V. v. nigrithorax that invaded France, South Korea and Japan. My software PGT, which generates geophylogenies and computes geographic areas for individual taxa, is useful for understanding biogeography in general and invasive species in particular.