Research Article
Research Article
Population genetics and diversity structure of an invasive earthworm in tropical and temperate pastures from Veracruz, Mexico
expand article infoDiana Ortíz-Gamino§, Josefat Gregorio|, Luis Cunha, Esperanza Martínez-Romero#, Carlos Fragoso¤, Ángel I. Ortíz-Ceballos§
‡ Subdirección de Recursos Naturales y Cambio Climático del H. Ayuntamiento de Xalapa, Xalapa, Mexico
§ Universidad Veracruzana, Xalapa, Mexico
| Instituto Politécnico Nacional, Mexico, Mexico
¶ University of South Wales, Pontypridd, United Kingdom
# Universidad Nacional Autónoma de México Campus Morelos, Cuernavaca, Mexico
¤ Instituto de Ecología, Xalapa, Mexico
Open Access


Pontoscolex corethrurus (Müller, 1857) is an invasive tropical earthworm, globally distributed. It reproduces through parthenogenesis, which theoretically results in low genetic diversity. The analysis of the population structure of P. corethrurus using molecular markers may significantly contribute to understanding the ecology and reproductive system of this earthworm species. This work assessed the genetic diversity and population structure of P. corethrurus with 34 polymorphic inter simple sequence repeat markers, covering four populations in tropical and temperate pastures from Veracruz State. Nuclear markers distinguished two genetic clusters, probably corresponding to two distinct genetic lineages. The number of clones detected in the AC population was lower than expected for a parthenogenetic species. Also, the apparent lack of differences in population structures related to the geographic region among the populations studied may indicate that human-mediated transference is prevalent in these areas. Still, most individuals apparently belong to lineage A, and only a few individuals seem to belong to the lineage B. Thus, the admixture signatures found among the four populations of P. corethrurus may have facilitated a successful invasion by directly increasing fitness. In summary, addressing the genetic variation of P. corethrurus with ISSR markers was a suitable approach, as it evidenced the genetic diversity and relationships in the populations evaluated.


Agroecosystems, asexual reproduction, exotic earthworm, peregrine species, Rhinodrilidae


Earthworms are not only ubiquitous ecological engineers of soil that create biogenic structures; they also sustain the functioning of the ecosystem through their fundamental actions (Lavelle et al. 1997; Brown et al. 2001; Barros et al. 2004). Despite the importance of soil organisms in ecosystem functioning, the soil ecosystem seems to be poorly studied (González et al. 2006). In fact, for a long time, little attention has been given to invasive soil organisms such as earthworms (Gates 1954), despite a rise in belowground invasion over the past 30 years (Hendrix et al. 2002; Craven et al. 2016; Cicconardi et al. 2017). In Mexico, the first records regarding exotic earthworms date from 1900–1906, being lumbricids (Dendrobaena octaedra Savigny, 1826; Lumbricus terrestris Linnaeus, 1758), megascolecids (Metaphire californica Kinberg, 1867), benhaminis (Dichogaster bolaui Michaelsen, 1891), and rhinodrilids (Pontoscolex corethrurus Müller, 1857) the earthworms registered (Eisen 1900; Michaelsen 1900; Beddard 1912). Since then, 51 exotic species have been described across the country by classical taxonomy (Fragoso and Rojas 2014). Several factors acting at different temporal and spatial scales are involved in earthworm invasion, but the overall picture is not yet understood (Marinissen and Van den Bosch 1992; Baker et al. 2006; Brown et al. 2006; González et al. 2006; Cameron and Bayne 2009; Marichal et al. 2012). In general, the genetic characteristics of invasive organisms have profound impacts on their establishment capacity, range expansions, and successful invasion (Lee 2002).

Pontoscolex corethrurus, formerly included in Glossoscolecidae but is now placed in Rhinodrilidae (James 2012), is an earthworm species originating in South America, in the Guiana Shield area of the Amazon (Gates 1954; Righi 1984). The global and local distribution of P. corethrurus have been addressed, as well as its interspecific interactions (Ortíz-Gamino et al. 2016; Taheri et al. 2018). It has been hypothesized that the successfulness of P. corethrurus is based on its genetic plasticity, which in turn is given by high genomic promiscuity associated to its reproductive strategies (Vitturi et al. 2002; Bengtsson 2009; Fernández et al. 2013; Cunha et al. 2014; Pavlíček et al. 2016).

Parthenogenesis is common in earthworms, usually associated with dispersal, where a single propagule is usually sufficient to stablish a new population (Terhivuo and Saura 2006). Thus, rapid adaptation of parthenogenetic clonal populations may be an essential mechanism for a successful colonization event, as is the case of P. corethrurus (Gates 1973; Lavelle et al. 1987; González et al. 2006; Hendrix et al. 2006; Terhivuo and Saura 2006; Buch et al. 2011). Under such a scenario, the genetic variation of P. corethrurus in regions far from its natural geographic range may be low if only a single or very few individuals were colonizers (Dupont et al. 2012). Through mitochondrial and nuclear markers, 792 earthworms collected recently over 25 countries demonstrated that P. corethrurus is a complex of cryptic species. This is represented by a monophyletic clade composed of four morphologically indistinguishable lineages named as L1, L2, L3, L4 (Taheri et al. 2018).

Around the world, P. corethrurus is distributed from 0 to 2000 m a.s.l., with an average altitude of 463 m (Fragoso et al. 1999). In Mexico, specifically in Veracruz State, its distribution ranges from 0 to around 1600 m a.s.l., living at an average temperature of 17 °C (Ortíz-Gamino et al. 2016). The distribution of P. corethrurus in Mexico seems to be strongly associated with human-mediated dispersal due to agricultural activities (González et al. 2006; Feijoo et al. 2007; Hendrix et al. 2008; Dupont et al. 2012). Moreover, although Taheri and co-workers have determined recently that P. corethrurus living in the State of Veracruz correspond to lineages L1 and L3 (Taheri et al. 2018), this finding was based on only a few specimens, which is likely not representing the entire population. For this reason, further research on population genetics may significantly contribute to understanding the ecology of P. corethrurus in Mexico. Thus, the objective of this work was to explore the genetic variation and population structure in P. corethrurus inhabiting tropical and temperate pastures in Mexico using a molecular approach based on ISSR markers.

Materials and methods

Sampling sites and animal collection

Sampling points were established according to the different attributes of the sites studied in Veracruz State, Mexico (Ortíz-Gamino et al. 2016). In brief, the sampling sites were Laguna Verde (LV), Actopan (AC), Ingenio La Concepción (LC) and Naolinco (NA) (Table 1 and Figure 1), each of them with characteristic ecological attributes. During September 2013, 40 mature (clitellate) individuals of P. corethrurus were collected (N = 10 per site). Earthworms were kept in plastic boxes with moistened soil and transported to the laboratory at INBIOTECA for taxonomical and anatomical identification (Moreno and Borges 2004). Specimens were rinsed in water to remove soil particles and were fixed with 96% ethanol. All samples were kept at -20 °C until further processing.

Figure 1. 

Pastures sampled in the central region of Veracruz State, Mexico. LV, Laguna verde; AC, Actopan; LC, La Concepción; NA, Naolinco. The digital elevation model was created using data provided by Instituto Nacional de Estadística y Geografía, México.

Table 1.

Attributes of earthworms sampling of four pastures in central Veracruz State, Mexico.

Sampling Site Altitude (m a.s.l.) Climate Grass species Soil texture (%)
Clay Silt Sand
Laguna Verde (LV) 24 Aw1(w)g Paspalum conjugatum, Cynodon nlemfuensis 26.6 28.1 45.3
Actopan (AC) 480 Aw0(w)gw” Saccharum officinarum L. 12.8 32.3 54.9
La Concepción (LC) 973–1036 (A)Ca(f)gw” Paspalum conjugatum, Cynodon nlemfuensis 26.6 28.2 45.3
Naolinco (NA) 1566–1667 Cb(fm)gw” Paspalum conjugatum, Cynodon nlemfuensis, Pennisetum clandestinum 12.8 32.3 54.9

DNA isolation and quantification

Tail-wall tissue was used for extraction of genomic DNA. Total DNA was extracted using the DNeasy Blood & Tissue kit (Qiagen, Mainz, Germany) following the manufacturer’s instructions. DNA was checked for quality by gel electrophoresis and quantified using a spectrophotometer (ND-2000, Nanodrop Technologies, Wilmington, DE).

Inter Simple Sequence Repeat (ISSR) protocol

Forty specimens of P. corethrurus (N = 10 per site) were used for ISSR screening. ISSR screening was based on five primers (Table 2) previously reported to produce polymorphic and reproducible DNA fingerprints for Eudrilus eugeniae (Kinberg, 1867) and Eisenia fetida (Savigny, 1826) (Sharma et al. 2011). Each PCR reaction contained 1X PCR buffer, 2 mM MgCl2, dNTPs at 100 µM each, primers at 0.8 µM each, 1.5 U of Taq DNA polymerase, 1ug/ul BSA and 30 ng template DNA. The PCR reaction mix was brought to a final volume of 10 µL with water. PCR amplifications were performed in an Eppendorf Mastercycler (Eppendorf, Hamburg, Germany) according to the following conditions: an initial step at 95 °C for 3 minutes, followed by 35 denaturation cycles at 95 °C for 30 seconds, annealing at primer-specific temperature for 30 seconds, and elongation at 72 °C for 1 minute. A final extension was performed at 72 °C for 10 minutes. The PCR products were visualized in 2% agarose gels (with ethidium bromide at 1ul/ml). Although the initial screening used a total of 18 primers, only those primers that were polymorphic and reproducible were selected for subsequent analysis (Table 2).

Table 2.

Primers used for PCR amplification of Pontoscolex corethrurus genomic DNA.

Primer Sequence Ta (°C) Maximum number of bands Estimated size (bp)
840 GAGAGAGAGAGAGAGAYT 59.5 8 2000-200
834 AGAGAGAGAGAGAGAGYT 61 7 2000-300
866 CTCCTCCTCCTCCTCCTC 70 6 2000-400
810 GAGAGAGAGAGAGAGAT 52.4 6 2000-400
807 AGAGAGAGAGAGAGAGT 54.4 7 2000-300

Data analysis

The amplified DNA fragments were transformed into a binary matrix (1 = presence, 0 = absence), as reported previously (Abbot 2001). A multilocus genotype (MLG) was constructed for each individual by pooling data of single ISSR fingerprints using the procedure available in the POPPR package in R for genetic analysis of populations (Kamvar et al. 2014). Isolates with the same MLG were considered clones, and some analyses were conducted for the original and clone-corrected datasets (Gramaje et al. 2014). The POPPR package in R was used to calculate dissimilarity distance matrices and generate a minimum spanning network from these matrices (Kamvar et al. 2014). To assess the potential evolutionary relationships among MLGs, a minimum spanning network was constructed using the genotypes of earthworm from each sampling location. Bootstrapping was performed with 1000 bootstrap resampling. Genotypic diversity, genetic richness and the evenness index were calculated for each population. The ‘rarecurve’ function from the VEGAN package in R (Oksanen et al. 2007) was used to generate rarefaction curves. Curves were calculated to determine whether the sampling intensity was adequate to detect most of the MGLs in each P. corethrurus population. Additionally, the minimum number of loci needed to distinguish all MGLs was also calculated. Given that sample size varied among populations, we employed rarefaction to explore the effect of sample size on observed MLGs.

Genetic differentiation, structure, and clustering analyses

The genetic variance for all MLGs was estimated through an analysis of molecular variance (AMOVA) using the GenAlEx v.6.5 software (Peakall and Smouse 2006). Genetic variance relative to total variance was calculated as PhiPT (analog of the Fst fixation index) for all populations, as well as regarding within-population genetic variance (Peakall and Smouse 2006). The significance was computed by using 9,999 permutations, and the confidence interval at 95%, by 10,000 re-samplings. For this analysis, only single copies of the different genotypes were used to give identical weight to MLGs. The Mantel test was used to explore the potential correlation between the matrix of genetic differentiation between pairs of MLGs and the matrix of spatial distances between populations, using Arlequin v.3.5 (Excoffier and Lischer 2010). The association between P. corethrurus individuals was assessed initially using a Principal Components Analysis (PCA) implemented in GenAlEx v. 6.51 (Peakall and Smouse 2006). As PCA is independent of any genetic hypotheses it is suitable for the analysis of partially clonal species. Additionally, Unweighted Pair Group Method with Arithmetic Mean (UPGMA) dendrograms were also created using the POPPR package (Kamvar et al. 2014). Bootstrapping was performed with the PVCLUST package in R using 10,000 bootstrap re-samplings (Suzuki and Shimodaira 2006).

Population structure was explored using the Bayesian clustering method implemented in STRUCTURE v.2.3.4 (Pritchard et al. 2000), as well as a distance-based approach using a discriminant analysis of principal components (DAPC) (Jombart et al. 2010). STRUCTURE v.2.3.4 was used to identify the number of genetic clusters within the dataset, and to assign individuals to the clusters defined using an admixture model. To this end and to confirm consistency, 15 replicate runs were carried out for each K (1–8). The most likely value of K was determined with “BestK” implemented in CLUMPAK (Evanno et al. 2005), which uses the ΔK method of Evanno et al. (2005). The results from the 15 replicate runs were pooled using CLUMPAK online version (Evanno et al. 2005). Relative dissimilarity distances were calculated according to the index of association (Brown et al. 1980; Smith et al. 1993). The approach returns a distance reflecting a ratio of the number of observed differences to the number of possible differences.

The hybridization status of individuals according to the Bayesian genetic clusters defined in Structure (defined as putative Lineage A and Lineage B) was further investigated using NEWHYBRIDS v1.1 (Anderson and Thompson 2002), which also uses a Bayesian assignment by implementing a multilocus allele frequency model-based approach. This approach clusters together MGLs without a-priori knowledge of parental allele frequencies, and also has the advantage of specifically assuming a mixture of parental and several hybrid classes (F1’s, F2’s, and various backcrosses as B1 and B2 hybrids) to assign them into categories. Individual posterior probabilities belonging to each hybrid category were estimated using the MCMC method in a Bayesian framework using Jeffreys-type priors and a burn-in period of 100000 iterations followed by 50000 sweeps from the posterior distribution sampling (Anderson and Thompson 2002). Linkage disequilibrium as an indication of random mating was calculated and tested for significance with 1,000 randomizations using the POPPR package in R (R core team 2004). The measures of gametic disequilibrium tested were the index of association (IA) (Brown et al. 1980; Smith et al. 1993) and a standardized alternative of the IA (r̄d) (Agapow and Burt 2001). The null hypothesis for this test is that there is a random association among alleles at different loci and IA = 0; the null hypothesis for random mating is rejected where if IA>0.


Beyond the ecological relevance of P. corethrurus, information on genetic variability is relevant to determine the selective forces that act on the reproductive system of this species. In that sense, the survey carried out in this study was aimed to reveal the population genetic structure of P. corethrurus in natural landscapes covering tropical and temperate pastures. For this, a set of ISSR markers were used to assess the genetic diversity and population structure of P. corethrurus genotypes from four locations in Veracruz State, Mexico.

Genotypic diversity of P. corethrurus populations

Following a PCR-based approach, ISSR primers produced bands on agarose gel that were suitable for assessing the genetic diversity and genetic relationship between and across populations of P. corethrurus. The PCR products ranged between ~200 and ~2000 bp from genomic DNAs of P. corethrurus. The total number of bands and polymorphism rates are shown in Table 3. A total of 33 MLGs among the 35 P. corethrurus individuals yielded reliable products. Overall, one MLG corresponding to AC was observed twice, whereas the rest of MLGs (31) were detected only once. As regards MLG diversity (H), this parameter varied across populations, with no correlation to any specific geographic location (Table 3). On the other hand, in contrast to AC, evenness values (E5) were higher for LV, LC, and NA, in agreement with the fact that all genotypes found were unique to these populations (Table 3). Nei’s unbiased gene diversity (Hexp) values varied from the highest (LC = 0.39) to the lowest (NA = 0.29). Most populations exhibit low genetic diversity (Table 3), except for the NA population, which displays a higher genetic diversity. This was supported by the rarefaction curves, which indicated that NA had a higher number of sampled loci, as well as higher MGLs, than the other three populations (Figure 2A). Interestingly, the minimum number of loci needed to define the total number of MLGs found reached a plateau after 18 loci (Figure 2B). Moreover, the variation among and within populations assessed with AMOVA resulted in values of 25% and 75%, respectively (Table 4). Altogether, all populations showed high genotypic diversity, and according to the PhiPT value (analog of Fst fixation index), there were significant differences between the LV, AC, and NA populations were found (Table 4).

Figure 2. 

Rarefaction curve of expected number of MLGs captured per earthworm of Pontoscolex corethrurus sampled (A), and a MLG accumulation curve according to the number of loci sampled (B).

Table 3.

Parameters of genetic variation in four Pontoscolex corethrurus populations living in central Veracruz State, Mexico.

Sampling Site N MLG eMLG Pb Tb H G E 5 Hexp Ia rbarD
LV 9 9 9 20 34 2.20 9 1 0.30 1.98* 0.10*
AC 10 8 8 22 34 2.03 7.14 0.93 0.30 3.23* 0.15*
LC 6 6 6 27 34 1.79 6 1 0.39 1.16* 0.04*
NA 10 10 10 24 34 2.30 10 1 0.29 4.08* 0.18*
Total 35 33 9.85 3.48 31.41 0.97 0.40 1.65 0.05
Table 4.

Analysis of molecular variance (AMOVA) testing for genetic differentiation between four populations of Pontoscolex corethrurus living in central Veracruz State, Mexico (A), and PhiPT pairwise comparisons (B).

Level of variation d.f. SS MS Est. Var. Proportion (%)
A) Among Populations 3 59.079 19.693 1.760 25%
Within Populations 29 154.497 5.327 5.327 75%
Total 32 213.576 7.088 100%
Population Laguna Verde Actopan La Concepcion Naolinco
B) Laguna Verde 0.000
Actopan 0.229** 0.000
La Concepcion 0.069 0.161** 0.000
Naolinco 0.335** 0.347** 0.199* 0.000

Clustering of MLGs: relationships between- and within P. corethrurus populations

Although some individuals cluster together according to site (e.g., animals in NA), most of the individuals scattered in a non-uniform clustering (Figure 3A). As shown in Figure 3A, axes 1 and 2 of the PCA scatter plot accounted for 26% and 15% of total genetic variability, respectively. On the other hand, the global minimum spanning network showed that all populations have MLGs that are closely related (Figure 3B). The comparison of the matrix of Euclidian genetic distance with the matrix of geographic distances using the Mantel test showed that there is no correlation between these two matrices. Thus, data in the genetic distance matrix is not explained by the geographic positioning of the populations (Figure 3). Moreover, the genetic distance between MLGs and between populations indicate no evident correlation with geographic locations, even though LC and LV appeared to be genetically related (Figure 4A, B, respectively). In summary, the clustering analysis shows no clear association with geographic distances.

Figure 3. 

A Principal Components Analysis, where colors indicate specimens of the population (A) and a Minimum Spanning Network where each node denotes a different MLG, with size matching the number of individuals. Edge thickness and color are proportional to absolute genetic distance. Edge lengths are arbitrary (B). Both analyses show the relationship between multilocus genotypes (MLGs) for four different earthworm populations of Pontoscolex corethrurus living in central Veracruz State, Mexico.

Figure 4. 

UPGMA dendrogram of genetic distance between MGLs (A) and between populations (B) observed in the distinct populations of Pontoscolex corethrurus collected in central Veracruz State, Mexico. Only bootstrap values higher than or equal to 70% are shown.

Estimation of population structure according to genetic clusters

The Bayesian analysis of population structure estimated two distinct genetic clusters (K = 2, Ln P (D) = -557.89 ± 0.31) distributed across the four geographic locations (Figure 5). Similar to the results obtained by the PCA and dendrogram analyses, the Bayesian analysis revealed no apparent structure that could be associated with geographic location (Figure 5). Nevertheless, the NA population seems to belong mostly to cluster 2 and the one from LV to cluster 1, whereas AC and LC appear to be strongly admixed (Figure 5). From this analysis, two distinct genetic lineages can be identified among populations, henceforth defined as lineages A (mostly LV individuals) and lineage B (mostly NA individuals). Besides, the DAPC analysis confirmed the above genetic relationship between LC and LV, even when they belong to distant geographic locations (Figure 6). On the other hand, similar to the results of the Bayesian analysis, AC (lineage A) and NA (lineage B) are detached (Figure 7). Notably, both clone-corrected (N = 33) and uncorrected (N = 35) reject the hypothesis of no linkage among markers (Table 3), supporting asexuality in all populations (Suppl. material 1: Figure S1). Altogether, the assessment of genetic diversity in the four populations of P. corethrurus suggests that they belong to two different lineages, with some relationships among them despite their distribution in different locations.

Figure 5. 

Estimated population genetic structure with a summary plot of Q estimates based on the ISSR data observed for four populations of Pontoscolex corethrurus in central Veracruz State, Mexico. Each individual is shown by a vertical line, which is partitioned into colored segments representing the fraction of the number of members in cluster K (%).

Figure 6. 

Genetic structure using ISSR data for 35 Pontoscolex corethrurus individuals based on discriminant analysis of principal components (DAPC). Proportion of eigenvalues in discriminant analysis (bottom left plot) and PCA eigenvalues (bottom right), with the first 12 significant principal components highlighted in black.

Figure 7. 

Classification of Pontoscolex corethrurus individuals according to a Bayesian assignment algorithm implemented in NEWHYBRIDS (Anderson and Thompson 2002) to detect gene flow. Each unit represents an individual corresponding to parental lineages (Lineage A and Lineage B), F1 generation, F2 (F1 x F1) and later generation or introgressive hybrids B1 (Lineage A x F1) and B2 (e.g., Lineage B x F1).


Soil texture and chemical composition are among the key environmental variables that play a role in the invasion process of earthworms (Hendrix and Bohlen 2006; Marichal et al. 2012; Craven et al. 2016). On the other hand, traits intrinsic to this species, such as dispersal capacity, life history, and reproduction mode, are known to influence the genetic structure of earthworms (González et al. 2006; Cameron et al. 2008; Marichal et al. 2012). All these features can contribute to the isolation of populations, favoring intra-specific relationships that determine the population structure (Baker 1998; Novo et al. 2010; Fernández et al. 2013). Thus, determining the genetic variability in asexual organisms, such as the earthworm P. corethrurus, is of interest, since it may influence population structure and therefore, putative evolutionary bottlenecks.

In this work, the study of genetic diversity and population structure of four populations living in Mexican tropical and temperate pasture was carried out through a molecular approach using ISSR markers. The use of ISSR markers is well supported by several studies since it produces a high percentage of polymorphic loci (Abbot 2001; Dušinský et al. 2006; Luan et al. 2006; Sharma et al. 2011; Seyedimoradi and Talebi 2014; Ng and Tan 2015; Buczkowska et al. 2016). Also, ISSR markers potentially discriminate isolated populations by geographic conditions and are able to differentiate cryptic species ((Dušinský et al. 2006; Buczkowska et al. 2016). Finally, since ISSR markers only amplify nuclear regions of eukaryotic genomes, it avoids the amplification of bacterial DNA fragments, i.e., members of the phylum Bacteroidetes, commonly associated with P. corethrurus (Bernard et al. 2012; Seyedimoradi and Talebi 2014).

According to PhiPT values, there were significant differences between the LV, AC, and NA populations, meaning that most populations possess high genotypic diversity (Figure 2 and Table 4). In this regard, the AMOVA analyses showed that most genetic diversity (75%) was found within populations, and only 25% among populations. This finding is similar to those reported by Cameron et al. (2008) and Cunha et al. (2014) for Dendrobaena octaedra and P. corethrurus, respectively. In both cases, most of the genetic variability was found within populations, and only a minor proportion occurred between populations (Cameron et al. 2008; Cunha et al. 2014). The high levels of genotypic diversity observed in this study could be expected because there is evidence of high diversity in populations of P. corethrurus using AFLPs (Cunha et al. 2014). It is also likely that the high genotypic diversity observed in the populations of P. corethrurus studied is due to the high evolutionary rate of change within ISSR regions compared to other nuclear regions tackled by AFLPs, RAPD or URPs (Sharma et al. 2011; Gramaje et al. 2014; Seyedimoradi and Talebi 2014; Ng and Tan 2015).

As regards the movement of MLGs across the different sites, a conspicuous behavior was observed, particularly in LV, AC, and LC. Such observation may be explained only by an intense human-mediated transference of P. corethrurus through road networks or activities associated with agriculture (Lapied and Lavelle 2003; Baker et al. 2006; González et al. 2006; Feijoo et al. 2007; Dupont et al. 2012; Ortíz-Gamino et al. 2016). On the other hand, clustering analyses through PCA and a dendrogram identified a significant admixture across all populations; however, it was possible to identify at least two divergent and well-differentiated genetic clusters (lineages A and B). It is likely that the lineages identified potentially correspond to those previously reported by Taheri and co-workers, namely lineages L1 and L3 (Taheri et al. 2018). It is important to highlight that both our sampling sites and the molecular approach, are different from the work done by Taheri et al. (2018). Their research involved specimens collected in different years, covering an extended period (from 1996 to 2016; Taheri et al. 2018, Suppl. material 1: Figure S1). Despite these differences, our results are consistent with those reported by Taheri et al. (2018), namely, the identification of two lineages in the P. corethrurus populations.

Lineage A seems to be widespread, covering LC and LV sites (Figure 6). Such distribution may indicate wide ecological tolerance, with populations probably well adapted to warm temperatures and poor soils. The distribution of lineage A may suggest that this lineage corresponds to L1 (Taheri et al. 2018), since it was found in most of the places sampled. Moreover, our results are similar to those reported for other species, such as Octolasion tyrtaeum (Savigny, 1826), for which a single haplotype was found in all sampling stations (Terhivuo and Saura 1993). Hence, for peregrine species, the evidence indicates that human activities are strongly shaping the dispersal pattern of MLGs through incidental transfer in crops soil or fishing (Cameron et al. 2008; Dupont et al. 2012; Ortíz-Gamino et al. 2016). However, there are also reports of the intentional introduction of earthworm species for commercial applications like waste management and land bioremediation (Hendrix 2006).

In contrast to lineage A, lineage B (mostly NA specimens) showed the best distinct cluster of individuals (Figure 5). This cluster is likely associated with the contrasting environmental conditions that predominate in NA, namely higher altitude, three types of grass (Paspalum conjugatum P.J. Bergius, Cynodon nlemfuensis Vanderyst, Pennisetum clandestinum Hochst. ex Chiov.), soil rich in organic matter, or even interactions with other soil organisms (Ortíz-Gamino et al. 2016). All these characteristics could be acting as an environmental screen that results in the clustering of NA individuals (Figure 3). Remarkably, a similar finding has been reported for Aporrectodea trapezoides (Dugés, 1828), in which clonal lineages seem to remain close to their original areas, indicative of some level of local adaptation or strong interspecific relationships (Fernández et al. 2013). Thus, the lineage described as L2 by Taheri et al. (2018) could correspond to lineage B in this work. If so, this lineage is likely well-adapted to micro conditions including habitat, feeding habits and biotic interactions. As proposed earlier, another barrier to the dispersal of lineage B may be temperature (Janzen 1967). This could be the case for the NA population, in which the mean annual temperature (17 °C) may be acting as the main barrier to MGL dispersion (Ortíz-Gamino et al. 2016). Importantly, under a global-change scenario, namely intensive land-use change and alarming global warming, this barrier could become weaker, thus enabling the invasion of pantropical earthworm species (Jiménez and Decaëns 2000; Eisenhauer et al. 2014; Gutiérrez and Cardona 2014).

On the other hand, sexual reproduction is a rare event in P. corethrurus (Gates 1973), as it is widely accepted that its reproduction occurs mainly by parthenogenesis (Gates 1973). The standardized index of association (r̄ d) supported the hypothesis of clonal population structure (Suppl. material 1: Figure S1). In this sense, r̄ d values suggest a widespread dispersal of MLGs across populations. This contrasts with the linkage disequilibrium tests, in which the null hypothesis of random mating was rejected for all populations. Nevertheless, these results should be interpreted with caution, as it is challenging to demonstrate the presence of linkage disequilibrium with small sample sizes (Hagenblad et al. 2006; Du et al. 2007; Gramaje et al. 2014). In earthworms, parthenogenesis has been associated with polyploidy, as well as with high levels of DNA methylation (Regev et al. 1998). Therefore, it is also plausible that methylation may be fostering the epigenetic control of phenotypic plasticity, which could be crucial for a successful colonization (Stürzenbaum et al. 2009). It is tempting to claim that temperature affects P. corethrurus and impacts its reproduction rate, which also may be regulated by polyploidy and epigenetic control (Ortíz-Gamino et al. 2016). Further studies regarding the number of chromosomes or genomic rearrangements are needed to address whether or not these features are linked to environmental features.

In summary, the screening of genetic diversity is helpful to monitor the dynamics of population structure and its relationships to ecological and environmental features and to contribute valuable information about the isolation of invasive earthworm species. In this sense, our work provides evidence of the existence of two lineages of P. corethrurus in Veracruz State, Mexico, showing different distribution patterns according to the prevailing environmental conditions found in regions studied. Therefore, our data represent an relevant contribution to know the movement dynamics and diversification of P. corethrurus, which will be useful information for planning successful strategies aimed to control or prevent the biological invasion of this species in Mexico.


Despite being random, biological invasions are intriguing events, mainly because they involve populations of organisms with certain features and particular habits. Intriguingly, a parthenogenic species such as P. corethrurus has been successful in colonizing areas all over the world. The interaction of P. corethrurus genetics with the environment should drive its selection and distribution pattern.

In this work, we assessed populations of P. corethrurus inhabiting tropical and temperate pastures of Veracruz, Mexico, in terms of genetic variation through ISSR markers. Our results revealed the existence of at least two well-differentiated genetic clusters, corresponding to different lineages (lineages A and B). The lineages identified in this work likely correspond to lineages L1 and L3 identified previously by Taheri et al. (2018). Although further research is needed to discern why P. corethrurus populations occur in some sites but not in others, our work suggests that genetic variation is playing a key role in the invasion process. The association between the genetic variability of P. corethrurus and its success in invading new sites is counter-intuitive due to its parthenogenic reproduction, i.e., clonal multiplication. Additional genotyping of P. corethrurus individuals inhabiting diverse environments or different States like Tabasco, Puebla, and Tamaulipas, will be necessary, not only to confirm our results, but to track the dynamics of dispersal and diversification of lineages, as well as to identify new dominant genotypes or newly introduced lineages. All this information should be gathered before using P. corethrurus in biotechnology research, remediation, or fishing, or before estimating its effects on local crops, plants, or organisms.


We would like to thank to Norma Flores-Estévez and Wilberth A. Poot-Poot for the support and help at INBIOTECA, Ernesto Ormeño-Orrillo for his help at Ciencias Genómicas, and to Rogelio Lara González and Katia J. Martínez Velázquez for their help during fieldwork. We are thankful to Alejandra Velázquez Lobo for valuable comments and careful revision of the manuscript and Alberto Ortíz Brito for proving the sampling map. Diana Ortíz-Gamino was supported by a scholarship (No. 251818) by Consejo Nacional de Ciencia y Tecnología (CONACYT).


  • Abbot P (2001) Individual and population variation in invertebrates revealed by Inter-simple Sequence Repeats (ISSRs). Journal of Insect Science 1: 1–8.
  • Anderson EC, Thompson EA (2002) A model-based method for identifying species hybrids using multilocus genetic data. Genetics 160(3): 1217–1229.
  • Baker GH (1998) The ecology, management, and benefits of earthworms in agricultural soils, with particular reference to southern Australia. In: Edwards CA (Ed.) Earthworm Ecology. St Lucie Press, Boca Raton, 1229–1257.
  • Baker GH, Brown G, Butt K, Curry JP, Scullion J (2006) Introduced earthworms in agricultural and reclaimed land: their ecology and influences on soil properties, plant production and other soil biota. Biological Invasions 8: 1301–1316.
  • Barros E, Grimaldi M, Sarrazin M, Chauvel A, Mitja D, Desjardins T, Lavelle P (2004) Soil physical degradation and changes in macrofaunal communities in Central Amazon. Applied Soil Ecology 26: 157–168.
  • Bernard L, Chapuis-Lardy L, Razafimbelo T, Razafindrakoto M, Pablo AL, Legname E, Brüls T, O’Donohue M, Brauman A, Chotte J-L, Blanchart E (2012) Endogeic earthworms shape bacterial functional communities and affect organic matter mineralization in a tropical soil. The ISME Journal 6(1): 213–222.
  • Bengtsson BO (2009) Asex and Evolution: A very Large-Scale Overview. In: Schön I, Martens K, Van Dijk P (Eds) Lost Sex: The Evolutionary Biology of Parthenogenesis. Springer, New York, 615 pp.
  • Brown GG, Fragoso C, Barois I, Rojas P, Patrón JC, Bueno J, Moreno AG, Lavelle P, Ordaz V, Rodríguez C (2001) Diversidad y rol funcional de la macrofauna edáfica en los ecosistemas tropicales mexicanos. Acta Zoológica Mexicana (n.s. ) 1: 79–110.
  • Brown AHD, Feldman MW, Nevo E (1980) Multilocus structure of natural populations of Hordeum spontaneum. Genetics 96(2): 523–536.
  • Brown GG, James SW, Pasini A, Nunes DH, Benito NP, Trigo P, Sautter KD (2006) Exotic, Peregrine, and Invasive Earthworms in Brazil: Diversity, Distribution, and Effects on Soils and Plants. Caribbean Journal of Science 42: 339–358.
  • Brown GG, Moreno A, Barois I, Fragoso C, Rojas P, Hernández B, Patrón JC (2004) Soil macrofauna in SE Mexican pastures and the effect of conversion from native to introduced pastures. Agriculture, Ecosystems & Environment 103(2): 313–327.
  • Buch A, Brown GG, Niva CC, Sautter KD, Lourençato LF (2011) . Life cycle of Pontoscolex corethrurus (Müller, 1857) in tropical artificial soil. Pedobiologia 54 Supplement: S19–S25.
  • Buczkowska K, Rabska M, Gonera P, Pawlaczyk EM, Wawrzyniak P, Czolpinska M, Bączkiewicz A (2016) Effectiveness of ISSR markers for determination of the Aneura pinguis cryptic species and Aneura maxima. Biochemical Systematics and Ecology 68: 27–35.
  • Cameron E, Bayne EM, Coltman DW, Coltman DW (2008) Genetic structure of invasive earthworms Dendrobaena octaedra in the boreal forest of Alberta: insights into introduction mechanisms. Molecular Ecology 17: 1189–1197.
  • Cicconardi F, Borges PA, Strasberg D, Oromí P, López P, Pérez-Delgado H, Casquet AJ, Caujapé-Castells J, Fernández-Palacios J, Thébaud B, Emerson BC (2016) MtDNA metagenomics reveals large-scale invasion of belowground arthropod communities by introduced species. Molecular Ecology.
  • Cuatrecasas J (1964) Cacao and its allies: a taxonomic revision of the genus Theobroma. Contributions from the United States National Herbarium 35: 379–614.
  • Cunha L, Montiel R, Novo M, Orozco-terWengel P, Rodrigues A, Morgan AJ, Kille P (2014) Living on a volcano’s edge: genetic isolation of an extremophile terrestrial metazoan. Heredity 112: 132–142.
  • Craven D, Thakur MP, Cameron EK, Frelich LE, Beauséjour R, Blair RB, Blossey B, Burtis J, Choi A, Dávalos A, Fahey TJ, Fisichelli NA, Gibson K, Handa IT, Hopfensperger K, Loss SR, Nuzzo V, Maerz JC, Sackett T, Scharenbroch BC, Smith SM, Vellend M, Umek LG, Eisenhauer N (2016) The unseen invaders: introduced earthworms as drivers of change in plant communities in North American forests (a meta-analysis). Global Change Biology 23(3): 1065–1074.
  • Du FX, Clutter AC, Lohuis MM (2007) Characterizing Linkage Disequilibrium in Pig Populations. International Journal Biological Sciences 3(3): 166–178.
  • Dupont L, Decaëns T, Lapied E, Chassany V, Marichal R, Dubs F, Maillot M, Ro V (2012) Genetic signature of accidental transfer of the peregrine earthworm Pontoscolex corethrurus (Clitellata, Glossoscolecidae) in French Guiana. European Journal of Soil Biology 53: 70–75.
  • Dušinský R, Kúdela M, Stloukalová V, Jedlička L (2006) Use of inter-simple sequence repeat (ISSR) markers for discrimination between and within species of blackflies (Diptera: Simuliidae). Biologia 61(3): 299–304.
  • Eisenhauer N, Stefanski A, Fisichelli NA, Rice K, Rich R, Reich PB (2015) Warming shifts ‘worming’: effects of experimental warming on invasive earthworms in northern North America. Scientific Reports 4: 1–6890.
  • Eisen G (1900) Researches in the American Oligochaeta, with special reference to those of the Pacific coast and adjacent islands. Proceedings of California Academy of Science 2: 85–276.
  • Excoffier L, Lischer HEL (2010) Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources 10(3): 564–567.
  • Feijoo A, Zúñiga M, Quintero H, Lavelle P (2007) Relaciones entre el uso de la tierra y las comunidades de lombrices en la cuenca del río La Vieja, Colombia. Pastos y Forrajes 30(2): 235–249.
  • Fernández R, Almodóvar A, Novo M, Gutiérrez M, Díaz Cosín DJ (2011) A vagrant clone in a peregrine species: Phylogeography, high clonal diversity and geographical distribution in the earthworm Aporrectodea trapezoides (Dugès, 1828). Soil Biology & Biochemistry 43(10): 2085–2093.
  • Fernández R, Almodóvar A, Novo M, Simancas B, Díaz Cosín DJ (2012) Adding complexity to the complex: New insights into the phylogeny, diversification and origin of parthenogenesis in the Aporrectodea caliginosa species complex (Oligochaeta, Lumbricidae). Molecular Phylogenetics and Evolution 64(2): 368–379.
  • Fernández R, Almodóvar A, Novo M, Gutiérrez M, Díaz Cosin DJ (2013) Earthworms, good indicators for palaeogeographical studies? Testing the genetic structure and demographic history in the peregrine earthworm Aporrectodea trapezoides (Dugès, 1828) in southern Europe. Soil Biology & Biochemistry 58: 127–135.
  • Fragoso C, Coria-Martínez ML, Camarena LM (2009) An update of the earthworm fauna of Los Tuxtlas Ver. and adjacent regions: are native species in the risk of extinction? In: Barois I, Huising EJ, Okoth P, Trejo D, De los Santos M (Eds) Below-Ground Biodiversity in Sierra San Marta, Los Tuxtlas, Veracruz, México. Instituto de Ecología, A.C., Xalapa, Mexico, 219–228.
  • Fragoso C, Rojas P (2014) Biodiversidad de lombrices de tierra (Annelida: Oligochaeta: Crassiclitellata). Revista Mexicana de Biodiversidad, Supl. 85: S197–S207.
  • Gates GE (1954) Exotic earthworms of the United States. Bulletin of the Museum of Comparative Zoology at Harvard College 111: 217–258.
  • Gates GE (1973) Contributions to North American Earthworms (Annelida) No. 6 Contributions to a revision of the earthworm family Glossoscolecidae. I. Pontoscolex corethrurus (Müller, 1857). Bulletin of Tall Timbers Research Station 14: 1–12.
  • Gramaje D, León M, Santana M, Crous PW, Armengol J (2014) Multilocus ISSR Markers Reveal Two Major Genetic Groups in Spanish and South African Populations of the Grapevine Fungal Pathogen Cadophora luteo-olivacea. PLoS ONE 9(10): e110417.
  • González G, Huang CY, Zou X, Rodríguez C (2006) Earthworm invasions in the tropics. In: Hendrit PF (Ed.) Biological Invasions Belowground: Earthworms as Invasions Species. Springer, Dordrecht, 47–56.
  • Hagenblad J, Bechsgaard J, Charlesworth D (2006) Linkage Disequilibrium Between Incompatibility Locus Region Genes in the Plant Arabidopsis lyrata. Genetics 173(2): 1057–1073.
  • Hendrix PF (2006) Biological invasions belowground-earthworms as invasive species. In: Hendrix PF (Ed.) Biological Invasions Belowground: Earthworms as Invasive Species. Springer, Netherlands, 4 pp.
  • Hendrix PF, Bohlen PJ (2002) Exotic Earthworm Invasions in North America: Ecological and Policy Implications Expanding global commerce may be increasing the likelihood of exotic earthworm invasions, which could have negative implications for soil processes, other animal and plant species, and importation of certain pathogens. Bioscience 52(9): 801–811.[0801:EEIINA]2.0.CO;2
  • Hendrix P, Baker G, Callaham M, Damof G, Fragoso C, Gonzalez G, James S, Lachnicht S, Winsome T, Zou X (2006) Invasion of exotic earthworm into inhabited by native earthworm. Biological Invasions 8: 1287–1300.
  • Hendrix PF, Callaham Jr. MA, Drake JM, Huang ChY, James SW, Snyder BA, Zhang W (2008) Pandora’s box contained bait: the global problem of introduced earthworm. Annual Review of Ecology Evolution, and Systematics 39: 593–613.
  • Jiménez JJ, Decaëns T (2000) Vertical distribution of earthworms in grassland soils of the Colombian Llanos. Biology and Fertility of Soils 32: 463–473.
  • Jombart T, Devillard S, Balloux F (2010) Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genetics 11: 1–94.
  • Kamvar ZN, Tabima JF, Grünwald NJ (2014) Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ 2: e281.
  • Lavelle P, Barois I, Cruz I, Fragoso C, Hernández A, Pineda A, Rangel P (1987) Adaptative strategies of Pontoscolex corethrurus (Glossoscolecidae, Oligocheta), a peregrine geophagous earthworm of the humid tropics. Biology and Fertility of Soils 5: 188–194.
  • Lavelle P, Bignell D, Lepage M, Wolters V, Roger P, Ineson P, Heal OW, Dhillion S (1997) Soil Function in a changing world: the role of invertebrate ecosystem engineers. European Journal of Soil Science 33: 159–193.
  • Luan S, Chiang T-Y, Gong X (2006) High Genetic Diversity vs. Low Genetic Differentiation in Nouelia insignis (Asteraceae), a Narrowly Distributed and Endemic Species in China, Revealed by ISSR Fingerprinting. Annals of Botany 98(3): 583–589.
  • Narayanan SP, Sathrumithra S, Christopher G, Thomas AP, Julka JM (2016) Current distribution of the invasive earthworm Pontoscolex corethrurus (Müller, 1857) after a century of its first report from Kerala state, India. Opuscula Zoologica (Budapest) 47(1): 101–107.
  • Ng WL, Tan SG (2015) Inter-Simple Sequence Repeat (ISSR) Markers: Are We Doing It Right? ASM Science Journal 9(1): 30–39.
  • Marichal R, Grimaldi M, Mathieu J, Brown GG, Desjardins T, Lopes da Silva Jr. M, Praxedes C, Martins MB, Velasquez E (2012) Is invasion of deforested Amazonia by earthworm Pontoscolex corethrurus driven by soil texture and chemical properties? Pedobiologia 55(5): 233–240.
  • Mantel N (1967) The detection of disease clustering and a generalized regression approach. Cancer Research 27(2): 209–220.
  • Moreno AG, Borges S (2004) Avances en taxonomía de lombrices de tierra (Annelida: Oligochaeta). Editorial Complutense. Madrid, 428 pp.
  • Nei M (1978) Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics 89(3): 583–590.
  • Novo M, Almodóvar A, Fernández R, Trigo D, Díaz Cosín DJ (2010) Cryptic speciation of hormogastrid earthworms revealed by mitocondrial and nuclear data. Molecular Phylogenetics and Evolution 56: 507–512.
  • Ortíz-Gamino D, Pérez-Rodríguez P, Ortiz-Ceballos AI (2016) Invasion of the tropical earthworm Pontoscolex corethrurus (Rhinodrilidae, Oligochaeta) in temperate grasslands. PeerJ 4: e2572.
  • Oksanen J, Kindt R, Legendre P, O’Hara RB, Simpson GL, Solymos P, Henry M, Stevens MHH, Szoecs E, Wagner H (2007) The vegan package. Community Ecology Package 10: 631–637.
  • Pavlíček T, Cohen T, Yadav S, Glasstetter M, Král P, Pearlson O (2016) Aneuploidy Occurrence in Oligochaeta. Ecology and Evolutionary Biology 1(3): 57–63.
  • R Core Team (2014) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. Available:
  • Schultes RE (1984) Amazonian cultigens and their northward and westward migrations in preColumbian times. In: Stone D (Ed.) Pre-Columbian Plant Migration. Cambridge: Papers of the Peabody Museum of Archaeology and Ethnology (Vol. 76). Harvard University Press, 19–38.
  • Seyedimoradi H, Talebi R (2014) Detecting DNA polymorphism and genetic diversity in Lentil (Lens culinaris Medik.) germplasm: comparison of ISSR and DAMD marker. Physiology and Molecular Biology of Plants 20(4): 495–500.
  • Sharma A, Sonah H, Deshmukh RK, Gupta NK, Singh NK, Sharma TR (2011) Analysis of genetic diversity in earthworms using DNA markers. Zoological Science 28: 25–31.
  • Stürzenbaum SR, Andre J, Kille P, Morgan AJ (2009) Earthworm genomes, genes and proteins: the (re)discovery of Darwin’s worms. Proceeding of the Royal Society B. 276: 789–797.
  • Taheri S, James S, Roy V, Decaëns T, Williams BW, Anderson F, Rougerie R, Chang CH, Brown G, Cunha L, Stanton DWG, Da Silva E, Chen JH, Lemmon AR, Moriarty Lemmon E, Bartz M, Baretta D, Barois I, Lapied E, Coulis M, Dupont L (2018) Complex taxonomy of the ‘brush tail’ peregrine earthworm Pontoscolex corethrurus. Molecular Phylogenetics and Evolution 124: 60–70.
  • Tamura K, Stecher G, Peterson D, Filipski A, Kumar S (2013) MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Molecular Biology and Evolution 30: 2725–2729.
  • Terhivuo J, Saura A (1993) Clonal and morphological variation in marginal populations of parthenogenetic earthworms Octolasion tyrtaeum and O. cyaneum (Oligochaeta, Lumbricidae) from eastern Fennoscandia. Bolletino di Zoologia 60: 87–96.
  • Uribe S, Huerta E, Geissen V, Mendoza M, Godoy R, Jarquin A (2012) Pontoscolex corethrurus (Annelida: Oligocheta) indicador de la calidad del suelo en sitios de Eucalyptus grandis (Myrtaceae) con manejo tumba y quema. Revista de Biología Tropical 60: 1543–552.
  • Vitturi R, Colomba MS, Pirrone AM, Mandrioli M (2002) rDNA (18S–28S and 5S) colocalization and linkage between ribosomal genes and (TTAGGG)(n) telomeric sequence in the earthworm, Octodrilus complanatus (Annelida: Oligochaeta: Lumbricidae), revealed by single-and double-color FISH. Journal of Heredity 93(4): 279–282.

Supplementary material

Supplementary material 1 

Figure S1

Diana Ortíz-Gamino, Josefat Gregorio, Luis Cunha, Esperanza Martínez-Romero, Carlos Fragoso, Ángel I. Ortíz-Ceballos

Data type: image

Explanation note: Linkage disequilibrium test using r̄ d as implemented in R package Poppr (Kamvar et al. 2014). Visualizations of tests for linkage disequilibrium, where observed values (blue dashed lines) of r̄ d are compared to histograms showing results of 999 permutations..

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (55.55 kb)