Mitogenomics and the Global Dispersion of Vespula germanica: A Case Study from South Africa Shows Evidence for Two Separate Invasion Events

Abstract

Vespula germanica is currently present in all major world regions outside its native Northern Hemisphere range and poses a biological threat to the invaded ecosystems. The genetic diversity of the species is poorly described in both the native and invaded ranges, thus hampering insights into possible mechanisms of invasion. In South Africa, V. germanica was first detected in 1972, and a recent study concluded that one large or several independent invasion events had occurred. However, the high number of low-frequency haplotypes reported therein raised doubts about the quality of the data. In this study, we reassessed the haplotype diversity of V. germanica in South Africa under improved methodological conditions. New mitochondrial markers were developed using complete mitochondrial genomes of V. germanica that allowed the identification of polymorphic regions and the design of robust species-specific primers. Contrary to two previous studies, only two mitochondrial haplotypes were found in South Africa despite almost doubling the number of sampled nests. It is likely that that the number of haplotypes previously reported was overestimated due to the miscalling of nucleotide positions in the electropherograms. Furthermore, the two haplotypes found have contrasting geographic distributions, which supports the known invasion history for this species. Availability of complete mitochondrial genomes for selection of polymorphic regions and design of robust species-specific primers improved the accuracy of the assessment of V. germanica diversity in South Africa. This approach will also be valuable for studying invasive wasp populations of this and other species globally.

1. Introduction

In a globally connected world with an ever-increasing risk of biological invasions [1], understanding mechanisms underlying the spread of invasive species is of utmost importance [2,3,4]. For example, knowledge on dispersal pathways, propagule pressure, and establishment success can be used to inform biosecurity protocols and to manage existing invasions more effectively [5,6]. Furthermore, investigating the origin of invasive species in non-native regions helps to determine the most likely routes of invasion, which potential source populations are best suited to receiving countries, and can clarify the number of invasion events that took place in a given region [7,8,9]. Notably, regions previously unsuitable may become habitable for invasive species such as social wasps due to global climate change [5,10], thus increasing the number of potential source populations from the native range or other invaded regions, i.e., bridge-head effects.

Vespula germanica L. (Hymenoptera: Vespidae) is native to the Palearctic region and currently represents the most widely distributed of all invasive social wasps [11]. Its current invaded geographic range is remarkably wide, spanning from Finland [12] and the Faroe Islands [13] in the Northern Hemisphere to South America as far south as Patagonia [14], and Oceania (New Zealand) [15]. Vespula germanica was first reported in South Africa approximately 50 years ago, in the Western Cape Province [16].

The five-decade-old presence of V. germanica in South Africa is a particularly interesting case study, as the expansion has been slow and geographically constrained. The wasp was confined to Table Mountain and Peninsula from 1972 until 2001, when it was suddenly found at several new localities 50 km or more from the original point of detection [17,18,19]. Subsequent to this jump dispersal event, there was evidence of further jump dispersal to outlying towns averaging 100 km away from the original invasive population. There is also evidence of ongoing invasion (2018–2020), as reproductive nests were found for the first time at new outlying localities [19]. Most nest records show an association with river courses and may be due to individual or combined effects of cooler temperatures more favorable for dispersal of queens and/or foraging wasps, and greater prey availability [19].

Recently, Rankin [20] highlighted the utility of mitochondrial, genomic, and methylation data for gaining insights into recent and established eusocial wasp invasions. Patterns of mitochondrial haplotypes are particularly useful for clarifying invasion history and estimating the minimum number of founder queens. Thus far, mitochondrial haplotype data that could inform patterns of worldwide dispersal of V. germanica have been restricted to a small portion of the mitochondrial genome (~1550 bp of cytochrome oxidase I and II, and cytochrome b), and the haplotype structure across the native and invaded geographic range has not been extensively described [9,21].

Two previous studies based on the number and diversity of mitochondrial haplotypes of South African specimens of V. germanica concluded that the species had established in the country as a result of a single large invasion event or multiple invasion events that brought single haplotypes [9,21]. Brenton-Rule et al. [21] and Eloff et al. [9] identified four and six haplotypes from nine and forty-two nests, respectively, and the results pointed to an origin of the invasion in central Europe. The haplotype accumulation curves presented by Brenton-Rule et al. [21] indicated that, compared to other invaded countries, South Africa had the highest number of unique haplotypes relative to the number of specimens analyzed, although Eloff et al. [9] showed that maximum diversity had been sampled in their study. However, the presence of multiple, closely related haplotypes reported in these two studies was unexpected given the recent presence of V. germanica in the country, thus justifying a follow-up study.

We report a reassessment of mitochondrial haplotypes of V. germanica across its present geographic distribution in South Africa using informative markers selected by comparison of complete mitogenomes.

2. Materials and Methods

2.1. Specimen Collection and DNA Extraction

Adult specimens of V. germanica (Figure 1) were collected at several locations in the Western Cape province of South Africa, and by collaborators in Argentina and New Zealand (Table S1). Specimens were stored individually in 100% ethanol at −20 °C before DNA extraction. DNA was extracted from a section of the thorax of each individual specimen using a standard phenol-chloroform method [22] and stored at −20 °C until downstream analyses.

2.2. Next-Generation Sequencing, Assembly and Annotation of Mitogenomes

Preliminary haplotyping for the selection of specimens from South Africa for next-generation sequencing (NGS) was based on DNA barcoding (data not shown) of 10 wasps from different nests that resulted in the identification of two haplotypes consistent with the two most frequent haplotypes found by Eloff et al. [9]. The mitogenomes of four wasps (one from each haplotype found in South Africa, one from Argentina, and one from New Zealand) were recovered for identification of polymorphic regions and conserved regions for design of species–specific primers. NGS was performed using the Ion Torrent ProtonTM sequencing platform (ThermoFisher Scientific, Waltham, MA, USA) available at the Central Analytical Facilities of Stellenbosch University, South Africa. In brief, sequence libraries were prepared using the NEXTflex™ DNA Sequencing Kit for Ion Platforms (PerkinElmer, Waltham, MA, USA) according to the BI00 Scientific v15.12 protocol. Libraries were pooled and sequenced using the Ion PI HiQ™ Sequencing Solutions Kit (Life Technologies, CA, USA). The reads for each wasp were mapped to the mitogenome of V. germanica (GenBank KR703583) using Geneious Prime v2021.1 (https://www.geneious.com, accessed on 25 January 2022) with the mapping parameters set for the medium sensitivity/fast option. The consensus sequences were annotated with the MITOS Web Server (http://mitos.bioinf.uni-leipzig.de/index.py, accessed on 25 January 2022) [23] and the ARWEN software (http://130.235.244.92/ARWEN/, accessed on 25 January 2022) [24], using the invertebrate mitochondrial genetic code. The positions of the transfer RNA genes (tRNAs) predicted by MITOS and ARWEN were compared, and the most probable predictions were manually selected. The positions of the ribosomal RNA genes (12s and 16s rRNAs) were determined by comparison with other Vespidae mitogenomes. The large noncoding region downstream of 12s rRNA was annotated as the AT-rich region. The new four mitogenome sequences were deposited on GenBank under the accession numbers OM735804 to OM735807.

2.3. Phylogenetic Tree

The phylogenetic relationships among the four new mitogenomes and one V. germanica sequence from Asia (GenBank KR703583) were reconstructed based on the 13 protein-coding genes using the maximum likelihood method. Sequences were aligned using the MAFFT algorithm [25] on Geneious Prime. The tree was run on the PHYML online server (http://www.atgc-montpellier.fr/phyml/, accessed on 30 January 2022) [26], using the Smart Model Selection and the fast likelihood-based test (aLERT) [27]. The final tree was drawn using FigTree v1.4.4 (http://tree.bio.ed.ac.uk/, accessed on 30 January 2022).

2.4. Identification of Mitochondrial Polymorphic Regions

Intergenic spaces, tRNA and rRNA genes, and the AT-rich region were manually excised from the alignment to exclude regions of lower sequence coverage that could include artefactual polymorphisms. Polymorphic regions were identified on the alignment of protein-coding genes (PCGs) using the DNA Polymorphism function of DNAsP6 with a 100-bp sliding window and a 25-bp step size [28]. As our study focused on the diversity of V. germanica in South Africa, regions that were only polymorphic in the wasps from Argentina and New Zealand were not included in the final panel of markers.

2.5. PCR Amplification and Sanger Sequencing of Polymorphic Regions

Primers for PCR amplification and Sanger sequencing were manually designed to have a perfect match with conserved regions upstream and downstream of the polymorphic regions across all mitogenomes, and to generate amplicons smaller than 700 bp for unidirectional sequencing or smaller than 1300 bp for bidirectional sequencing (

Table 1. Primers for PCR and Sanger sequencing of five polymorphic mitochondrial regions in Vespula germanica specimens collected in South Africa. * Sequencing primers.

Region	Location	Primer	Sequence (5′-3′)	Size (bp)
Ves1	COX1	Ves1-F	GATTTTGATTATTACCTCCATCC	686
		Ves1-R *	ACCGTAAATTGTTGCTAATCATC
Ves2	COX3	Ves2-F	ACCAATGATGACGAGACGTAG	550
		Ves2-R *	GCAGCTTCGTATCCAATATGG
Ves3	ND5	Ves3-F *	CCTCATATTCTAGAAAGATATAAC	1145
		Ves3-R *	TGTCCAAATATTATGGGTTTATTG
Ves5	CytB	Ves5-F *	TCTGCTATTCCTTATATTGGCC	608
		Ves5-R *	TTTGAGTTAATTCAATGAAAGGTG
Ves6	COX1-COX2	Ves6-F *	CTCGTGCATATTTTACATCTGC	1189
		Ves6-R *	ATCTGTTGATGTTGTTAGAATTC

). PCR amplifications were performed in a total volume of 5 µL containing 1x of QIAGEN Multiplex PCR Kit (QIAGEN), 0.2 μM of each primer, 0.5 μL of MilliQ H₂O, and 1.0 μL of template DNA (~100 ng/µL). Thermal cycling was performed as follows: 15 min at 95 °C; 35 cycles of 30 s at 94 °C, 90 s at 54 °C, 90 s at 72 °C; and 10 min at 72 °C.

Sanger sequencing was performed unidirectionally or bidirectionally (Table 1), following the manufacturer’s protocol for the BigDye v3.1 Cycle Sequencing Kit (Applied Biosystems). Sequencing reaction products were run on an ABI3170xl Sequencer (Applied Biosystems) at the Central Analytical Facilities of Stellenbosch University, South Africa. Sequences were manually inspected, curated for ambiguities, and translated into amino acid sequences for assessment of the presence of premature stop codons and/or frameshift mutations suggestive of spurious amplification of noncoding regions. The sequences of the polymorphic regions for each individual wasp were concatenated using Geneious Prime and treated as a single marker in downstream analyses. Individual sequences for all specimens were deposited on GenBank under the accession numbers OM719750 to OM719878.

2.6. Haplotype Analyses

The dataset generated by Eloff et al. [9] was downloaded from GenBank, and the three sequences of each specimen from South Africa (n = 58) were concatenated in a single haplotype (1597 bp) for re-analyses. The new dataset generated in our study included a total of 78 individual wasps (South Africa, n = 76; Argentina, n = 1; New Zealand, n = 1) sequenced for a combination of mitochondrial regions (Figure 2). First, the complete mitogenomes of two wasps from South Africa, one from Argentina, and one from New Zealand were sequenced for identification of polymorphic regions and primer design, as described above. Following this step, 16 wasps were sequenced for the polymorphic regions Ves1, Ves2, Ves3, and Ves5 in addition to the monomorphic Ves6 (total of 3798 bp) covered in Eloff et al. [9] to ensure maximum overlap between their sequences and the new dataset. Then, 11 wasps from the same nests reported in Eloff et al. [9] were sequenced for Ves1 and Ves6 (total of 1086 bp) for confirmation of the suspected artefactual haplotypes, which were not found in the new specimens. In light of these results, an additional 47 wasps were only sequenced for Ves1 (618 bp) because this short polymorphic region, which contains four SNPs, is sufficient for discriminating between the two haplotypes found in South Africa. The geographic distribution of haplotypes was represented in a map generated in ArcGIS.

3. Results

3.1. Mitochondrial Genomes, Phylogeny of Vespula germanica, and Polymorphic Regions

The NGS runs generated an average of 14 million reads with an average size of 162 bp for each specimen, and an average of 1.8 million reads were mapped to the reference sequence. The average sequence coverage was high (115x), but four tRNA genes located downstream of the A+T-rich region were not recovered: tRNA^Ile (I), tRNA^Gln (Q), tRNA^Met (M), and tRNA^Tyr (Y). These genes were also not identified in the mitogenome of V. germanica from Asia available on GenBank (KR703583) or in Vespula vulgaris [29], but were annotated in Vespula flaviceps [30]. The new mitogenomes of V. germanica had the same gene rearrangement as V. vulgaris and V. flaviceps, with two differences compared with the hypothetical Arthropoda ancestor: a translocation or remote inversion of tRNA^Tyr from the W-C-Y cluster to the region of the I-Q-M cluster and the translocation of tRNA^L1 from downstream to upstream of ND1 (Figure 3).

The phylogenetic reconstruction of the relationships among V. germanica using the currently available mitogenomes evidenced the divergence between the Asian sequence and the cluster formed by South Africa, New Zealand, and Argentina, with maximum nodal support (Figure 4). Although the data available for phylogenetic analyses were very limited, the tree indicated that V. germanica is phylogeographically structured between East (Asia) and West (Europe) and supported a European origin for the invasions in South Africa, New Zealand, and Argentina, as opposed to an Asian origin.

The four mitogenomes were highly similar with a total of 38 SNPs in the 13 PCGs. The number of SNPs between all pairs of sequences ranged between 17 (Argentina/South Africa VG6) and 20 (Argentina/New Zealand). The regions with the highest number of polymorphisms (Pi ≥ 0.01) were located in the COX1 (Ves1; 686 bp), COX3 (Ves2; 550 bp), ND5 (Ves3; 1145 bp), ND4 (Ves4; 438 bp), and CYTB (Ves5; 608 bp) genes (Figure 5). The region spanning COX1-tRNA^Leu-COX (Ves6; 1189 bp) was sequenced to ensure full coverage of the region surveyed by Eloff et al. [9], although it contained no SNPs in the the complete South African mitogenomes sequenced in our study.

3.2. Haplotypes of Vespula germanica in South Africa

A total of 76 wasps from 72 nests were haplotyped, including 88% (36/41) of the nests sampled in Eloff et al. [9], and represent a much larger dataset and a higher coverage of the geographic distribution of V. germanica in South Africa compared to the previous study. Overall, the sequencing of a combination of selected mitochondrial regions showed two haplotypes, of which Hap1 was present at a higher frequency (0.76) than Hap2 (0.24) (Figure 2a;

Table 2. Haplotyping results of Vespula germanica in South Africa where 88% (36/41) of the nests surveyed in the previous study [9] were also analyzed in the present work.

Eloff et al. [9]	n	Freq
Hap1	43	0.74
Hap2	8	0.14
Hap3	2	0.03
Hap4	2	0.03
Hap5	2	0.03
Hap6	1	0.02
Total number of individuals	58	1.00
Total number of nests	41	-
This study	n	Freq
Hap1	58	0.76
Hap2	18	0.24
Total number of individuals	76	1.00
Total number of nests	72	-

). This result disagrees with Eloff et al. [9], who also found Hap1 and Hap2 but reported four additional haplotypes at low frequency. The mitochondrial regions sequenced in Eloff et al. [9] were covered in our study, and the polymorphisms in COX2 and CytB that defined the low-frequency haplotypes in a total of seven wasps reported in the previous study (SA22, SA23, SA271, SA482, SA461, SA481, and SA491) were not found (Figure 2a), although wasps from the same nests were analyzed. Furthermore, most of the wasps with low-frequency haplotypes reported in Eloff et al. [9] were found in the same nest along with wasps with other haplotypes; SA271 (Hap3) and SA272 (Hap1) were from the same nest, as well as SA481 (Hap5) and SA482 (Hap3), and SA491 (Hap6) and SA492 (Hap1). This unexpected result was not mentioned by the authors. The SNPs that define the low-frequency haplotypes reported in Eloff et al. [9] are located in regions of the sequences that are prone to low quality nucleotide calling; some polymorphisms are adjacent to a homopolymeric T-stretch in COX2, and others at the end of the sequence in CytB (Figure 2a). The low-frequency haplotypes reported in Eloff et al. [9] are closely related to Hap1, and Hap1 and Hap2 appear separated by two mutational steps. The network of the haplotypes found in our study, including Argentina and New Zealand for context, showed a closer relationship between Hap2 and Argentina than between Hap2 and Hap1 (Figure 6).

3.3. Geographic Distribution of Haplotypes in the Western Cape

The geographic distribution of the two haplotypes found in the nests sampled across the distribution of V. germanica in the Western Cape was highly structured (Figure 7). Hap1 was the most frequent haplotype overall, and it was found across the whole range, but at a very low frequency in the Table Mountain complex where Hap2 dominated. Hap2 had a lower frequency outside the Table Mountain complex and was only found in four nests in the greater Stellenbosch area. This phylogeographic structure results in an East–West split with Hap1 dominating in the East and Hap2 in the West. It was also evident that nests outside these two main distribution clusters were Hap1, except for the two nests sampled in the northern Cape Peninsula, which were Hap2.

4. Discussion

Due to their generally uniparental mode of inheritance without recombination, animal mitochondrial genomes are effectively haploid and can offer insights into the genetic diversity and phylogeographic structure of populations [31,32]. This application is particularly useful in the reconstruction of biological invasions, as these typically start with a small number of individuals that represent a subset of the diversity of the species in the original geographic range. Therefore, tracing back the source of biological invasions relies on robust genetic markers, well-characterized populations in the native range, and comprehensive sampling across the invaded range. Although V. germanica has invaded all major world regions and is a well-documented invasive species, these conditions have not been met and the origins of populations in invaded regions worldwide are poorly understood.

Prior to this work, only two studies have attempted to gain insights into the phylogeographic structure of V. germanica. In 2018, Brenton-Rule et al. [21] presented a spatially structured mitochondrial haplotype network of wasps collected in the native range (mainland Europe and the United Kingdom, but omitting Asia), and four regions in the invaded range (Argentina, Australia, New Zealand and South Africa). As expected, the native range had higher haplotype diversity than the invaded regions, but Europe and the UK did not share haplotypes. Europe shared one haplotype with Argentina and another with Australia, and the UK shared one haplotype with Australia and another with New Zealand. South Africa had four exclusive haplotypes: one dominant haplotype from which two singletons were derived by one mutational step, and a second haplotype which diverged from the rest. More recently, the genetic diversity of V. germanica in South Africa was assessed using mitochondrial sequences and microsatellite markers [9]. While microsatellite-based analyses were largely inconclusive, mitochondrial analyses indicated that V. germanica in South Africa had originated from two to six queens. The authors found six haplotypes among 42 nests: one haplotype at high frequency (0.74) from which four haplotypes represented by one or two specimens were derived, and one diverged haplotype at a lower frequency (0.14). Given the recent presence of V. germanica in South Africa, the four haplotypes separated from the most frequent haplotype by one- or two-step mutations could represent sequencing artifacts. As Eloff et al. [9] and Brenton-Rule et al. [21] followed similar methodological approaches in terms of PCR and sequencing, it is likely that the dataset generated by the latter also included artefactual haplotypes. Unfortunately, it was not possible to re-analyze the dataset generated by Breton-Rule et al. [21] because the sequences are not publicly available, which is contrary to good practices in reporting of sequence data; therefore, our comparisons can only include the dataset presented by Eloff et al. [9].

Overall, our re-analyses indicated that the haplotype diversity reported by Eloff et al. [9] was overestimated due to the presence of sequencing artefacts interpreted as legitimate polymorphisms. These sequencing artefacts likely resulted from the use of nonspecific PCR primers developed in previous studies for amplification of mitochondrial genes in a wide range of species [33] and primers specific for V. vulgaris [29] at low annealing temperatures (45 °C and 47 °C) [6]. Nonspecific PCR primers are employed when the sequence of the species of interest is not known, and this strategy often produces spurious amplification of nontarget regions and low-quality sequences with ambiguous nucleotide positions [34]. The SNPs that Eloff et al. [9] reported approximately in the middle of COX2 are located close to a stretch of eight T nucleotides, which also creates opportunity for Taq polymerase stuttering and generation of PCR artefacts. In our analyses, the marker that covered the COX2 region analyzed by Eloff et al. [9] (Ves6) was monomorphic in all wasps sequenced for this region (n = 29), including specimens from the nests where the low-frequency haplotypes were reported in the previous study. Eloff et al. [9] also found SNPs in CytB that we did not find in the wasps sequenced for this region (Ves5; n = 18), which included specimens from the nests where the low-frequency haplotypes were reported in the previous study. In this case, the artefactual SNPs likely resulted from deficient curation of low-quality nucleotide calling at the 5′ end of the sequences. Eloff et al. [9] also did not mention that they found three cases of different haplotypes from wasps collected from the same nest, which suggests that this peculiar result evaded their attention. All three cases involved one wasp with Hap1 and another with a secondary/artefactual haplotype derived from Hap1 by one or two mutational steps. This again suggests that experimental errors were responsible for overestimating the number of haplotypes found in South Africa.

In summary, Eloff et al. [9] reported four low-frequency haplotypes (0.03–0.02) in seven specimens (12.1% of the total number of specimens), which we believe were defined by artefactual nucleotide variants. We did not find these haplotypes in our study, although we covered 72 nests, including 88% of the same nests sequenced in Eloff et al. [9]. The high number of low-frequency haplotypes can be explained by the reliance of the previous study on nonspecific PCR primers used at low annealing temperatures. This strategy may have resulted in low-quality sequences that were not correctly curated for ambiguous positions. While manual inspection of Sanger sequences is time-consuming and requires some level of expertise in the interpretation of chromatograms, sequence curation is vital for generating high-quality data. Therefore, nucleotides that differ from the most frequently found variant must be carefully checked for ambiguities and may have to be confirmed by sequencing in the opposite direction.

The use of nonspecific PCR primers also contributes to low-quality sequence data because they tend to generate spurious amplicons along with the legitimate product. It would have been possible to design robust species-specific primers at the time of the studies by Brenton-Rule et al. [21] and Eloff et al. [9] because the first mitogenome for V. germanica was deposited on GenBank in May 2015 [35]. Our work contribute four new mitogenomes (South Africa, Argentina, and New Zealand) to the publicly available data, but this dataset represents a small fraction of the total diversity of V. germanica. Therefore, global surveys of V. germanica haplotypes may require the use of other mitochondrial markers, as the polymorphic regions used in this study may not be the most informative in other populations found across the globe.

Based on the geographic distribution of the two haplotypes found in our study, the data suggest that Hap2 was the original haplotype that was first recorded in 1972 in South Africa [16]. Hap2 is the dominant haplotype around the original invaded range, in sharp contrast to the rest of the species’ distribution, where Hap1 is dominant. Furthermore, it is likely that the arrival of Hap1 in South Africa launched a second invasion centered in the Boland region of the Western Cape, given the invasion history of a sudden expansion around this area in 2000 [18]. The low co-occurrence of Hap1 and Hap2 also suggests that the movement of different populations between the two regions in the Western Cape is limited, and that Hap1 is the more invasive haplotype and is responsible for recent invasive spread of V. germanica to the north, east, and south [19]. Veldtman et al. [19] showed that environmental variables related to moisture and temperature stress are responsible for the slow spread and geographic constraints on V. germanica in South Africa. Therefore, we hypothesize that Hap1 arrived after Hap2 and originated from a native population that may be more suited to conditions in the Western Cape, whereas Hap2 would have arrived first, but seems to have originated from a native population climatically less similar to this invaded region. This proposed scenario is supported by data in Eloff et al. [9] that found Hap1 to match a specimen from Austria and Hap2 to group with samples from Belgium, France, and Portugal, although exact haplotype matches were not found and further sampling in the native range is required to draw definite conclusions.

5. Conclusions

In light of the current paucity of haplotype data for V. germanica, we recommend that future assessments of phylogeographic patterns of this widespread invasive species include: (1) comprehensive sampling across the native and invaded ranges, including North America and Asia; (2) sequencing of complete mitogenomes of the most frequent lineages in each geographic area for selecting polymorphic regions that offer optimal haplotype resolution; and (3) the deposit of sequences and specimen collection data on public repositories (e.g., GenBank) for availability to other researchers. Sequencing of complete mitogenomes is presently inexpensive, and user-friendly software platforms for assembly, annotation, and subsequent analyses are widely available. Therefore, informative mitochondrial markers can be easily developed and used for surveys in other wasp species, such as Vespula vulgaris, Polistes dominula, and Vespa velutina, in their native and invasive ranges.