Multiple introductions and environmental factors affecting the establishment of invasive species on a volcanic island

Invasive species pose signi ﬁ cant challenges to local biodiversity and ecosystem function, especially on islands. Understanding the factors affecting the establishment of invasive species and how these relate to their genetic background is crucial to improve our ability to manage biological invasions. Here, we performed a phylogeographic study of two cosmopolitan megascolecid earthworms of Asian origin: Amynthas gracilis and Amynthas corticis at 38 localities on S ~ ao Miguel Island in the Azores archipelago (Portugal). Samples from putative source populations in China, Taiwan, Malaysia, as well as ‘ outlier ’ populations in USA, Mexico, Brazil and Spain were also included, resulting in a total of 565 earthworms genotyped at the mitochondrial cytochrome oxidase I (COI) and 16S ribosomal RNA genes. Soils were characterised for elemental composition, water holding capacity, organic matter content, texture and pH, and some habitat features were recorded. Both species showed a wide distribution across S ~ ao Miguel and their abundances were negatively associated, suggesting spatial segregation/competition, with the parthenogenetic A. corticis being relatively more successful. The presence of multiple mitochondrial lineages within each species, one of them found exclusively in the Azores, suggests a complex invasion history. Environmental factors affected the establishment of the different lineages, with metal concentrations, topographical elevation and the degree of human in ﬂ uence being differently linked to their abundances. Lineage diversity was negatively correlated with metal concentrations. These results emphasise the importance of genetically characterising invasive species to better understand their invasion patterns. © 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Invasive species can adversely affect local biodiversity due to alterations of recipient ecosystems, impacts on native species, such as competition, predation or hybridisation, or as carriers of disease (Pejchar and Mooney, 2009).Global environmental change and globalisation of trade networks mean that introductions are more likely to occur and their proper prediction, prevention and management is crucial.Therefore understanding the influence of the environment, geography and genetic features on colonisation and invasiveness of species should be a priority.Islands are particularly vulnerable to invasions, which are the main cause of population declines and species extinctions as well as having substantial socioeconomic impacts within these areas (Reaser et al., 2007).Invasive species have shown rapid adaptations to new biotic and abiotic environments and much of this evidence has been observed on islands, highlighting them as evolutionary hotspots (Mooney and Cleland, 2001).
Terrestrial invertebrates have been identified as one of the groups with the most profound ecological and economic impact through their roles in environmental processes that give rise to ecosystem services (Vila et al., 2010).However, attention has been focused on invasions above ground while soil organisms have been largely unexplored, even though their invasions may profoundly affect soil ecosystems (Hendrix et al., 2008).Alterations of the soil structure by earthworms can cause a cascade of ecological effects (Frelich et al., 2006).In the last decade several authors have studied the effects of introduced earthworms on temperate and tropical regions, showing their capacity to greatly modify nutrient and organic matter dynamics, above-and below-ground community structures, and soil structural properties; the role of humans in earthworm dispersal has also been shown to be significant (Hendrix et al., 2008 and references therein).A well-documented case concerns the invasion of European earthworms (Lumbricidae) in North America, which has resulted in a change in forest floor litter dynamics with resulting effects on ecosystem processes (e.g., Bohlen et al., 2004;Hale et al., 2005;Frelich et al., 2006).Less understood are the effects of megascolecid species, although Burtelow et al. (1998) showed that Amynthas gracilis in particular has impacted C and N fluxes in soils of the northeastern United States.
The Azores archipelago is located in the North Atlantic Ocean and comprises nine islands.Of these, São Miguel is the largest and is made up of active volcanic areas including fumarolic fields and cold and thermal springs with soil-diffuse degassing (Viveiros et al., 2009).Humans first populated the island in the late 1430's (Santos et al., 2005) and dramatic changes in land-use started after colonization.Native vegetation at low and middle altitudes became extinct or was modified and exotic plants and animals were introduced (Borges et al., 2006 and references therein).Among the Macaronesian archipelagos, the Azores contain fewer single-island endemic species (Amorim et al., 2012), which has been referred to as the 'Azores diversity enigma' (Carine and Schaefer, 2010).This relatively depauperate native species diversity could make the Azores archipelago more susceptible to successful invasions because more niches are unoccupied and, for example, 70% of the current flora and 58% of the arthropod fauna are exotic (Borges et al., 2006).
Most phylogeographic studies of earthworms to date have been carried out within the natural ranges of the species, including megascolecids (e.g., Chang and Chen, 2005;Chang et al., 2008;Buckley et al., 2011;Fernandez et al., 2011;Novo et al., 2011).Cameron et al. (2008) and Cunha et al. (2014) investigated the population genetic structure of invasive earthworms in the US and the Azores respectively, in relation to human-mediated dispersal and landscape features.Porco et al. (2013) and Shekhovtsov et al. (2014) have also applied DNA barcoding to study invasions of earthworms in the US and Kamchatka.More studies are needed to better understand the dynamic forces affecting earthworm invasions given the potential consequences for native species and ecosystems.
The Megascolecidae is the largest earthworm family and has a suggested Pangean origin (Jamieson et al., 2002).The genus Amynthas is thought to be native to the eastern Palearctic, with species being described as cosmopolitan, peregrine and invasive.Amynthas corticis and A. gracilis are amongst the most invasive earthworms on earth, mainly due to their inherent environmental plasticity, rapid growth, high reproductive performance and relatively large adult body size (Burtelow et al., 1998;Garcia and Fragoso, 2002).Additionally, some lineages of A. corticis including those from the Azores (unpublished observations by the authors) are parthenogenetic.This mode of reproduction can facilitate invasiveness (Terhivuo and Saura, 2006).For example, parthenogenetic species can found stable populations initiated by a single individual in novel environments.The benefit of uniparental reproduction is highest after long-distance dispersal (Baker, 1965(Baker, , 1967)), which would be the case following human introductions.
It is not presently known whether the introduction of megascolecid worms to the Azores was a single or multiple event(s); neither is it known whether the invaders originated from one or more populations in the Eastern Palearctic (Hendrix et al., 2008).Given that the establishment on São Miguel Island by both A. corticis and A. gracilis is probably a relatively recent humanmediated event, they constitute an ideal model to test how the genetic background of an invader responds to local geography and the chemical characteristics of a new environment, in this case an island characterised by a volcanic-driven topography and geology.Physical features are known to affect earthworm distributions even on a very small scale (Rossi et al., 1997;Nuutinen et al., 1998;Hernandez et al., 2007;Jim enez et al., 2014).
Our aim was to study the phylogeography of megascolecids in São Miguel Island in the Azores (Portugal), placing their human introduction into a global context and exploring the likely factors affecting their establishment on this island.Specifically we aimed to: i) compare the haplotypes of Azorean A. gracilis and A. corticis to those of both species from other global sources (including Asian) in an attempt to detect different introduction events; ii) study the extent of genetic variability at the landscape scale in São Miguel by identifying genetic lineages within both species and their spatial distributions; iii) explore the relationships between the abundances of different lineages and environmental characteristics.Addressing these issues has allowed us to ascertain whether the invasion success of Amynthas on São Miguel is attributed to just one highly adapted genetic lineage per species, or to multiple lineages with potentially different adaptative capacities to spatially heterogeneous environmental conditions.

Sampling
Samples were collected from throughout São Miguel Island (Azores, Portugal).Its area is 757 km 2 and includes three active volcanoes: Furnas (eastern end, ca. 1 million years old), Sete Cidades (western end, 550e750 thousand years old) and Fogo (center, ca.350 thousand years old) (Gomes et al., 2006;Cunha et al., 2014).The oldest area of the island, to the east is ca. 4 million years old (Harris et al., 2013).The climate in São Miguel is oceanic and temperate being strongly influenced by altitude (Ricardo et al., 1977).The summits are higher in the eastern part (Pico da Vara, 1080 m) than in the western part (Eguas, 873 m), and the central part is the lowest (maximum of 400 m) (Louvat and All egre, 1998).Mean annual temperatures range between 11.5 C at the summit of Agua de Pau volcano (800e900 m) and 17.3 Cin Ponta Delgada (sea level).Annual rainfall and average relative humidity increase with elevation: rainfall ranges from 1000e1500 mm/yr in the driest part of the island (the south coast) to 3000 mm/ yr at Lagoa do Fogo; relative humidity ranges from 77e78% on the coast to 87e88% at the highest elevation (Malucelli et al., 2006).Around 500 megascolecid earthworms were collected from 37 localities (Supplementary Table 1, Fig. 1B, C).The sampling effort was standardized in all the locations in order to record relative (but not absolute) abundances and to compare these between sites.These relative numbers are named abundances hereafter within the text.We sampled at least 8 quadrats (50 cm Â 50 cm) randomly distributed in the sampling area that were separated by at least 3 m.Sampling effort was limited by time and people making the combination uniform among sites (generally 4 people digging for 60 min) with the search targeted to Amynthas specimens.All individuals were hand collected, washed in distilled water, weighted in vivo and subsequently fixed and preserved in ca.96% EtOH.A portion of the integument (ca. 25 mg) was cleansed and preserved at ¡20 C for DNA extraction.Samples from other countries in the species' global distribution, including from the presumed native range (Malaysia: MAL, China: YN, HB; Taiwan: TW), were donated by colleagues (see Supplementary

Genetic data analyses
Sequences were aligned in CLUSTALX v. 2.0.12 (Thompson et al., 1997) using default settings.Haplotypes were retrieved in DNAcollapser (FaBox 1.41) and a Bayesian phylogeny was estimated for each gene with MrBayes v. 3.2.2(Ronquist and Huelsenbeck, 2003).GTR þ I þ G, the best-fit substitution model implemented in MrBayes, as indicated by jModelTest2 (Darriba et al., 2012), was used.Parameters in MrBayes were set to forty million generations and a frequency of sampling of 4,000, retrieving 10,000 trees.The analysis was run twice and 25% of the trees were discarded as burnin.The remaining trees were combined to find the maximum a posteriori probability estimate for each clade.Sequences from the earthworm Aporrectodea trapezoides (JF313607; HQ621864), which belongs to a different family, Lumbricidae, were used as an outgroup.We performed the same analysis for the COI and 16S datasets separately and concatenated.A COI haplotype network was constructed using TCS version 1.21 with statistical parsimony and a 50-step connection limit (Clement et al., 2000).
An ultrametric tree was generated for COI using BEAST v. 1.7.5 (Drummond et al., 2012) with the evolutionary rate calculated by Chang et al. (2008) for Megascolecidae (2.4 substitutions/MY).The analysis was conducted under a constant size coalescent model, an uncorrelated lognormal relaxed clock and by setting GTR þ I þ G as the evolutionary model.We performed 30 parallel runs, each of which included 10 million generations, sampling each 1,000th generation.We combined tree and log files in Logcombiner v.1.7 by resampling at a lower frequency (30,000) and visualized the results using Tracer v. 1.5 (Rambaut and Drummond, 2007).The final tree was generated by TreeAnnotator v.1.7 with a burnin of 2000.
Island lineages have been observed to increase the rate of nonsynonymous substitutions in mitochondrial protein coding genes (Johnson and Seger, 2001;Woolfit and Bromham, 2005).The percentage of non-synonymous substitutions was calculated after translating COI sequences in DNAsp v 5.10.1 (Librado and Rozas, 2009) using the invertebrate mitochondrial code.

Environmental characterisation
Sampling sites were characterised by elemental analysis together with environmental variables, such as moisture, organic matter, texture and pH.Codes from 0 to 4 (five categories) were assigned for covering percentage of grass, moss, ferns and stones ( 0 ¼ absence; 1 ¼ 1e25%; 2 ¼ 26e50%; 3 ¼ 51e75%; 4 ¼ 76e100%).Codes were assigned as well to the degree of human interference at the sites (0 ¼ absence, natural areas; 1 ¼ scarce, grasslands and groves; 2 ¼ moderate, edges of unpaved roads or forests paths; 3 ¼ high, foraged with grazing livestock, areas near houses; 4 ¼ very high, active crops and gardens).
Soils were air-dried prior to return to the UK for analysis.In the UK, samples were again air-dried and sieved to <2 mm.Subsamples were oven dried at 105 C and moisture content determined.Concentrations are expressed on an oven dried soil basis.For pH measurement 10 g of air-dried soil was added to 25 ml deionised water and shaken for 15 min; pH was then measured (ISO, 2005).Loss on ignition was used as a proxy for organic matter content.Oven dried soil was weighed into a crucible and ignited overnight at 500 C in a muffle furnace after which mass loss was determined (Rowell, 1994).Bulk metal content was determined by aqua regia digestion (Arnold et al., 2008) followed by analysis using a Perkin Elmer Optima 7300 DV inductively coupled plasma optical emission spectrometry instrument.Detection limits (Supplementary Table 1) were calculated as the mean plus 6 times the standard deviation on ten replicate analyses of the blank calibration standard (Gill, 1997).Method blanks were run and results were blank corrected where appropriate (for Ca, Mg, Na and Zn).For quality control we digested NIST Certified reference material Soil SRM 2709 and compared our results to those obtained after a US EPA 3050 acid digest.Recoveries were between 90 and 110% for Ca, Cr, Cu, Fe, Mg, Mn, Na, Ni, Pb, Sr and Zn.
Soil particle size distributions were determined using a Malvern MasterSizer 2000 with a Hydro2000MU wet dispersion unit.One to two grammes of air-dried, <2 mm sieved-soil were analysed with an obscuration of between 5 and 25%.Accurate instrument performance was confirmed using a Malvern 15e150 mm quality audit standard and "general purpose sand, 40e100 mesh" purchased from Fisher.Rather than measuring soil moisture content in the field during sampling, water holding capacity of the soil samples was determined following the method in ISO guideline 11274 (ISO, 1992) to give an indication of the relative potential moisture content of the different samples.

Statistical analyses
All statistical analyses were performed in R version 3.1.1(R Development Core Team, 2014), using stats, car, MASS, psych (Revelle, 2011) and vegan (Oksanen et al., 2007) libraries.We performed a preliminary analysis using simple correlations between the initial set of 31 variables (see Supplementary Table 1) and species and lineage abundances.We removed those variables not showing significant relationships and performed subsequent analyses with a set of 19.Principal component analysis (PCA) was then applied to describe the main sources of variation and relationships between the retained habitat features.The 'varimax' rotation method was used to increase the interpretation of axes and the number of PCA axes examined was determined by sequential eigenvalue reduction (Legendre and Legendre, 1998).The relationship between the abundance of each earthworm species or lineage as well as lineage diversity in a location was examined using generalised linear models (GLM) with altitude, the abundance of the other earthworm species, and the two stressor gradients from PCA analysis (i.e.PC1 and PC2) as fixed factors.As previously validated in rivers (Maceda-Veiga et al., 2013) and highlighted by K€ orner (2007), altitude was used as a surrogate for the position of the sampling sites on the mountain, and summarised the role of natural spatial gradients in earthworm abundance.Logtransformation was applied to altitude and earthworm abundances to reduce heterogeneity of variances and increase model fitting.A quasi-poisson error distribution was used for abundances, whereas Gaussian errors were assumed for lineage diversity.Significance of explanatory variables (factor effects) in GLMs was assessed using a likelihood ratio test.Best models were selected using a manual stepwise backward deletion of non-significant terms from the full global models containing all four factors and interactions.Models were validated with qeq plots of residuals and plotting fitted vs. predicted values, and spatial correlation was examined using Durbin Watson tests.As the modelling approach employed lacks a true variation coefficient (i.e.R 2 ), we calculated a pseudo-R 2 coefficient as follows: (null deviance e residual deviance)/null deviance.To complement the results of GLMs and to test the robustness of our results, we performed a hierarchical partitioning analysis ("hier.part"function in R) (Walsh and Mac Nally, 2011) using the error distributions validated in the GLM approach.An advantage of this approach is that it controls for colinearity amongst explanatory variables that, even at low levels, can cause variance inflation and lead to erroneous conclusions (Graham, 2003).Other modelling criteria such as AIC and model averaging are discouraged because co-linearity results in biased parameter estimates (Freckleton, 2011).We assessed the significance of HP models using a randomization test for hierarchical partitioning analysis (function "rand.hp" in R).Significance in HP analysis was based on the upper 0.95 confidence interval, but it was reached at p < 0.05 in the remaining statistical procedures.The mapping software SURFER © (Golden Software, Golden, Colorado) was used to convert the total metal load (mol$kg À1 ), altitude and earthworm abundance data sets into a series of 2D contour maps.

Species and haplotype geographic distributions
Most of the sampled megascolecids that showed a wide distribution within São Miguel Island (Fig. 1B), belonged to A. gracilis (n ¼ 110) or A. corticis (n ¼ 315), whereas a small group of unclassified megascolecids (n ¼ 22) occurred in four localities, and were dominant in only one.In total, 565 COI sequences were generated, including those samples received from the various global sources.Some localities were dominated by a single species or even a single haplotype in the case of A. corticis.The distribution of COI haplotypes for A. corticis is shown in Fig. 1C.Assignment of COI haplotypes is shown in Supplementary Tables 3 and 4 (Azorean and non-Azorean samples, respectively), distribution of species and lineages per site in Supplementary Table 1 and assignment of 16S haplotypes in Supplementary Table 5.

Haplotype network
The networks generated by TCS with the COI and 16S genes are shown in Fig. 2.F orA.gracilis a common COI haplotype was found (HAP-1, n ¼ 152), which was shared with non-Azorean specimens (Brazil, Malaysia, Mexico, Taiwan, USA and China).Some individuals from China, Taiwan and two populations in São Miguel showed COI haplotypes separated by a single nucleotide polymorphism (SNP) when compared to HAP-1, but they all shared a 16S haplotype implying a common origin.This lineage was notated a, and was separated by 26 SNPs (0.043 uncorrected p-distance) from lineage b, composed of two COI haplotypes (HAP-8 and HAP-9) but sharing a 16S sequence.The last group comprised three A. gracilis populations from São Miguel (S25, S26, S27) and two non-Azorean specimens (USA and Brazil).The latter did not amplify for COI but showed the same 16S sequence.
For A. corticis three main haplogroups were found.Lineage A comprised 24 individuals from São Miguel with COI haplotype HAP-29 and specimens from Spain and USA with different COI haplotypes but sharing a 16S sequence.Lineage B included 155 individuals from the Azores, Brazil, USA, and China.Interestingly specimens from Brazil and USA showed a different 16S sequence (Fig. 2).Also a single individual at site S22 possessed different COI and 16S haplotypes.Lineages A and B were the most distinct (43 SNPs for COI, 0.069 uncorrected p-distance) but an intermediate haplotype from a Chinese sample (HAP-37) was found.Finally, lineage C possessed only 5 SNPs different to lineage B (i.e.0.008 COI uncorrected p-distance) and comprised only individuals collected from the Azores (n ¼ 164).Lineage B seems to be the origin of the other groups according to TCS, although intermediates may be missing from our sampling.1 and 2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Phylogenetic trees and dating
The topologies of the phylogenetic trees were consistent for all analyses (separate or concatenated datasets; Fig. 3) and provided concordant groupings with the network analyses.Dates recovered using BEAST imply that the divergences between the identified lineages are very old and certainly pre-dated human introductory events.For A. gracilis, the a and b split was dated to around 1.14 Mya.
For A. corticis, lineage A was dated as having diverged at 1.93 Mya from B and C, with the latter pair diverging approximately 267 Kya.  1 and 2 and haplotype codes in Supplementary Tables 4, 5 and 7.

Non-synonymous substitutions
Relatively low numbers of non-synonymous substitutions were generally found between pairs of sequences ranging from 0 to 16.7%, even when comparing the different megascolecid species analysed using up to 126 substitutions.However, percentages of non-synonymous substitutions between 28.6 and 100% were found between A. corticis haplotypes involving the singly observed haplotypes 28, 31, 32 and 37. Nevertheless it is worth highlighting that between the A. corticis haplotypes 30 and 33, which seem to be well-established because of their geographical distribution and abundances, we found 40% non-synonymous substitutions (two out of five).

Relationship between earthworm abundance and environmental variables
Elemental concentrations and habitat features were summarised using PCA, producing two significant axes explaining 72% of the total variation (Table 1).PC1 was found to be loaded on element concentrations and accounted for 52% of the variation, whereas PC2 included habitat features, representing the degree of human influence, and explained 20% of the variation.As the environmental gradients built with the PCA and the abundance of earthworms were significantly correlated with altitude (Spearman's rho range ¼ 0.3e0.7), it was necessary to disentangle this relationship and analyse the independent and joint effects of these predictors on the abundance of each species of earthworm or clade groups.The results of GLM and HP models were mostly concordant indicating that co-linearity between predictors played a minor role in our data-set.Discrepancies were only observed between GLM and HP when the former retained as significant the PC2 gradient due to its highly significant and positive relationship with altitude (r ¼ 0.7, p < 0.01).None of the GLMs exhibited spatial autocorrelation (i.e.dependency amongst observations in a geographic space) in the residuals according to Durbin Watson tests (D-W>1.7,All p > 0.10).
At least one of the environmental gradients from the PCA and/or the abundance of the other earthworm species was retained as significant, and made the largest contribution to the abundance of the earthworm species included as a response variable (Tables 2  and 3).Interestingly, the abundance of A. gracilis was negatively associated with the abundance of A. corticis (À1.43 ± 0.43; c 2 ¼ 10.72, df ¼ 1, P < 0.001), suggesting spatial segregation/ competition between the two invasive species (Fig. 4).The abundance of A. gracilis also made the largest contribution to explain the abundance of A. corticis (À0.69 ± 0.21; c 2 ¼ 13.06, df ¼ 1, P < 0.001) followed by the metal concentrations (PC1) (À0.20 ± 0.02; c 2 ¼ 5.84, df ¼ 1, P ¼ 0.02).When the abundance of each lineage was considered in the analysis, metal concentration was also identified as an important factor explaining the abundance of A. corticis from lineage C and A. gracilis from lineage b (Table 2).
A. corticis lineage C worms were, however, mostly associated with higher altitude (1.96 ± 0.56; c 2 ¼ 15.22, df ¼ 1, P < 0.001) and according to HP, with localities with a high degree of habitat naturalness (PC2) (Table 3, Fig. 5).The opposite relationship was found between habitat naturalness and the abundance of A. corticis from lineage A (1.78 ± 0.70; c 2 ¼ 10.73, df ¼ 1, P ¼ 0.001) and A. gracilis from lineage b (1.56 ± 0.65; c 2 ¼ 8.12, df ¼ 1, P ¼ 0.004), both being positively related with human influence.This illustrates that habitat requirements may differ between lineages from the same invasive species.Finally, we found that lineage diversity was affected negatively by the metal load (PC1) (Table 2).

Discussion
Our results show that Amynthas species are successful invaders of São Miguel Island.They appear to be present in all habitable soils across the island.Lineage analysis based on mitochondrial genes, and compared with the genotypes of non-Azorean populations, strongly indicates that the extant populations result from multiple introductions.The distributions of the two main megascolecid species A. gracilis and A. corticis are negatively correlated, indicating spatial segregation of the invaders.A. corticis seems to be more widespread and therefore more successful in occupying the diverse habitats found on the island; this was particularly the case with A. corticis lineage C, which was the only lineage which we were unable to identify in any non-Azorean population.We found that besides the negative correlation among the two invasive species, altitude, soil parameters and habitat were also linked to the distribution patterns of the different mitochondrial lineages and to lineage diversity.

Amynthas species as successful invaders of São Miguel Island
Amynthas species are widely distributed across São Miguel as shown by our widespread sampling.Earthworms of this genus have previously shown high invasive potential, and are now widely spread over tropical and temperate geographic regions around the globe although with particular habitat preferences mainly associated with anthropogenic disturbance (Hendrix et al., 2008).Different characteristics of certain earthworm species have been noted to contribute to successful invasion, such as the ability to aestivate, physiological plasticity, and parthenogenesis (Lavelle et al., 1999;Edwards, 2004).Amynthas could be particularly aided by its high mobility and reproduction rate (Burtelow et al., 1998;Garcia and Fragoso, 2002), and parthenogenesis in the case of A. corticis (see Section 4.3).They are also epigeic species (surface dwellers), which disperse actively through the soil but are also prone to physical dispersal forces such as flooding, wind or phoresy (Eijsackers, 2011).Success of Amynthas species as invaders in the particular setting of São Miguel may be due to several factors identified here.Firstly, multiple introduction events could promote adaptation to the new environment by enhancing genetic diversity (Cameron et al., 2008).Such events would counteract the genetic bottleneck resulting from a singular introduction (Caron et al., 2013).We have inferred at least four distinct Amynthas introduction events in São Miguel, two for each species (i.e. A. corticis A, B; A. gracilis a, b).Lineage C in A. corticis remained of unknown origin (see below).More COI haplotypes were found but 16S clustered them into those five groups.The number of times that a species arrives has been listed as an important factor for a successful colonizer, together with endurance of environmental conditions and the range of habitats in which a species can establish (Eijsackers, 2011).Habitat matching facilitates invasion (Hendrix et al., 2006).As stated by Beddard (1912), temperate species tend to invade temperate regions and tropical species tend to invade tropical regions.Consequently, the similarity of the Azores subtropical climate to that of the Amynthas source-regions in Asia has likely benefitted the worms, enabling the pioneers to survive and subsequently establish sustainable populations.Also, being epigeic species they are expected to have found the shallow but relatively fertile Azorean volcanic soils, with ample and widespread associated surface litter (pers.obs.), compatible with their ecological requirements.

Negative correlation of abundances of both invaders
The negative relationship between the abundances of A. gracilis and A. corticis was found in this study to make the largest independent contribution to explain their respective abundances in São Miguel.Many studies have analysed the effects of invasive species on native relatives, including niche displacement or competitive exclusion (e.g., Mooney and Cleland, 2001;Johnson et al., 2009), although this effect is not as clear in earthworms (Hendrix et al., 2006;Gonz alez et al., 2006).However, the interaction between two invasive species from the same taxon has not been studied extensively, although some studies of multiple invaders from different animal groups provide evidence of facilitation (meltdown hypothesis) among invasive species (Ricciardi, 2001;O'Dowd et al., 2003;Grosholz, 2005), and only a few show negative interactions (Ross et al., 2004;Johnson et al., 2009).Other studies have shown spatial segregation of introduced related species, indicating a similar use of resources by both (Allen and Shea, 2006;Rauschert, 2006;Monroy et al., 2014).Hitherto, the laboratory-based study on lumbricids by Butt (1998) appears to be the only one that has examined interactions between potential invasive earthworm species; he concluded that the interaction showed the hallmarks of  1 in the same units.facilitation.To our knowledge the present field-based study of two closely related megascolecids is the first record of a negative correlation in the abundances of the invaders.We cannot discount the possibility of an initial facilitation between the species, but their similar resource requirements along with different adaptation capabilities (see Section 4.4) may have caused their subsequent exclusion.Uvarov (2009) concluded that, unlike deep-burrowing anecic earthworms, which mainly show positive interactions, epigeic earthworms are mostly competitive, and especially so in the acquisition of food (Abbott, 1980).However, the extent of competition among Amynthas species in São Miguel is uncertain, and possibly confounded by edaphic preferences, because they have been found together in some sites.Richard et al. (2012) postulated that interspecific interactions are important factors driving the spatial distribution of earthworms at a local scale.Future studies on the distribution of the earthworm fauna of São Miguel Island could possibly benefit from focusing a more highly resolved spatial scale in order to shed light on the inter-play between interspecific and abiotic environmental parameters.

A. corticis as a more successful colonizer than A. gracilis
A. corticis showed a broader distribution across São Miguel Island and has colonised a wider variety of habitats than A. gracilis.Our results suggest that the different lineages of A. corticis have colonised different soil types (see Section 4.4), which may indicate higher physiological plasticity and adaptation potential.A. corticis is recognized to be widely tolerant to environmental factors (Fragoso et al., 1999).Between the two species, A. corticis has greater mobility (pers.observ.)and also reproduces parthenogenetically in the Azores, unlike the sexually reproducing A. gracilis (author's work in progress).Parthenogenetic earthworms, most of them also known to be polyploid (Muldal, 1952;Jaenike and Selander, 1979;Viktorov, 1997;Lavelle et al., 1999;Terhivuo and Saura, 2006), have several fitness advantages that would promote colonization and an ability to sustain stable populations in new environments since a single organism can establish an entire population.Higher levels of genetic variability and exceptionally plastic genomes that permit rapid adaptation, supported by higher reproductive rates enables the production of a disproportionate number of offspring that, in turn, repeat the parental dispersion behaviour and reproductive performance, thus extending and sustaining niche occupation (Hughes, 1989;Terhivuo and Saura, 2006;Díaz Cosín. et al., 2011).As such, these differences in biology between the two species may contribute to the relative success of A. corticis in colonising new habitats.

Environmental requirements of different lineages
Individuals belonging to the five main mitochondrial lineages identified for A. gracilis and A. corticis have different habitat requirements that would have been overlooked by the analysis at species level.These results stress the importance of the genetic characterization of invasive species to better understand their invasion patterns.Other studies have also found that genetic lineages of the same or related species have different ecological requirements (VanDyke et al., 2004;Mozdzer and Zieman, 2010;Gul, 2013).Such genetically-linked differences probably modulate their invasive potential.In the current study, metal concentration in soil, as broadly defined by PC1, was negatively associated with the abundance of lineage C in A. corticis and b in A. gracilis.Soil features, such as texture and elemental composition, seem to be more important for colonization by earthworms than inherent ecological characteristics (such as r and K strategies) and they also affect dispersal (Eijsackers, 2011).Invasive earthworms can avoid adverse conditions, such as extremely acidic, extremely alkaline or contaminated soils (Eijsackers, 2011).In the current study altitude was a significant determinant of the abundance of lineage C of A. corticis.A positive effect of altitude on this lineage, plus its negative relationship with human influence and soil metal concentrations, seems to indicate that this lineage is related to the most natural conditions in the island.Suarez et al. (2006) showed that exotic earthworms were more likely to occur on low and flat sites in comparison to more remote, steep sites; lineage C of A. corticis may be an exception that proves the rule.However, according to the geographic parthenogenesis theory (Vandel, 1928), parthenogens can colonize higher altitudes than their nonparthenogenetic conspecifics.They also have broader distribution areas (Cosendai et al., 2013 and references therein), which could explain why at least one lineage of A. corticis inhabits higher altitudes than other A. corticis lineages and each of the A. gracilis lineages.
Species invasions are frequently underpinned by anthropogenic drivers and hence the establishment of invasive species may frequently be related to measures of human influence.(Hendrix et al., 2008).Abundance of lineage b of A. gracilis and lineage A of A. corticis showed a positive relationship with the degree of human influence supporting this statement only for those lineages.Habitats altered by humans are also characterized by increased species richness and dominance by alien species (Borges et al., 2006).We found no relationship between the number of lineages present at a site and the degree of human influence but our results show that more lineages were present in places with lower metal concentrations.
The fact that different lineages can tolerate different conditions could positively influence the success of the invasion of Amynthas in the Azores.Thus representatives of the two species and their different lineages appear to inhabit the full range of available volcanic soil types on São Miguel Island.

Lineage C of A. corticis has no known origin
All the haplotypes, except for lineage C (haplotype 33 for COI) in A. corticis, could be detected outside the Azores and therefore must have been introduced independently to the island, having diverged beforehand.Another exception was haplotype 31, but it was found in only one individual.Lineage C is the most widespread in São Miguel and is related to the most natural conditions in the island (see above).It could be argued that lineage C has evolved within the island and is thus related to the island's natural environmental characteristics.In fact, the percentage of non-synonymous substitutions of COI between lineages B and C (haplotype 30 and haplotype 33) of A. corticis in comparison with the remainder pairwise values, supports an island origin for the latter lineage.It has been shown that island lineages have increased rates of nonsynonymous substitutions in mitochondrial protein coding genes, apparently starting soon after the isolation (Johnson and Seger, 2001).This phenomenon contradicts the dates given by the calibrated tree, which indicates an ancient separation of these two lineages and, therefore, a high probability that lineage C evolved outside the Azores and was subsequently introduced.It is conceivable that the tree dates could be valid for natural colonization, but such an event seems improbable due to the oceanic isolation of the entire Azorean archipelago (but see Sharma and Giribet, 2012).However evolutionary rates cannot always be extrapolated from one group to another since different rates occur even within the same family (Novo et al., 2012).Also, rapid evolution has been proposed in island environments (Aleixandre et al., 2013), even over relatively short time scales from a few decades up to several thousands of years (Millien, 2006).This could justify the old dates in the tree for the split between lineages B and C when applying non-specific rates.Other studies show the lack of evidence of an increase in evolutionary rates in island taxa (Bromham and Woolfit, 2004), but still demonstrate the existence of higher ratios of non-synonymous to synonymous substitutions (Woolfit and Bromham, 2005).It is prudent to conclude that neither the evolution of lineage C within the island, nor evolution outside the island and coupled with posterior human-mediated introduction, can be discarded after this study.Wider collection in the home range of the species may be useful to try and discover the existence of lineage C outside the island.

Conclusion
The studied Amynthas species are successful invaders in the Azores due to their inherent biological and ecological characteristics and enhanced by the fact that their presence is the result of multiple introductions.The relatively higher success of A. corticis may be due to parthenogenesis and higher plasticity.One of its lineages has an unknown origin and seems to inhabit more natural, undisturbed, areas on the island.Abundances of A. gracilis and A. corticis are negatively correlated, suggesting segregation or competition between both invaders as opposed to the 'meltdown' hypothesis.Distinct genetic lineages showed different environmental preferences, with metal concentrations, degree of human influence and altitude linked to their abundances.Lineage diversity was higher in places with lower metal concentrations.Our results highlight the need for genetic studies in order to better understand invasion patterns since a single species can present complex responses to new environments.The inclusion of nuclear markers such as microsatellites or RAD-tag analysis of Amynthas species found in the Azores populations could further develop understanding about the more recent adaptation mechanisms to further explain observed distribution patterns and possible competition of Amynthas species in São Miguel.Finally, we envisage that a suite of genotyping and bioinformatics approaches such as those used in the present study, combined with the mentioned nuclear markers, could be invaluable tools to help characterise and manage invasive species in diverse aquatic as well as terrestrial ecosystems worldwide.

Fig. 1 .
Fig. 1.Geographic distribution of the samples analysed.A: Sample sites outside Azores: São Miguel location is represented by a square.B: Sample sites within São Miguel Island (Azores).Proportion of different species found is shown.C: Distribution of haplotypes of the COI gene found within Amynthas corticis.

Fig. 2 .
Fig. 2. COI network for all individuals.Coloured circles represent haplotypes present in the Azores; empty circles represent haplotypes present only outside the Azores.Redbrownish indicates Amynthas gracilis haplotypes; Green indicates A. corticis haplotypes.Blue indicates other Megascolecidae found in São Miguel.Names indicate origin of non-Azorean specimens (italics and underlined) or localities in São Miguel for the minority haplotypes.Geographic distribution of main haplotypes is found in Fig. 1.Size of the circles is proportional to the number of haplotypes.Number of SNPs separating haplotypes is shown on the connecting lines.Groups recovered using the less variable gene 16S are grouped within a dashed line.Only Amynthas individuals outside the Azores, and only in three instances (YN for A. gracilis; US and BR for A. corticis) showed different 16S haplotypes.See locality codes in Supplementary Tables1 and 2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Bayesian phylogenies for the mitochondrial markers analysed for the Amynthas dataset.A: COI gene e posterior probabilities (above) and the dates generated by Beast for relevant ancestors originating the defined groups (below) are shown on the nodes.B: concatenated data set, including only A. gracilis and A. corticis.Posterior probabilities are shown on the nodes.Clades identified with a vertical bar represent haplotypes present in Azores.Codes for non-Azorean specimens are shown in italics and underlined.See locality codes in Supplementary Tables1 and 2 and haplotype codes in Supplementary Tables 4, 5 and 7.

Fig. 4 .
Fig. 4. Surface maps generated with SURFER © showing the overlaid levels of abundance of individuals of Amynthas gracilis, and A. corticis derived from 32 collection sites in São Miguel Island, Azores.Geographical coordinates in decimal degrees are represented by x and y and the location of the corresponding sampling points are shown in Supplementary Table1in the same units.

Fig. 5 .
Fig. 5. Surface maps generated with SURFER © showing the overlaid levels of abundance of individuals of Amynthas corticis lineage C, total metal load in soils (mol$kg À1 ), and altitude (m) derived from 32 collection sites in São Miguel island, Azores.Geographical coordinates in decimal degrees are represented by x and y and the location of the corresponding sampling points are shown in Supplementary Table1in the same units.

Table 1
Loadings for axes 1 and 2 according to PCA built with soil characteristics and other habitat features measured in São Miguel Island (Azores).Bold values are considered high !0.7.

Table 2
Results of the final GLM models for the abundance of each earthworm species and lineage, and the lineage diversity that include the significant variables highlighted in previous full models with interactions.Pseudo-R 2 indicates the proportion of variation explained by the model (see Methods).
Table 1 in the same units.