Bridgehead Effect in the Worldwide Invasion of the Biocontrol Harlequin Ladybird

1Elucidating the routes and modalities of introductions of undesirable organisms is crucial for managers who wish to prevent new invasions.,. It is also a prerequisite to generating and testing useful hypotheses regarding the environmental and evolutionary factors responsible for biological invasions. In the large body of literature on biological invasions, a number of studies suggest that successful invasions involve a particular invasive population, which serves as the source of colonists for remote new territories,,,,. We call here this phenomenon the invasive bridgehead effect, likening the successful population of a biological invader to a military force that establishes a foothold, historically at the far side of a bridge, prior to further incursions into hostile territories. Convincing demonstrations of such an invasive bridgehead effect are however scarce because of the lack of appropriate methods to confidently reconstruct the routes of invasion and hence formally test this scenario against alternative ones. Moreover, too few invasive populations have been studied to capture the global picture of the worldwide invasion process for most species.

2The harlequin ladybird Harmonia axyridis (HA), or multicolored Asian lady beetle, is an appropriate biological model to test for the existence of an invasive bridgehead effect at a worldwide scale. Native to Asia, HA has been introduced repeatedly as a biocontrol agent against aphids since 1916 in North America,, since 1982 in Europe and since 1986 in South America. Despite recurrent intentional releases of ladybirds originating from various source populations in its native range for acclimation attempts, the species did not establish for decades. However, for unknown reasons it recently and suddenly became invasive on four different continents. Invasive populations were first recorded in eastern (Louisiana, USA) and western (Oregon, USA) North America in 1988 and 1991, respectively,. They were then recorded in Europe (Belgium) and South America (Argentina) in 2001 and in Africa (South Africa) in 2004. The species has spread widely in these areas where it has become a harmful predator of non-target arthropods, a household invader, and a pest of fruit production.

3Most of our knowledge about introduction pathways of invasive species relies on historical and observational data, which are often sparse, incomplete and sometimes misleading. Population genetics has proven to be a useful approach to reconstruct routes of introduction, highlighting how complex and counter-intuitive the real story can be,. This is especially true since the recent development of a new approach termed approximate Bayesian computation (ABC),,,,. Using molecular and historical data, this method allows performing model-based inference in a Bayesian setting for complex demographic or evolutionary scenarios such as those related to the introduction histories of invasive species, where bottleneck, multiple introductions and/or genetic admixture events are often suspected. Here, we retraced the routes of all five worldwide invasive HA populations using a rigorous quantitative analysis of microsatellite genetic variation relying on ABC methods. We compared large sets of HA invasion scenarios covering all invaded areas, taking into account historical data (e.g. dates of first observation of the outbreaks and dates of initial collection of biocontrol strains) and all potential sources (native, older outbreaks and biocontrol), accounting for the possibility of genetic admixture among them.

4Our results unambiguously point to a surprising scenario (Fig. 1). The outbreaks in eastern and western North America originated from two independent introductions from the native range (either from biocontrol or accidental introductions). The South American and African outbreaks both originated independently from eastern North America. The European outbreak also originated from eastern North America, but with substantial genetic admixture with individuals of the European biocontrol strain (estimated at 43%, 95%CI: [18%–83%]; Fig. 2a). The establishment of new invasive populations is characterized by low (South America) to minute (all other outbreaks) bottleneck severities suggesting a substantial number of founding individuals and/or a quick demographic recovery (Fig. 2b). Each choice of scenario is supported by very high posterior probabilities using two different sets of prior distributions of demographic, historic and mutation parameters (Fig. 1, Table S2 and Table S3). In all analyses the 95% CI of the most likely scenario never overlapped with those of competing scenarios. We also computed the type I and type II errors from control simulated data sets for each of the five nested analyses and found that our method selected the true scenario with high confidence and markedly low type II errors (Table 1).

5For the scenarios not complicated by substantial admixture, the raw classical statistics measuring genetic variation between populations such as pairwise F ST and assignment likelihood agree with our conclusions based on ABC methods (Table S4). On the other hand, such classical statistics are unable to detect North America and the European biocontrol strain as the sources of the admixed European invasive population. Rather, these simpler methods suggest that the latter outbreak originated from the native area only. We show with simulated data that, in case of admixture between two populations deriving from the same source, pairwise F ST and assignment likelihood values tend to select the ancestral population (here the native area) as the source of the admixed population (Fig. S2). An admixed origin of the European invasive population involving eastern North America and the European biocontrol strain as inferred from ABC is also supported by the raw allelic distributions observed in these samples (see Fig. S3 for an illustration). The ABC methods used here have hence three advantages over more standard methods based on raw measures of genetic distance: they use all the data simultaneously in inference, allow an assessment of uncertainty in all inferences and therefore provide confidence in the choice of invasion routes, and they avoid misleading biases such as those due to genetic admixture (if included in the set of compared scenarios), an increasingly acknowledged common feature of species invasions,.

6The complete panorama of the invasion history of HA strikingly points to the worldwide dissemination of a single successful invasive population (eastern North America) into remote new territories. Although ladybirds were repeatedly introduced in North America, Europe and South America for biocontrol from native Asian populations and laboratory strains, eastern North America is the proximate origin of the worldwide invasion of HA. This pattern illustrates what we call the invasive bridgehead effect, whereby a particular invasive population serves as the source of colonists for other areas. Our results have important implications for the understanding and the management of biological invasions. The history of HA reveals that prevention of new invasions should include minimizing accidental dissemination from invasive bridgehead populations (here eastern North America) rather than focusing on the native range of invasive species. For instance, we found that live HA individuals intercepted in Europe (Norway) in 2007 on imported timber most likely originated from eastern North America based on the sample being significantly differentiated genetically from all tested populations (p<1×10−6) except that of eastern North America (p = 0.060). The implication of eastern North America being an invasive bridgehead, together with the absence of remote establishment and slower local spread of the population in western North America and the long history of unsuccessful biocontrol introductions of HA from its native range, jointly suggest that an evolutionary shift triggering invasion occurred in eastern North America. It is worth pointing out that the invasive bridgehead effect is evolutionarily parsimonious: a single evolutionary shift in a single introduced population (the bridgehead) is required whereas multiple changes are required if introduced populations had become invasive independently. In Europe, the potential role of admixture with the European biocontrol strain is unknown, but the single eastern North American origin of the South American and South African outbreaks suggests that the genetic admixture observed in Europe is not required for an eastern North American propagule to establish and start an invasive population in diverse ecological contexts.

7In conclusion, our results provide a convincing demonstration of an invasive bridgehead scenario in a worldwide emblematic invasive species. That invasion can proceed via a bridgehead has important implications both for invasion theory (i.e. the number and nature of evolutionary shift(s) involved in large scale invasions) and for ongoing efforts to manage invasions by alien organisms (i.e. heightened vigilance against invasive bridgehead populations). Our study highlights the interest of new model-based methods such as approximate Bayesian computation to rigorously compare complex invasion scenarios, including those with genetic admixture, a feature of species invasions that is increasingly acknowledged to be common,,. It would be useful to revisit previously published genetic data sets on other worldwide invasive species using similar ABC approaches to examine whether bridgeheads and/or admixtures are a common mechanism driving invasion. Our ability to confidently reconstruct the routes of introduction of worldwide invaders is crucial to generating sensible hypotheses regarding the environmental and evolutionary factors responsible for biological invasions. In HA, our results suggest that an evolutionary shift triggering invasion likely occurred in the bridgehead population in eastern North America. Forthcoming studies based on quantitative genetics approaches, will shade light on the evolutionary basis of HA invasiveness in eastern North America and into the potential role of genetic admixture in Europe.

8Population samples were genotyped at 18 microsatellite markers. The samples of populations from the five invaded areas (Fig. 1 and Table S1) were collected near the sites where the outbreaks were first observed. The three samples from the native range were collected in eastern Asia, where individuals historically used for biocontrol purposes had been collected,,. The three native population samples were genetically homogeneous (mean pairwise F ST = 0.0057), and were therefore pooled into a single sample. Using native population samples individually did not change our conclusions. We genotyped numerous additional population samples collected from all invaded areas (30 additional samples) and from the native range (5 additional samples); analysis of these samples indicated that the subset of samples used in the study provides an adequate description of the main invasive and native populations (i.e. none significant to low level of genetic differentiation within the native range and within each invaded area). The European biocontrol sample (Ebc) was from the INRA rearing stock of 1987 initially collected from the native area in 1982 and used for biocontrol in Europe and South America. Genotyping of the biocontrol strains produced in European biofactories from the late 90's to the present (4 additional samples) confirmed that they are indeed derived from the original INRA strain we genotyped.

9Genetic variation within and between populations was summarized using a set of statistics traditionally employed in approximate Bayesian computation analyses (ABC),. For each population and each population pair we used the mean number of alleles per locus, the mean expected heterozygosity and the mean allelic size variance. The other statistics used were the mean ratio of the number of alleles over the range of allele sizes, pairwise F ST values, mean individual assignment likelihoods of population i assigned to population j and the maximum likelihood estimate of admixture proportion.

10We performed five serial nested ABC analyses of invasion scenarios involving successive HA outbreaks (Table S1). In analysis 1, we dealt with the introduction pathway for the first recorded outbreak in eastern North America in 1988, with native populations (either from biocontrol or accidental introduction) and/or the European biocontrol population as potential sources, thereby defining three competing scenarios (see Fig. S1). The second outbreak in western North America in 1991 was examined in analysis 2, taking into account the scenario selected in analysis 1. For this analysis, there were three possible sources: the first outbreak (eastern North America), the native range, and the European biological control strain (Table S1). With the potential for admixture between them, this gave 6 competing scenarios. The European and South American outbreaks in 2001 were addressed in analyses 3 and 4, respectively (10 scenarios for each outbreak), taking into account the scenario selected in analysis 1 and 2. The year of first observation in Europe and South America was the same, thus we assumed that one could not be a potential source of the other. Nevertheless, we performed a specific ABC analysis to test this assumption which confirmed the independence of these two outbreaks (P = 0.773 with CI = 0.555–0.991). Finally, the African outbreak in 2004 was considered in analysis 5 (21 scenarios), taking into account the scenarios selected in analyses 1, 2, 3 and 4.

11The ABC analyses were performed using parameter values drawn from the prior distributions described in Table S2 (prior set 1), and by simulating 106 microsatellite data sets for each competing scenario in the first four analyses and 5×105 data sets in analysis 5 because of the high number of scenarios (21) and summary statistics (170) which made a larger analysis computationally unfeasible. For each of the five analyses, we estimated the posterior probabilities of the competing scenarios using a polychotomous logistic regression on the 0.1% of simulated data sets closest to the observed data set. The selected scenario is that with the highest significant probability value with a non-overlapping 95% confidence interval. We estimated the posterior distributions of demographic parameters under the final HA invasion scenario presented in Fig. 1 using a local linear regression on the 1% closest of 5×106 simulated data sets,.

12We evaluated the robustness of our inferences by running a second analysis with an alternative set of priors (prior set 2 in Table S2) and by estimating the posterior probabilities of scenarios using the 0.1% or 1% closest simulated data sets for both sets of priors (Table S2). For each of our five ABC analyses, we evaluated the ability of the methodology to correctly select the true scenario by analyzing test data sets simulated from the various competing scenarios with the same number of loci and individuals as in the real data set. For each scenario, one hundred such data sets were simulated using parameter values drawn from the same probability distributions as the priors (prior set 1, Table S2). Posterior probabilities of each competing scenario were estimated for each simulated test data set using the 0.1% closest data sets. These probabilities were used to compute type I and II errors in the selection of scenarios.

13Most ABC computations were processed using a modified version of the Windows package DIYABC. Some of the computations for assessing confidence in scenario selection (i.e. those for analysis 5) necessitated the development of specific Linux programs and scripts running on a computer cluster (available under request from A.E.).

14We ran computer simulations to tackle the question of a high frequency of erroneous selection of the source population in a situation of genetic admixture when considering only the raw values of genetic differentiation statistics (i.e. pairwise F st and assignment likelihood values). We used the program DIYABC to simulate genetic data sets under a scenario of invasion with admixture similar to the one considered in the case of the invasion of Europe by HA. We simulated 6×104 data sets drawing the parameter values in the prior set 1 (Table S2), except the admixture rate (ar) which was drawn in a discrete {0, 0.1, …,1.0} distribution instead of a uniform [0.1;0.9] distribution. The number of sampled diploid individuals was fixed to 30 in all populations. The deduced origin of the admixed introduced population is the population with which the F ST-value is the smallest or the assignment likelihood value is the highest. We then calculated the proportion of simulated data sets, as a function of the admixture rate, for which the deduced origin was the ancestral population of the two actual source populations (i.e. the native area in the case of the invasion of Europe by HA), considering either F ST or assignment likelihood values (see Fig. S2).

15Finally, we also genotyped 200 live HA individuals intercepted in Europe (Åndalsnes, Norway) in 2007 on timber imported from North America. We tested for the genetic differentiation, using Fisher's exact tests, between this interception HA sample and each population sample used in the ABC analysis to elucidate its most likely origin.