Species delimitation in asexual insects of economic importance: The case of black scale (Parasaissetia nigra), a cosmopolitan parthenogenetic pest scale insect

Asexual lineages provide a challenge to species delimitation because species concepts either have little biological meaning for them or are arbitrary, since every individual is monophyletic and reproductively isolated from all other individuals. However, recognition and naming of asexual species is important to conservation and economic applications. Some scale insects are widespread and polyphagous pests of plants, and several species have been found to comprise cryptic species complexes. Parasaissetia nigra (Nietner, 1861) (Hemiptera: Coccidae) is a parthenogenetic, cosmopolitan and polyphagous pest that feeds on plant species from more than 80 families. Here, we implement multiple approaches to assess the species status of P. nigra, including coalescence-based analyses of mitochondrial and nuclear genes, and ecological niche modelling. Our results indicate that the sampled specimens of P. nigra should be considered to comprise at least two ecotypes (or "species") that are ecologically differentiated, particularly in relation to temperature and moisture. The presence of more than one ecotype under the current concept of P. nigra has implications for biosecurity because the geographic extent of each type is not fully known: some countries may currently have only one of the biotypes. Introduction of additional lineages could expand the geographic extent of damage by the pest in some countries.


Introduction
Delineation of asexual lineages remains a challenge for taxonomists. The biological species concept [1] and other approaches that require consideration of mating (e.g. specific mate recognition species concept; [2]), gene pools (e.g. genetic species concept; [3][4]) and independent evolutionary trajectories (e.g. evolutionary species concept; [5]) are not meaningful for a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 delimiting lineages in which every individual is reproductively isolated from every other and thus on its own evolutionary path. Other concepts, such as the phylogenetic species concept [6] and the mtDNA barcode concept [7], are also unsatisfactory for applying to asexual organisms because they require arbitrary cut-offs to determine the level at which a clade is considered a species. Whether there is a barcoding gap [8] or a coalescent point [9] is irrelevant for delimitation of asexual lineages, as neither have biological meaning. They simply represent patterns on phylogenies that could be the result of a plethora of causes, including extinction or under-sampling.
Some authors (e.g. [10]) have argued that the taxonomic rank we call "species" has no greater biological meaning than higher ranks such as "family" or "genus". "Species" might simply be groups of individuals that are more closely related to one another than they are to other such groups, and this is probably more the case in asexual organisms than in obligatorily sexual ones. Nevertheless, species remain the fundamental units in systematic biology and ecology, and are essential for communication and application of conservation strategies [11] and quarantine decisions [12]. Confusion over species identity can have serious negative impacts on pest management programs [13] and international trade [14]. For these reasons, it is just as important to define "species" in asexual taxa as it is for sexual taxa, even though the biological meaning might differ.
Some studies dealing with species delimitation in asexual taxa, such as bdelloid rotifers (Rotifera: Bdelloidea: Philodinavidae) (e.g. [15][16][17]) and darwinulid ostracods (Crustacea: Ostracoda) (e.g. [18]), have defined species on an ecological basis, i.e., if two or more lineages have sorted into different niches, then each lineage may be considered to be a distinct species. The authors of these studies suggested that the main drivers of speciation are the same in sexual and asexual organisms (population isolation and adaptive selection), but the boundaries between asexual species are formed by processes other than sexual reproduction (an idea shared by De Queiroz [19]).
Here we investigate species delimitation in Parasaissetia nigra (Nietner, 1861) (Hemiptera: Sternorrhyncha: Coccoidea: Coccidae), the "nigra scale" or "black scale insect", which is a polyphagous scale insect that feeds on plants from more than 80 families [20][21]. The species is thought to have originated in Africa [22] but is now widespread throughout the world. It is an important pest of ornamental and greenhouse plants [23], fruits and other plantations in tropical and subtropical countries [24][25][26][27], and is a moderate pest of several crops in temperate regions [22]. Parasaissetia nigra has also been identified as a potentially invasive species in areas of environmental concern, such as the Galápagos Islands [28]. Infestation by the insect can reduce the vigor of host plants by depleting their sap reserves, and the honeydew produced by the insects forms a suitable medium for the growth of sooty molds, which can lower the market value of fruits and ornamental plants [22].
It is generally accepted that females of P. nigra reproduce parthenogenetically [24, 29-31]: males have been reported twice [32][33], but there has been no evidence presented that show a direct connection between these specimens and P. nigra. Information that could confirm the existence of males in P. nigra (such as the existence of sperm bundles in the spermatheca or heterochromatic chromosomes in embryos, e.g. [34]) has not been documented in the species. If males do exist they must be rare: P. nigra is common, frequently collected and intensively studied worldwide, yet these two reports remain the only suggestion of the existence of males.
De Lotto [35][36] reported that there was considerable variation in the size, form and color of specimens of P. nigra from different areas of Africa, but could find no consistent morphological differences to warrant recognition of distinct species. Ben-Dov [30] examined specimens of P. nigra from across its range and reported differences in numbers of marginal setae, preopercular pores, submarginal tubercles and pocket-like sclerotisations, and different lengths of antennae and legs. Similarly, Hodgson [37] reported that adult females collected from Chad and Laos differed from those of other countries. Both Hodgson and Ben-Dov thought morphological differences might be linked to cryptic species diversity within P. nigra, but neither thought the evidence conclusive.
In this study, we aim to address the following questions: 1. Is P. nigra a single species with broad host-use and a cosmopolitan distribution, or should the clade be considered to represent more than one species?
2. If considered multiple species, are those species ecologically distinct, e.g., associated with specific host plants, geographic locations or climatic variables?
We use both single (COI) and multi-locus (18S, 28S, COI, Dynamin and EF-1α) sequence data and multiple coalescent approaches [38-39] to determine the number of genetic clusters within P. nigra that could be recognized as species under coalescence-based species delimitation [40]. Alternative delineation hypotheses were tested under coalescent approaches and recent speciation events were reconstructed by applying probabilistic models. However, because there is no (or little) possibility of genetic exchange between individuals of a strictly (or largely) parthenogenetic organism, other than from mother to daughter, we need additional criteria for determining species boundaries that are testable and have some meaning for end users of the taxonomy.
Mallet [41] argues that conspecific individuals have common gene pools (sensu Dobzhansky [3][4]), and treats species as separately identifiable genotypic clusters. However, reproductive isolation is not viewed as the main mechanism stopping the exchange of genetic material among clusters. Under this concept, species boundaries are maintained by selection, and different species should have different ecological and/or geographic distributions-isolating mechanisms that Van Valen [42] and Andersson [43] used when applying the ecological species concept to delineate plants. We use a modification of this approach, defining species within parthenogenetic lineages as genetic clusters that are ecologically differentiated from other such clusters. The possibility of this resulting in paraphyletic species is discussed. have allowed for an easier taxonomic interpretation. In the present study, 65 samples identified as P. nigra were included in molecular analyses. All were collected from outdoors and they represent populations from at least 50 different host plant species (25 families) and 44 localities across four continents and multiple islands ( Table 1). The four other described species of Parasaissetia could not be included as specimens were not available for DNA extraction, but these are all morphologically different and clearly distinguishable from P. nigra [22]. Four species of Saissetia Déplanche, including the type species S. coffeae (Walker), were chosen as outgroups because the genus is closely related to P. nigra [46]. We included C. longulus (Douglas) because it is also closely related to P. nigra [47][48]. A specimen of C. hesperidum L. was used to root the trees because it represents a different but closely related tribe, Coccini [49].

Taxon sampling and DNA extraction
Insects collected in the field were preserved in absolute ethanol (> 99.5%; p.a.) and then stored at 4˚C. Genomic DNA was extracted from young adult females or post-reproductive  Table 2). Three primer pairs were used for amplifying two regions of COI (Table 2). CI-J-2183 (Jerry) and C1-N-2568 (Ben) ( Table 2) was used for a 3' region of COI with a step-down PCR program ( Table 2). The barcode region of COI was amplified with primer pair PcoF1 and HCO for the six outgroup taxa but, because this combination did not work well for P. nigra, HCO was replaced with newly designed primer (nigra_Ben) ( Table 2) modified from C1-N-2568. A negative control was used for all PCR reactions and each 25 μL PCR mixture comprised 5 μL 5x PCR buffer, 2 μL dNTP (2mM), 1.  Species delimitation in Parasaissetia nigra PCR products were cleaned using Exonuclease I and Antarctic Phosphatase (cat. no. M0293S and M0289S, New England BioLabs, Australia), or with the ammonium acetate/ethanol precipitation method used by Lin et al. [47]. In amplicons with multiple bands, the target band was excised from a 1% agarose gel under UV illumination and purified using the Wizard SV Gel and PCR Clean-up System (cat. no. #A9281, Promega, Madison, USA) following the manufacturer's instructions. All PCR products were sequenced using Sanger sequencing at Macrogen Inc. (Republic of Korea).

Sequence editing and alignment
Sequences were edited using MEGA5 [63], then imported and aligned visually in Se-Al v.2.0 [64]. Translations to amino acids were used to assist alignment of the three protein coding regions (COI, Dynamin and EF-1α) and to check for stop codons. Intron-exon boundaries of Dynamin and EF-1α were assigned using the GT-AG rule [65] and comparison with annotated sequences in GenBank.

Identification of clades/lineages
We used phylogenetic methods to identify clades common to individual gene trees and analyses of concatenated data that could represent putative species, and thus form the basis for additional analyses. Most methods of phylogeny estimation assume that base composition among taxa is homogeneous [66] so we firstly checked for base frequency bias among taxa (non-stationarity) in all datasets using PAUP Ã 4.0b10 [67]. Each codon position for the three proteincoding regions was tested independently, with and without invariant sites.
We used Bayesian inference (BI) implemented in MrBayes v.3.2.1 [68] to analyze a concatenated dataset and each gene region separately, except that we grouped the nrDNA loci 18S and 28S given that they are physically linked in the ribosomal array. Outgroups were used to root the phylogenies and introns were excluded because they could not be unambiguously aligned across the distantly related taxa.
Substitution models for each partition were selected using MrModeltest 2.3 [69] ( Table 3). The COI, Dynamin and EF-1α datasets were each assigned two partitions: first plus second codon positions, and third codon positions. The GTR [70] + I + G model (nst = 6, rates = invgamma) was chosen for each partition of COI, EF-1α and 28S. A K2P [71] + I + G model was applied to 18S, and a JC69 [72] model and K2P + I + G model was selected for the two partitions of Dynamin. Each analysis comprised two independent runs of 90 million (nrDNA), 60 million (COI, Dynamin and EF-1α) or 40 million (concatenated dataset) generations with the default setting of four Markov chains (three heated and one cold), temperature = 0.10, starting from a random tree and sampling every 1000th generation.
The performance of each Bayesian analysis was checked by examining the average standard deviation of split frequencies (should be less than 0.01) [73], PSRF values (should be close to 1.00) [68], the absolute value of the difference between the harmonic means of the two runs Table 3. DNA substitution models, the length of runs and the number of burn-in used for each gene region in BI analyses. (should be less than 2) [74] and the ESS (Effective Sample Size, should be > 200) [75] of each statistic determined using Tracer v.1.6 [76]. The number of trees discarded as the burn-in period varied with each analysis, depending on when stationarity was reached (60 (18S + 28S), 50 (COI), 24 (Dynamin), 20 (EF-1α) and 30 (concatenated) million generations respectively.). A maximum clade credibility (MCC) tree with posterior probability values from the two runs of each analysis was selected using TreeAnnotator [77] and the post-burn-in trees.
As another check of clade membership, we also analyzed each dataset using maximum parsimony (MP) using PAUP Ã 4.0b10 [67] (S1 Appendix) because it has different underlying assumptions from BI and congruence among them indicates results are not sensitive to the method of phylogeny estimation.
Because introns and hypervariable regions of rDNA had to be excluded in the outgrouprooted analyses, we used BEAST 1.8.0 [77] to analyze a concatenated dataset that included only sequences from P. nigra but which included these regions. Outgroups are unnecessary in BEAST because this software roots trees by applying a (calibrated or implicit) molecular clock [75].
XML files for BEAST were generated using BEAUTI 1.8.0. [77]. Partitioning was the same as for the MrBayes analyses except for the addition of intron partitions for Dynamin and EF-1α. A HKY [78] + I + G was used for all partitions except the exons of Dynamin, for which a K2P + I model was applied.
A gamma distribution with initial value, shape and scale set to 1.0, 0.001 and 1000 was used as the prior for these analyses. The exon and intron regions of Dynamin and EF-1α were treated separately and each had the same, but unlinked, clock model. The tree prior was set to the pure birth Yule speciation process [79], which assumes that all lineages speciate at the same rate in any given time period and have uniform birth rates (prior distribution low = 0.0; upper = 1.0E100; initial = 1.0). This is the simplest branching model in species-level processes [80].
Each run comprised 50 million generations starting from a random tree and was sampled every 1000th generation. The performance of each BEAST run was checked by examining the ESS values of each statistic shown by Tracer v.1.6. A maximum clade credibility tree was chosen from the posterior set using TreeAnnotator and after removing samples from the burn-in period.

Multi-locus coalescent species delimitation (*BEAST)
Using the same partitioning schemes and models as those applied for BEAST, we tested multiple species hypotheses based on supported clades recovered in analyses of individual and concatenated datasets that had a COI divergence >2%. Eight species hypotheses, ranging from two to six species, were tested (Table 4) by assigning terminals to different taxa and comparing the marginal likelihood values of the different hypotheses after the species tree was generated.
Each run started from a random tree and was sampled every 5000 generations. We implemented two models of population size, "piecewise linear and constant root" (PLCR) and "piecewise linear" (PL), rather than the piecewise constant population size model because P. nigra is extremely widespread [20] and is likely to have undergone a population expansion. Each hypothesis was run under both Yule [79] and birth-death [81] speciation processes separately. Convergence of each Ã BEAST run was assessed in the same manner as for BEAST. The number of generations/run (200 to 600 million) and the burn-in (5 to 200 million) required for stationarity varied for each analysis.
The harmonic means [82] of runs for each species hypothesis were estimated in BEAST 1.8.0. The favored hypothesis was that with the lowest harmonic mean that was significantly different from other means, with significant difference assessed using Bayes Factors [74]. Because the harmonic mean method might not be sufficiently sensitive [83], we also used a posterior simulation-based approach (Akaike's information criterion through Markov chain Monte Carlo (AICM) comparisons [84]) to identify the best species hypothesis. The hypothesis with the lowest AICM was assumed to be best, with the P values from calculating the exponential value ((minimum AICM-maximum AICM)/2) indicating whether any two hypotheses were significantly different from each other (P < 0.05, d.f. = 1).

Single-locus coalescent species delimitation (GMYC)
We used the single-threshold General Mixed Yule Coalescent (GMYC) method [85] on the COI sequences of P. nigra. GMYC estimates the number of clusters ("species") by recognizing the transitions from between-to within-species branching patterns on an ultrametric and fully dichotomous tree. A likelihood ratio (LR) test was applied to determine whether the mixed coalescent and Yule stochastic lineage growth (GMYC) model provided a better fit to the data than the null model, which hypothesizes that the entire sample belongs to a single species with no shift in branching processes. The ultrametric tree needed as input was estimated using BEAST with a strict clock, the substitution rate set to 1.0 (as suggested by Gamble et al. [86]), and two partitions (first + second codon positions: third codon position) with an HKY+I+G substitution model for each. The GMYC analysis was implemented using the APE [87] and SPLITS [88] packages in R [89].

Environmental niche modelling
To determine whether there is evidence of ecological differentiation among the identified lineages within P. nigra, which could be interpreted to indicate the presence of distinct ecological species, we compared the geographic distribution and host-use of each clade.
Poor representation of most lineages precluded a thorough assessment of geographical and ecological limits to distribution, but ANZ (from 6 localities) and W3 (23 localities) lineages have wide distributions and relatively good sampling from Australia. This allowed a comparison of the lineages across a continental-scale rainfall and temperature gradient. The two collecting localities of P. nigra (W3 lineage) from Rakimov et al. [48] were included in our species distribution modelling to increase sample size. We used the geocode for each collection and initially modelled the distribution of each lineage in MaxEnt v. 3.3.3k [90] using the "Bioclim" dataset [91] of 19 high resolution (c. 1 km at the equator) layers (S2 Appendix). We then used the six top performing variables (those contributing ! 10% to the model) (S2 Appendix) in an identity test run using ENMTools (perl script version 1.0) [92]. We limited the area considered to a background to the geographic extent determined in the original full Bioclim MaxEnt model: therefore, the question addressed by the identity test was "In eastern and southern Australia, do ANZ and W3 occupy environments that are less similar than expected by chance?".

Results
Sequence data for all five DNA regions were obtained for all specimens except that Dynamin could not be amplified for three (S3 Appendix). GenBank accession numbers of sequences are given in S3 Appendix. No base composition bias (non-stationarity) was detected among taxa in any of the datasets (P = 1.00 in all). Stationarity and convergence between runs was reached in all Bayesian analyses.

Identification of lineages/clades
There was strong support for the sampled specimens of P. nigra forming a clade in all analyses  (Fig 1). The COI uncorrected genetic distances (pdistances) within the six clades ranged from 0.0-0.9%, and between the clades ranged from 3.0-10.0% (Table 5).

Tests of species hypotheses using *BEAST
Because clades ANZ and W1 were well supported in the BEAST-generated chronograms of all datasets, the two clades were treated as distinct species except in two hypotheses (Sp2 and 3Sp2) that combined the two clades into a single species (Table 4). Because sequences of the two specimens from Ghana (YPL00238 and YPL00239) did not cluster together in BEAST analyses of EF-1α (S4 Fig) or rDNA (S5 Fig), five alternative scenarios for their relationships were developed (3Sp1, 4Sp, 5Sp1, 5Sp2 and 5Sp3) ( Table 4). Finally, the 6Sp he hypothesis treated the six lineages as different species.
Except for hypothesis 5Sp3, which had low ESS (< 50) in all priors linked to EF-1α even after 600 million generations, all Ã BEAST runs reached stationarity. With the exception of the 2Sp, 3Sp2 and 5Sp2 hypotheses, there were no significant differences between results using the PLCR and PL population models under both the Yule and Birth-death speciation processes, as assessed by AICM (P < 0.05) and harmonic means (absolute value of difference > 6) ( Table 6). Both Bayes Factors and AICM indicated that the 6Sp hypothesis was the best, regardless of speciation and population size models (Table 6), but this model was only significantly better than the second best hypothesis (5Sp2) in the AICM from the combination of Birth-Death and PL models ( Table 6). Species delimitation using GMYC Six ML entities ("species") (confidence interval: 6-6) (S6 Fig) were recognized after running the single threshold GMYC model (likelihood ratio = 31.03, P < 0.01). These six ML entities correspond to the six major clades found in the BEAST analyses of COI, Dynamin and concatenated datasets (Fig 1; S2 and S3 Figs), which also formed the basis of the 6Sp hypothesis in Ã BEAST analyses.

Ecological differentiation
None of the lineages of P. nigra were clearly host specific and some hosts were shared across lineages (Table 1). For example, Cordia dichotoma was host to YPL00361 (clade W1) and YPL00426 (clade W2), and Plumeria obtusa (YPL00492 and YPL00085) and Psidium guajava (YPL00495 and YPL00522) were host to some specimens from clades W2 and W3.
Clades ANZ and W3 occupy different climatic zones in Australia (Identity test: Shoener's D and I statistic, 0.03 < P 0.04) (Fig 2). Temperature in the warmest period (bio 5) was the most important variable in the model for ANZ, whereas precipitation coldest quarter, precipitation warmest quarter, and annual mean temperature (bio19, 18 and 1 respectively) were most important for the model for W3. In general, ANZ is distributed in temperate areas Species delimitation in Parasaissetia nigra (Fig 3A), whereas W3 is mostly restricted to coastal subtropical (eastern areas) and Mediterranean (regions around Adelaide and Perth) areas in Australia (Fig 3B).

Discussion
Overall, we found at six deeply divergent lineages (>3% COI divergence; defined in Fig 1) within P. nigra that do not represent host-specific or geographically distinct clades. In sexual organisms, these divergent lineages could be considered distinct species under biological and evolutionary genetic species concepts. Given that they are parthenogenetic, if all lineages within P. nigra were equivalent in biology and ecology, recognizing only a single species could be justified given there is also no clear morphological differentiation. Here, however, we find that two  Species delimitation in Parasaissetia nigra lineages occurring in Australia and elsewhere appear to be ecologically distinct, occupying climatically different regions. We argue that ecological differentiation warrants the recognition of two species among these parthenogenetic lineages (ANZ and the rest). Below we discuss the practicalities of this, such as needing to recognize paraphyletic species-those lineages that have the same morphology and ecology but which do not form a monophyletic group.

Species delimitation for biosecurity
How species are delimited can affect biosecurity, particularly in the areas of international and intra-national quarantine and trade. For example, products infested with invasive alien species are prohibited entry to some countries-the reason agricultural and other goods are inspected at state and national borders (e.g. [12]). If an invasive pest species is known to already occur in a region, products carrying that species might be allowed, whereas they might be prohibited or destroyed if the pest is not known to already occur there. Taxonomic "over-splitting" could be an unnecessary hindrance to trade, whereas "under-splitting" or "over-lumping" might lead to greater spread of invasive alien species and pose a threat to global agriculture. For asexual lineages, the delimitation of species is far more arbitrary than it is for obligatorily sexual ones -a criterion of no or little gene flow cannot be applied since every individual would be recognized as a species, and each species would have a short life span (several months in the case of P. nigra). Application of so-called "objective" methods (e.g. [93]), such as multi-locus coalescence, has been proposed for use in delimiting asexual lineages (e.g. 4x rule, [16]). Coalescence is strongly influenced by the effective population size (Ne), time, mutation rate and population structure of the target organisms (e.g. [94]) such that, if populations are large and exchange genes (migrants), coalescence will be deep. In coalescence-based methods of species delimitation using multi-locus data on sexual species, the number of species can be overestimated because population structure, rather than species boundaries, might be recovered [95]. Also, uneven sampling from across a species' range, such as we have for P. nigra, can also lead to coalescence that represents population structuring other than species boundaries (as in many barcoding studies reviewed by Lohse [96]).
Time to coalescence increases dramatically with reduction of gene flow or with increasing number of populations [97], and therefore coalescence is expected to be much deeper for asexual lineages (which have no gene flow and all individuals are "populations"). This means that depth of coalescence in asexual lineages is much deeper than that in sexual species, even when the number of individuals is the same. There is no "population structure" in asexual lineages, in the sense of differential likelihood of genes being exchanged, and no "speciation" in the sense of genetic isolation as a consequence of cessation of gene flow, thus coalescence is largely driven by mutation rate, Ne and extinction. Consequently, we argue that coalescence-based species delimitation methods, despite their occasional use for asexual organisms (e.g. [17,98]), have little useful meaning for delimiting parthenogenetic species for biosecurity other than to apply the same method of delimitation as those applied to sexual species. Indeed, given that phylogeny alone provides an explanation for the pattern of traits in parthenogenetic organisms, any division of parthenogenetic lineages into species is arbitrary [99].
For parthenogenetic organisms of potential biosecurity concern, we suggest that the priority should be to identify ecologically distinct lineages as species, while minimizing any arbitrary species cut-offs, i.e. walking a tightrope between minimizing potential environmental or agricultural harm and minimizing potential disruption to trade. Morphological distinctions might sometimes be considered surrogates for ecological differences.
Is P. nigra a cryptic species complex?
The current concept of P. nigra is that of a phenotypic cluster (sensu Sokal & Crovello [100])a group of individuals that are morphologically similar to each other but phenotypically different from all other scale insects. If P. nigra was largely sexual, reciprocal monophyletic across mitochondrial and multiple nuclear genes (Fig 1; S2-S5 Figs), combined with the coalescence results found here (GYMC and Ã BEAST), would indicate that there are six species present within the currently recognized species. There would probably be little dissent from scale insect taxonomists if we were to recognize and describe five additional species because the levels of genetic differentiation are high and the generally accepted tests of alternative species delimitation using Bayes Factors are decisive. However, as argued above, such application of species delimitation is arbitrary for parthenogenetic species and, given that P. nigra is considered a pest, the naming of new species could have implications for international trade.
In general, most insect herbivores are relatively host specific, feeding on plants of only one plant family or genus [21, [101][102][103][104] leading to the question of whether some polyphagous species are, in fact, multiple host-specific lineages that are morphologically cryptic. This has been found to be the case in some scale insects (e.g. [105][106][107]) and in other insects (e.g. [108][109][110][111][112]) in which host-specific lineages have been identified in what was originally considered to be a single species. Here, however, we found no evidence that P. nigra comprises host-specific lineages because clades are not host-specific and there was sharing of hosts across lineages in allopatry and sympatry. The lineages of P. nigra might be truly polyphagous, like some other generalist insects (e.g. [113][114][115][116][117]).
In Australia, there are two lineages that appear to be ecologically distinct: one in temperate regions and the other in subtropical regions. The ecological distinctiveness indicated by the niche identity test for Australian collections is supported by the geographical distribution of other members of those lineages: the lineage restricted to temperate regions of Australia (ANZ) occurs also in New Zealand, another temperate region, and the lineage shown to be subtropical in Australia (W3) occurs in some other tropical areas (Colombia, Thailand and Malaysia). The ecological distinctiveness of the two lineages in Australia (ANZ and W3) is striking, especially given that they overlap in latitude and longitude (Fig 3). At first glance, three members of W3 seem to be outliers in the species distribution model (Fig 3B), but these individuals were collected from highly human-modified and maintained environments (e.g. supermarket carpark and vineyard).
In some localities, multiple genetic lineages of P. nigra occur in sympatry, such as W1 and W2 in Taiwan, indicating no clear ecological distinction. All lineages except ANZ appear to occur in warm regions (subtropical and tropical) and there is no clear geographic structuring.
Overall, there appears to be two ecotypes of P. nigra-temperate (ANZ) and subtropical/tropical (the rest). Only lineage ANZ appears to be ecologically differentiated from other lineages and it is therefore the only one that warrants consideration as a distinct species. The type collection of P. nigra was from Sri Lanka, a subtropical/tropical region, and thus the ANZ clade should be described and named as a new recognized taxon if taxonomic changes are to be made. With better geographic sampling, greater distinction might be discovered among the other genetic lineages that indicate ecological or phenotypic differentiation not evident here.
If lineage ANZ is to be recognized as a distinct species and the remaining lineages are to be kept under P. nigra, we will have a situation where ANZ is monophyletic but P. nigra is not (because ANZ is nested within). Paraphyly is probably inherent during early speciation in many sexual species [118] and among bacterial species [119], and a criterion of reciprocal monophyly should not be viewed as necessary for delimiting species [120][121], even though it might be favored for higher-level classifications. Unless ecological or phenotypic differences are found for the other main lineages of P. nigra (G, ICB and W1-3), we think it prudent to consider P. nigra as comprising at least two biotypes (ecotypes), one of which (ANZ) should probably be described and named as distinct species.

Consequences for quarantine
A quarantine pest is a species, biotype or strain of plant, animal or pathogen that has the potential to cause economic or other damage in an area where it does not currently occur, and which is not widespread or being officially controlled already [122]. Because of the presumed cosmopolitan distribution of P. nigra, it does not attract significant attention from quarantine authorities and is not included in the quarantine strategies of most areas (e.g. [12,22,123]). We have shown that the currently recognized P. nigra is likely a species complex, with each lineage being polyphagous: there is potential for further cross-border incursions. However, given the limited sampling available from outside Australia, it is not known whether the ANZ lineage is more widespread than just Australia and New Zealand. In particular, we have no DNA sequences for specimens from temperate regions of the northern hemisphere. This is relevant because many of these countries import hosts of P. nigra (e.g. apples, citrus and grapes) from Australia and New Zealand [124][125] where lineage ANZ is present. More generally, a greater awareness of the existence of cryptic diversity within P. nigra is warranted.  Table 1