Multiple invasions of a generalist herbivore—Secondary contact between two divergent lineages of Nezara viridula Linnaeus in Australia

Abstract The presence of distinct evolutionary lineages within herbivorous pest insect taxa requires close attention. Scientific understanding, biosecurity planning and practice, and pest management decision‐making each suffer when such situations remain poorly understood. The pest bug Nezara viridula Linnaeus has been recorded from numerous host plants and has two globally distributed mitochondrial (mtDNA) lineages. These mtDNA lineages co‐occur in few locations globally, and the consequences of their divergence and recent secondary contact have not been assessed. We present evidence that both mtDNA lineages of N. viridula are present in Australia and their haplotype groups have a mostly separate distribution from one another. The north‐western population has only Asian mtDNA haplotypes, and the population with an eastern distribution is characterized mostly by European mtDNA haplotypes. Haplotypes of both lineages were detected together at only one site in the north of eastern Australia, and microsatellite data indicate that this secondary contact has resulted in mating across the lineages. Admixture and the movement of mtDNA haplotypes outside of this limited area of overlap has not, however, been extensive. Some degree of mating incompatibility or differences in the climatic requirements and tolerances of the two lineages, and perhaps a combination of these influences, might limit introgression and the movement of individuals, but this needs to be tested. This work provides the foundation for further ecological investigation of the lineages of N. viridula, particularly the consequences of admixture on the ecology of this widespread pest. We propose that for now, the Asian and European lineages of N. viridula would best be investigated as subspecies, so that “pure” and admixed populations of this bug can each be considered directly with respect to management and research priorities.


| INTRODUC TI ON
The association of insect pests with agriculture has facilitated the invasion of many of them outside their native range (Anderson, Tay, McGaughran, Gordon, & Walsh, 2016;Cristescu, 2015;Paini et al., 2016). Generalist host relationships, in particular, may facilitate invasion because of the greater number of host plants suitable for such organisms. Interpreting the invasion history and ecological features that underpin the geographical expansion of widely distributed insect species with varied ecology can be difficult. Understanding how the populations and ecology of an invasive species relate to populations found elsewhere in its distribution, however, is crucial for setting appropriate management and research priorities, as well as biosecurity measures.
Invasive populations may be comprised of once-restricted geographic lineages that have come into secondary contact and may have experienced genetic admixture (Blackburn et al., 2011;Dlugosch, Anderson, Braasch, Cang, & Gillette, 2015;Garnas et al., 2016;Reil et al., 2018;Walther et al., 2009). Some invasive insects, once seen as one generalist species, have been found to represent a complex of cryptic species, with each one relatively more specialized in its host use (Dumas et al., 2015;Malka et al., 2018;Rafter, Hereward, & Walter, 2013). But genetic analysis has also revealed cryptic species complexes comprised of different generalist species (Gikonyo et al., 2017;Hereward, Hutchinson, McCulloch, Silva, & Walter, 2017;Vyskočilová, Seal, & Colvin, 2019). All of this implies that close attention should be paid when divergent evolutionary lineages have already been documented in a particular species. Quantifying the respective invasive distributions of such lineages, and their species status, is an important first step towards understanding the ecology of each and their further invasive potential and pest status.
Scientific understanding, biosecurity planning and management decisions may be compromised when multiple invasions are not clearly understood, especially when unresolved cryptic species complexes are involved (Walter, 2003). The practical implications may be substantial. For example, knowing that the Scirtothrips aurantii Faure (Thysanoptera) population (sensu lato) that established in Australia does not represent the economically damaging African citrus pest of that name, but a host specialist species within that complex (Rafter et al., 2013), means that quarantine precautions must continue to be strictly maintained. Insects that exhibit differences in host use across their invasive and native ranges, as seen in S. aurantii (Rafter et al., 2013), provide a clear signal that more than one species may be involved. However, when considerable overlap is evident in host use across evolutionary lineages (Hereward et al., 2017;Vyskočilová et al., 2019), detecting the presence of cryptic species becomes more difficult (Paterson, 1991;Walter, 2003). This is particularly relevant to generalist insect herbivores because the diversity of their host interactions tends to obscure any differences across species.
Multiple invasions of such generalists are therefore more likely to remain undetected and uninvestigated, and so our understanding of these situations, and how to deal with them in practical terms, is likely to remain limited.
The biology and ecology of N. viridula, a globally distributed generalist bug, has been investigated for decades. The broad host plant range of N. viridula and the association of these insects with agriculture help explain its almost global distribution (Panizzi, 1997;Todd, 1989). However, three mitochondrial (mtDNA) lineages of N. viridula are known and their global distributions are mostly separate (Kavar, Pavlovčič, Sušnik, Meglič, & Virant-Doberlet, 2006;Li et al., 2010).
The European and Asian mtDNA lineages of N. viridula are distributed relatively widely but are reported to be sympatric only in Japan (Kavar et al., 2006). A highly divergent African lineage is restricted to Africa, but this lineage has been characterized genetically from only a single individual (Kavar et al., 2006;Li et al., 2010) and almost certainly represents a species distinct from the Asian and European lineages of N. viridula (Li et al., 2010). The consequences of secondary contact between the Asian and European lineages has not been investigated comprehensively.
Without investigating the biological characteristics mentioned above in relation to the mitochondrial lineages known to exist in N. viridula, this issue will remain unresolved. Such studies need to be conducted with a sound understanding of the genetic background of local N. viridula populations, particularly in countries where both mtDNA lineages are found as any admixture between the lineages will further complicate these investigations.
Nezara viridula feeds on plants from over 30 families, including many species of agricultural or horticultural significance, and all life stages of this insect feed preferentially on the fruit or seeds of their host plants (Todd, 1989). The bugs move through a sequence of host plants annually, as those hosts become more or less available and suitable for feeding and oviposition (Olson, Ruberson, Zeilinger, & Andow, 2011;Panizzi, 1997;Panizzi & Meneguim, 1989;Panizzi, Vivan, Corrêa-Ferreira, & Foerster, 1996;Todd, 1989;Velasco & Walter, 1992), with regional populations and different generations using only a small subset of the entire host plant range of the species.
At least one host plant species, Ricinus communis, does not appear to be equally suitable for N. viridula bugs of different provenance (Panizzi, 1997;Panizzi & Meneguim, 1989). Nezara viridula is invasive to Australia where it feeds mostly on crops and weeds (Velasco & Walter, 1992). In Australia, N. viridula has two principal generations (Velasco, Walter, & Harris, 1995), two fewer than is typical in the United States (Todd, 1989), and its abundance is frequently low during periods where few suitable host plants are available (Velasco & Walter, 1992;Velasco et al., 1995).
Nezara viridula has been present in Australia since at least 1911 (Clarke, 1992), and we present sampling and genetic analyses showing that both Asian and European mtDNA haplotypes are present on the continent. We investigated the distribution of the Asian and European N. viridula lineages in Australia by assessing the geographical distribution of their mtDNA haplotype groups. Microsatellites were then used to assess the amount of gene flow among individuals with mtDNA haplotypes of each lineage and among different geographic regions. We aimed to determine whether mating occurs between these lineages after secondary contact. We also conducted a global phylogeographic analysis of N. viridula using publicly available mtDNA sequences to compare Australian individuals with globally distributed samples of this species. Climatic factors affecting the distribution of each lineage in Australia are also explored. The implications for understanding and anticipating the outcome of multiple invasions and secondary contact in N. viridula are discussed. The results are then discussed within the context of the ecology and pest management of N. viridula, and their importance for generalist pest insects more broadly.

| Sampling and DNA extraction
A total of 649 individuals of N. viridula were sampled across Australia from 2014 to 2016, from 33 host plant species and across about 3,000 km (north to south) in eastern Australia, and also in two locations over 400 km apart in the north-west of the continent. The colour morph (Kiritani, 1970)  Note: The population site codes represent the specific sites of collections, and they are grouped together into particular populations for analysis (see text for explanation and Table S1 for actual site details). The geographical region labels include state and territory abbreviations as follows: QLD, Queensland and NSW, New South Wales. Populations marked with an asterisk (*) were excluded from all analyses that did not group populations into regions because of small sample size. Latitude and longitude here represent the mid-point between the sample sites that make up each designated population. The designated populations within geographical regions were at most 30 km from one another.
after hatching from the egg mass (Todd, 1989), so adult insects were used to avoid introducing bias. When the sample size was small for a site and host combination, however, a maximum of one 4th or 5th instar nymph was included to bolster sample size.

| Gene sequencing and analyses
The cytochrome C oxidase subunit I (COI) mitochondrial gene region was PCR-amplified using primers LCO1490 and HCO21980 (Folmer, Black, Hoeh, Lutz, & Vrijenhoek, 1994) and sequenced to determine whether the Australian insects belonged to the Asian or European mtDNA haplotype groups of N. viridula. Mitochondrial sequences used in previous phylogenetic analyses (Kavar et al., 2006;Li et al., 2010) and from various published and unpublished sources were obtained as of 18 April 2018 (Table S2). Errors in the published COI sequences were assessed visually after aligning all sequences, and poor-quality sequences were excluded where they could not be trimmed adequately. All sequences were aligned and trimmed for quality.
A phylogenetic analysis was performed on unique sequences from a 557-bp fragment of the more usual COI barcoding region for 519 sequences, which included trimmed sequences from this study.
A comparison was also made between sequences from this study and shorter sequences from the Kavar et al. (2006) data set to be certain no African lineage haplotypes were present in our sequences.
These analyses were used to determine the relationship between Australian N. viridula and those found elsewhere, and not to supersede more comprehensive phylogenetic analyses that have been performed by other authors using more genetic data (Kavar et al., 2006;Li et al., 2010). To compare evolutionary models, Jmodeltest 2.1.4 (Darriba, Taboada, Doallo, & Posada, 2012;Guindon & Gascuel, 2003) was used together with corrected Akaike information criteria (AICc). The comparison showed the HKY (Hasegawa, Kishino, & Yano, 1985) evolutionary model with a Gamma (+G) distribution was the most appropriate under all likelihood criteria so it was used for the phylogenetic analysis. Maximum-likelihood phylogenetic analysis was performed using the PhyML (Guindon et al., 2010) plugin as implemented in Geneious 9.0.5. Branch support was estimated using 5,000 bootstraps. Only sequences that were unique within each data set were used, and the COI sequence of Eurydema gebleri (Pentatomidae) (GenBank accession: KP207595.1) was included as the outgroup for the analysis.
Two nuclear DNA (nDNA) gene regions were also sequenced as previous studies used only mtDNA gene regions, and the evolutionary history of mtDNA is frequently not the same as nDNA (Toews & Brelsford, 2012). The tubulin alpha 1 (Tubα1) and the elongation factor 1 alpha (EF1α) gene region were sequenced for up to six Australian N. viridula from each host plant and site combination. Primers for both genes were developed de novo using Primer3 (Untergasser et al., 2012) based on the few sequences that amplified successfully using previously published primers (EF1α: Shirley and Prowler (Cho et al., 1995); Tubα1: TH_TubA forward and reverse (Buckman, Mound, & Whiting, 2013)). The presence of a microsatellite within the targeted EF1α gene region caused quality issues for many of the reverse primer sequences, so only the forward primer sequences could be used for many individuals (Table S1). To estimate the frequency of the nDNA haplotypes from diploid sequences,

| Microsatellite development and analyses
Microsatellite loci were used to assess the population genetic struc-  Table 1) were grouped into aggregate populations (i.e. regardless of host plant) for microsatellite analysis. These aggregate populations represent individuals collected at the same general location and time (Table 1).
The 16 aggregate populations that had 19 or more genotyped individuals (see Table 1) were used to assess the suitability and quality of the microsatellite markers using summary statistics (  (Earl & vonHoldt, 2012) and CLUMPAK (Kopelman, Mayzel, Jakobsson, Rosenberg, & Mayrose, 2015). Population genetic structure could thus be assessed independent of any a priori population designations. STRUCTURE analyses were performed with the "admixture" and with the "correlated alleles" model. Values of K (assumed number of populations) from 1 to 10 were tested with 10 replicates for each K value. A burn-in of 100,000 iterations, followed by 2,000,000 iterations, was used. The most appropriate K values were estimated using the Evanno method (Evanno, Regnaut, & Goudet, 2005 (Raymond & Rousset, 1995;Rousset, 2008), with and without high null allele loci (~10%), to test whether these loci influenced the analyses unduly. Differences in allelic richness across the two gene pools that were identified by STRUCTURE were calculated manually, to compare their genetic diversity. Pairwise F ST was used in Mantel tests to estimate isolation by distance (IBD) using the mantel.
randtest function in the ade4 package 1.7 (Dray & Dufour, 2007) in R, with 1,000 permutations and using F ST /1-F ST as the measure of genetic distance and the log of geographic distance as per Rousset (1997). The IBD analyses were performed with (a) all populations, and (b) with only populations from eastern regions (excluding those from north-western Australia). Eastern Australian populations were analysed separately because the sites are spatially more continuous, and the north-western region is a geographic outlier. When two populations had been collected from the same sample location at different times (Table 1), the population sampled at the earliest time was used and so 89 individuals from populations EMRb, NARb and GRIb were excluded from IBD analyses.

| Modelling climatic suitability
Climatic limits to the distribution of the lineages of N. viridula were modelled using MaxEnt version 3.4.1 (Phillips, Dudík, & Schapire, 2020) as implemented in BCCVL (Hallgren et al., 2016). MaxEnt was chosen because it predicts the largest possible range size from the available data (Merow, Smith, & Silander, 2013 (Merow et al., 2013;Saupe et al., 2012). Receiver operating characteristic (ROC) curve was used to assess the suitability of the model.    Figure 1) showed a change from significant to nonsignificant differences after these loci had been excluded (Tables 3 and Table S5).

| Gene flow across mtDNA lineages in Australian Nezara viridula
For the STRUCTURE analysis of the microsatellite data, the most likely inferred value of K was 2 (see Figure S2 for the likelihood of K and ΔK). The eastern populations of N. viridula cluster together in this K = 2 analysis, and are separate from those populations further west in northern Australia (from Kununurra and Darwin). The analysis, most importantly, shows the discordance between the microsatellite and mtDNA data in eastern Australia, with gene flow between individuals F I G U R E 1 (a) Map of northern and eastern Australia with Nezara viridula sampling sites represented by black dots. The circles that surround sets of black dots amalgamate these nearby sampling sites into labelled regions (Table 1). Superimposed on the map is a haplotype network showing the relationship between the COI haplotypes of 480 Australian N. viridula individuals (all sequences), and also a K = 2 STRUCTURE analysis of the microsatellite data from 571 individuals. Horizontal bars above these locality-specific inset boxes indicate the mitochondrial haplotype of the corresponding individual in the STRUCTURE analysis boxes below the bar. Those individuals that were genotyped, but which did not have their COI gene sequenced, are represented above the boxes by a thin black line. Bowen is the only location in which bugs with both mitochondrial haplotypes were found together, with a lighter shade of either orange or blue to distinguish bugs from that location from others. (b) A K = 3 STRUCTURE analysis is also shown below the map. (c) A global inset map is used to show the placement of all publicly available N. viridula sequences ( Figure S1) with mtDNA belonging to both lineages (Figures 1 and 3). The private allele analysis shows the same pattern (Table S6). The north-western region had 15 private alleles from 185 individuals when compared to one private allele in a single individual for northern QLD and no private alleles for southern NSW. Other regions were excluded from the analysis so that an equivalent spatial could be made. Although eastern N. viridula individuals represent a single interbreeding population (STRUCTURE analysis in Figure 1), in the K = 2 STRCUTURE analysis, individuals from northern Queensland and central Queensland were more frequently assigned to the same population as north-western F I G U R E 2 Haplotype networks for the EF1α (440bp) and Tubα1 (423bp) nuclear genes sequenced for individuals of Nezara viridula collected in Australia. Haplotype shade (ordered stepwise geographically from black (north) to white (south)) indicates the number of haplotypes obtained from each region (Table 1). For each region that contributes to the number of bugs with a particular EF1α or Tubα1 haplotype, the colour (blue or orange) of the circle around that haplotype corresponds to the individuals associated with each mitochondrial lineage (see key). Each individual is represented by two haplotypes as estimated by DnaSP 5.1 (Librado & Rozas, 2009). EF1α and Tubα1 haplotypes associated with both mtDNA lineages in the same region are indicated with a lighter shade of either orange or blue to distinguish th`em from those with a single regional mtDNA lineage  (Table 1). The results of pairwise exact G tests for genotypic differentiation are shown above the diagonal, with levels of significance: *0.01-0.05, **0.01-0.001, ***≥0.001. If more than one nearby location is included in a single population, the average pairwise geographical distance was used.
individuals than were individuals from southern Queensland and New South Wales (Figure 1). This was the same for the K = 3 analysis, but fewer eastern individuals were assigned to the same population as north-western individuals ( Figure 1). All 12 loci were included in this analysis as the presence of null alleles should have only a minor impact on assignment tests (Carlsson, 2008).  (Kiritani, 1970) of N. viridula were found in the north-western Australian populations, while only the G colour morph (Kiritani, 1970) was found in eastern Australia.

| Climatic suitability
The north-western region, where only bugs with Asian mtDNA are found, was predicted to be unsuitable for individuals with European lineage mtDNA in the warmest annual quarter (Figure 5a). Much of F I G U R E 3 A principal components analysis (PCA) using microsatellite genotype data from the 10 loci without high null allele estimates, and from all Nezara viridula populations from Australia. Individuals are coloured according to: (a) their region of origin and (b) in relation to their mitochondrial haplotype group. Regional abbreviations are: NSW, New South Wales; NT, Northern Territory; QLD, Queensland; WA, Western Australia  Figure S5). Mean temperature for the warmest quarter was the best predictor for climatic suitability of N. viridula with the Asian mtDNA haplotype and mean humidity contributed little ( Figure S6).

N. viridula populations
Nezara viridula populations in Australia, as well as globally, mostly belong to the European mtDNA haplotype group (Figure 1 and Figure S1). The mitochondrial data presented here reveal that the Asian and European mtDNA lineages of N. viridula are known to cooccur only in Australia and Japan, though they may also be found together elsewhere. Secondary contact between the Asian and European lineages of N. viridula in Australia has resulted in crossmating and the production of offspring, as evidenced most clearly by the discordance between the mtDNA and nuclear (microsatellite) data from insects collected at Bowen (Figure 1). Mating between the Asian and European lineages in Australia has, however, not resulted in complete admixture or the widespread distribution of both mtDNA haplotype groups, as has been found other invasive insects in Australia (Toon et al., 2016). Importantly, the admixture across  (Tables S6 and   S7). Individuals from southern NSW and northern QLD, by contrast, had almost identical allelic diversity (Tables S6 and S7). Three of the four colour morphs found in Australian N. viridula are associated only with bugs from the north-western regions, providing further evidence for limited movement of insects between these regions. Some amount of mating probably occurs between these two populations, but further tests are required to determine whether mating between the lineages occurs at a lower frequency that within-lineage mating.
Where mtDNA haplotype groups of both lineages are found together in northern Queensland, a bias towards the European lineage in the microsatellite data is clear, even when mtDNA haplotypes of the Asian lineage are more common (Figures 1 and   3), and even when geographic distance has been accounted for ( Figure 4). Asian lineage mtDNA haplotypes are not found south of the northern QLD region, despite mating occurring between bugs of both lineages in sympatry, and despite N. viridula being a capable flier (Gu & Walter, 1989;Kester & Smith, 1984;Kiritani & Sasaba, 1969). The lack of Asian mtDNA haplotypes further south must be explained, as there are no obvious impediments to movement between the regions in which bugs of both mtDNA haplotype groups co-occur and regions further south which have only European mtDNA haplotypes (northern QLD and central QLD, respectively; Figure 1). The genetic resolution of the nuclear gene fragments was low overall and because of this generally uninformative. The most common haplotypes are likely shared because of incomplete lineage sorting given that few substitutions separate them, although admixture between the lineages might also have play a role.
Biogeographical barriers to gene flow such as regional host availability, and possible ecological differences across the Asian and  Walter, 1992), and these are likely to be scarce in dry regions, including the gaps between sampled regions, and so not sustain N. viridula populations year-round. Host availability may be higher in wet years and facilitate the movement of individuals between populations that are regularly isolated from one another by dry conditions. Host plant distribution does not explain completely why the distribution of each mtDNA haplotype group is restricted in eastern Australia. Some degree mating incompatibility (Jeraj & Walter, 1998;Ryan et al., 1996), population ephemerality or sex biased dispersal or might also play a role. Experimental tests to investigate each of these possibilities need to be made directly.

| Ecological differentiation in Nezara viridula and speciation in generalist insects
The extent to which the Asian and European lineages of N. viridula are ecologically differentiated with respect to one another is still unclear and requires experimental investigation. The genetics results we present provide the framework from which sound ecological tests can be designed and conducted in Australia. Some of the differences documented across N. viridula populations and across localities (Aldrich et al., 1987;Jeraj & Walter, 1998;Kon et al., 1988;Miklas et al., 2000;Panizzi & Meneguim, 1989;Ryan et al., 1995Ryan et al., , 1996Todd, 1989;Virant-Doberlet et al., 2000)  Australia (Jeraj & Walter, 1998;Ryan et al., 1996).   (Velasco & Walter, 1993 (Chanthy et al., 2015). Experiments investigating interactions between life-history parameters and climatic factors need to be conducted for primarily Asian lineage bugs from north-western Australia and for admixed populations in northern Queensland.
Comparison across N. viridula and other Nezara species offers further support for climatic differences across lineages within this genus.
Nezara antennata, the nearest relative of N. viridula, is adapted to cooler climates than N. viridula (Musolin, 2012;Tougou, Musolin, & Fujisaki, 2009;Yukawa et al., 2007), and given the similarities in host range and mating between these two species, adaptation to different climates may have been the common precursor to speciation in this group.
Climatic tolerances might be one important component of speciation for host-plant generalists more broadly (Brunner & Frey, 2010;Gikonyo et al., 2017;Hereward et al., 2017). In Australia, N. viridula from eastern regions experience much cooler conditions and more variable photoperiods than those in the north-western region (Table S8) (Musolin, Tougou, & Fujisaki, 2011;Tougou et al., 2009;Yukawa et al., 2007) and this pattern may, in part, be because European lineage bugs invaded that country.
With respect to the host plant associations of N. viridula in Australia, no host-associated genetic structure was evident within any of the regions sampled (STRUCTURE, Figure 1, and PCA,

| Implications for management
Two genetic populations of N. viridula are present in Australia ( Figure 1) that should be treated separately ecologically and with respect to pest management and research. However, populations of N. viridula from eastern regions of Australia can be treated as a single interbreeding population, even though the biology of individuals in northern regions may be affected more strongly by introgression between the Asian and European lineages (Figure 1 and Table 3). Any ecological differences between Asian and European lineage N. viridula will influence the local pest ecology of this bug in Australia and internationally, but thorough experimental and genetic investigation of N. viridula is required to establish whether there are any consistent ecological differences between them.
The possibility of climatic differences between the lineages has serious consequences for the geographic range of N. viridula globally.
It is unlikely to be coincidental that the two genetically different pop-

| General conclusions
Secondary contact and admixture between the lineages indicates that they do not represent cryptic species. However, the incomplete mixing of the two lineages in Australia, and the potential ecological differences across them suggest that they could be subspecies. Detecting cryptic diversity and multiple introductions of generalist insects is challenging because the features of their ecology that would normally be seen to signal their presence (such as interacting with novel plant species in the case of host plant specialists) can be difficult to discern. This point is illustrated by the situation that we have uncovered in N. viridula, where two divergent and previously allopatric lineages mate when in sympatry but are not completely admixed within their invasive range. Over time N. viridula in Australia may even become a single fully admixed population with traits from both lineages, but this has not happened yet despite probably over a century of both being in Australia (Clarke, 1992). Clarification of the behavioural and ecological process responsible for the genetic structure of Australian N. viridula is a priority for our understanding of this pest, and resolving this point will also contribute to our understanding of how independent lineages of generalist herbivores should be in general treated. In the meantime, treating the lineages of N. viridula independently of one another with respect to quarantine risk, pest management and ecological research seems the most prudent approach to dealing with the uncertainty still surrounding our understanding of this species.

ACK N OWLED G M ENTS
We are grateful to Tanya

CO N FLI C T O F I NTE R E S T
There are no competing interests to declare.

AUTH O R CO NTR I B UTI O N S
All authors designed the research. DRB conducted the sampling, designed the primers, performed the laboratory work and analysed the data. JPH performed the next-generation sequencing. All authors contributed to the interpretation and writing the paper.

DATA A N D M ATE R I A L S AVA I L A B I LIT Y
GenBank accession numbers were assigned to all sequences gen-  (1987). Pheromone strains of the cosmopolitan pest, Nezara viridula