Diversity of introduced terrestrial flatworms in the Iberian Peninsula: a cautionary tale

Many tropical terrestrial planarians (Platyhelminthes, Geoplanidae) have been introduced around the globe. One of these species is known to cause significant decline in earthworm populations, resulting in a reduction of ecological functions that earthworms provide. Flatworms, additionally, are a potential risk to other species that have the same dietary needs. Hence, the planarian invasion might cause significant economic losses in agriculture and damage to the ecosystem. In the Iberian Peninsula only Bipalium kewense Moseley, 1878 had been cited till 2007. From that year on, four more species have been cited, and several reports of the presence of these animals in particular gardens have been received. In the present study we have: (1) analyzed the animals sent by non-specialists and also the presence of terrestrial planarians in plant nurseries and garden centers; (2) identified their species through morphological and phylogenetic molecular analyses, including representatives of their areas of origin; (3) revised their dietary sources and (4) used Species Distribution Modeling (SDM) for one species to evaluate the risk of its introduction to natural areas. The results have shown the presence of at least ten species of alien terrestrial planarians, from all its phylogenetic range. International plant trade is the source of these animals, and many garden centers are acting as reservoirs. Also, landscape restoration to reintroduce autochthonous plants has facilitated their introduction close to natural forests and agricultural fields. In conclusion, there is a need to take measures on plant trade and to have special care in the treatment of restored habitats.


INTRODUCTION
Most animal invasive species detected in Europe are terrestrial invertebrates (Roques et al., 2009). Invading edaphic organisms can have dramatic effects on the environment, due to the direct effects on native soil organisms, and through their interactions with the environment aboveground. However, overall, their impact in human health and economy is greater than their ecological impact (Vilà et al., 2010). Among these organisms, land planarians are becoming an important and diversified group of introduced species in Europe.
Terrestrial planarians (Platyhelminthes, Geoplanidae) are divided into four subfamilies (Bipaliinae, Microplaninae, Geoplaninae and Rhynchodeminae) with a cosmopolitan distribution (Winsor, Johns & Yeates, 1998); however, most species are found in the southern hemisphere. Bipaliinae (Fig. 1A) is absent from the American and European continents, Geoplaninae ( Fig. 1B) have an exclusively Central and South American distribution, while Microplaninae ( Fig. 2A) and Rhynchodeminae (Fig. 2B) are the subfamilies with the most northerly distribution, including Europe. Terrestrial planarians are the only free-living Platyhelminthes that do not live in an aquatic habitat. However, they have not developed the capacity to prevent water loss and are thus strongly dependent on environmental moisture levels (Froehlich, 1956;McDonald & Jones, 2007). They seem to withstand this limitation through behavioral strategies such as hiding in damp refuges during the day and becoming active during the night. Due to these characteristics, these animals are considered to have a low capacity to disperse. In fact, in their areas of origin, although a few species are well-adapted to open and human-transformed lands (Baptista & Leal-Zanchet, 2010), most species are restricted to humid forest areas.
A total of 36 species of terrestrial planarians are known to have been introduced in different countries around the globe. Most of these species have a big effect on terrestriaĺ ecosystem processes because they prey on soil invertebrates (see references in Winsor, Johns & Barker, 2004). So far, five of these species are considered to be either invasive and cause problems with local biodiversity (Platydemus manokwari De Beauchamp, 1963), or horticultural pests (Arthurdendyus triangulatus (Dendy, 1894)) or earthworm farm pests (Bipalium adventitium Hyman, 1943; Bipalium kewense Moseley, 1878; Dolichoplana striata Moseley, 1877, see Winsor, Johns & Barker, 2004).
After receiving multiple reports from non-scientists on the presence of "large and colored" terrestrial flatworms in several localities in the IP, and given their observed locations, particularly in private gardens, we decided to analyze their presence in garden centers and plant nurseries.
The aims of this work were to: (1) estimate the number of terrestrial flatworm species introduced in the IP, and find their region of origin; (2) check whether plant nurseries and garden centers are acting as entrance gates and reservoirs; (3) estimate the invasive potential of some introduced species by considering their diet and by using Species Distribution Modeling (SDM); (4) propose measures to prevent their becoming invasive and to prevent further introductions and spread.

Specimen collection
Specimens were sampled from four sources (Tables 1 and 2): (1) gardens, (2) nurseries and plantations, (3) semi natural areas, and (4) from other countries (either the original area of distribution or other invaded areas). Specimens from sources 1 and 2 were either sent by people who knew our work through the information in social networks, or sampled by us (all the localities reported by non-scientist collaborators correspond to gardens). Specimens from source 3 were sampled by us. Specimens from source 4 were sent by colleagues, specialists of the group, to whom we requested material for comparison with the Iberian populations.
Data from a total of 13 domestic gardens, seven nurseries, two plantations (all confined, humanized locations), and three semi natural areas (humanized environments that are not confined and in direct contact with agricultural and forest areas) have been analyzed (Table 1). The three "semi natural areas", located in North-eastern Iberian Peninsula, were: (1) Cal Tet, Parc Natural Delta del Llobregat, Barcelona (Fig. 3, Loc-code O); (2) Can Cabanyes, Granollers, Barcelona (Fig. 3, Loc-code M); (3) Viaducte de Rubiò, Vall d'en Bas, Girona (Fig. 3, Loc-code P). In all three places recent habitat restoration activities have been performed, including the transplantation of autochthonous plant species from commercial nurseries.
Amateur collaborators photographed the animals alive and fixed them in absolute ethanol. Specimens we collected were also photographed and external morphological characters recorded. Subsequently, animals were subjected to two different procedures to proceed to the species identification: (1) specimens for molecular analyses were fixed in 100% ethanol and (2) specimens for histological studies were killed with boiling water, fixed with 10% formalin for 24 h, and then preserved in 70% ethanol.

Morphological studies
Preserved specimens were examined under a stereo microscope and notes of their dimensions, appearance, color (though this is affected by preservation), eyes, any stripes or pattern, the position of the pharyngeal aperture (mouth) and gonopore, if present, were taken. Specimens with no visible gonopore were considered to be immature. It was possible to identify some specimens, even immature ones, to species level without further examination. For unrecognized specimens, or where identity was uncertain and required confirmation, a mature specimen (evidenced by an open gonopore) was selected and divided into various portions, being embedded in wax. The copulatory apparatuś   Table 2).
(gonopore) and a small anterior region were sagittally and transversely sectioned at 10 or 15 µm, respectively, stained in Harris' haematoxylin and eosin and mounted in Canada balsam.

DNA extraction, gene amplification and sequencing
A small piece of tissue fixed in absolute alcohol was digested with Wizard Genomic DNA Purification lysis Buffer (Promega, Madison, WI, USA) and Proteinase K overnight at 37 • C, following manufacturer's instructions. The rest of the tissue is kept as voucher in the Genetics Department (Universitat de Barcelona). We amplified an approximately 1 kb fragment of the mitochondrial cytochrome c oxidase I (Cox1 gene) and a fragment of approximately 1,500 bp of the 28S rRNA gene (28S) by PCR reaction. PCRs were carried out in a volume reaction mixture of 25 µl. For Cox1 we used primers BarS (Álvarez-Presas et al., 2011) and COIR (Lázaro et al., 2009) and conditions were as inÁlvarez- Presas et al. (2011); 28S rDNA gene was amplified in two different overlapping fragments using the primers 28S1F, 28S4R, 28S2F and 28S6R, and conditions as inÁlvarez-Presas, Baguñà & Riutort (2008). Amplification products were purified with a vacuum manifold (Multiscreen HTS Vacuum Manifold; Millipore Corporation, Billerica, MA, USA). DNA sequences were determined from both strands using Big-Dye Terminator (3.1, Applied Biosystems, Foster City, CA, USA) and the reaction products were separated on the ABI Prism 3730 automated sequencer (Unitat de Genòmica dels Centres Científics i Tecnològics de la UB).
PCR products of the 28S gene for some individuals, that yielded double bands in the direct sequences, were cloned using HTP TOPO TA Cloning Kit for Sequencing (Invitrogen) in order to be sure that only one type of sequence was recovered (since the existence of a duplication of the ribosomal cluster is known in terrestrial planarians, Carranza et al., 1996). The sequences of the clones showed that these bands corresponded Table 2 Sequenced specimens. To each new sequence a three digit numeric code was assigned. Sequences from the GenBank database do not have specimen code numbers, only when there are more specimens from the same species in the same locality was a specimen code assigned (three letters + one number). Loc codes are as described in Table 1 Table 1.
to polymorphisms of one of the types. Seqman (v. 4.2.2, Gene Codes) was used to revise the chromatograms and obtain the definitive sequences.

Molecular data analyses
Ribosomal sequences were aligned using MAFFT v. 7 (Katoh & Standley, 2013) with the G-INS-i iterative refinement method and 1000 cycles. Mitochondrial coding DNA sequences were translated into aminoacids and aligned manually in Bioedit v.7.0.9.0. (Hall, 1999). All sequences were unambiguously aligned. We estimated the DNA sequence evolution model that best fits the data for both molecules using jModelTest 2.1.4. (Darriba et al., 2012), applying the Akaike information criterion (AIC). Phylogenetic relationships were estimated by Maximum Likelihood (ML) using RAxML 7.0.0 software (Stamatakis, 2006) and Bayesian inference (BI) using MrBayes v. 3.2. (Ronquist et al., 2012). Bootstrap support (BS) values were obtained for ML trees from 10,000 replicates. In the BI analyses we ran four chains to allow heating and used default priors, three million generations were run using the Markov Chain Monte Carlo (MCMC) analysis in two independent runs. Sampling was every 1,000 generations. The stationarity and convergence of the runs were checked by plotting Log likelihood values vs. number of generations and inspecting when the standard deviation of split frequencies had reached <0.01, respectively.

Potential distribution modeling
Using data describing the known distribution of C. coerulea in Australia, we estimated the potential distribution of this species in the Iberian Peninsula, as an exercise to find out whether climatic variables could detect potentially at risk areas where the establishment of the introduced species will be favored if only affected by climate. This could be a tool to help limit potential activities in order to avoid the introduced animals becoming invasive in the most likely areas for them to be successful.
For the SDM, a total of 179 Australian geographical coordinates of presence observations extracted from the literature, internet sources and personal communications (L Winsor, 2013) were used for calibration of models (training dataset). To avoid overparameterization and loss of predictive power, we discarded the climatic variables that were highly correlated. To do this we extracted environmental information from 10,000 randomly generated points and determined the linear relationships among them using Spearman and Pearson correlations. Although all correlations were significant they show low correlation coefficients (r ≤ 0.12). According to this analysis we used the 9 bioclimatic variables from the WorldClim database v. 1.4. (http://www.worldclim.org/, Hijmans et al., 2005) with less dependence, to form the present climatic dataset at a scale of 30 arc s. Those variables were: annual mean temperature; mean diurnal range; isothermality; maximum temperature of warmest month; minimum temperature of coldest month; precipitation of wettest month; precipitation seasonality; precipitation of wettest quarter; and precipitation of warmest quarter. The maximum entropy model, a presence-only algorithm that requires known species occurrence points and environmental variables (Maxent v.3.3.3k;Phillips, Anderson & Schapire, 2006), was applied. We selected the software default values for the convergence threshold, regularization values, and features. The maximum number of iterations was set to 1,000 and 1,000 bootstrap replicates were used. All possible geographic locations were partitioned between training and test samples (75% and 25%, respectively) in order to achieve higher predictive accuracy (Phillips & Dudík, 2008). Once the models were trained, we projected the results using the IP climatic dataset, to study the possible expansion of C. coerulea in the region. Model performance was evaluated using the AUC test (area under the receiver operating characteristic curve (ROC)) and the binomial test of the omission-dependent threshold was calculated by Maxent. Finally, binary maps of the outcome of the models were overlapped in the geographic information system, ArcMap v.10 (ESRI 2011. ArcGIS Desktop: Release 10. Redlands, CA: Environmental Systems Research Institute).

Morphological identification of the specimens
Based on the external appearance of the flatworms we initially grouped the specimens into nine morphotypes. We classified four of them at the species level due to their characteristic shapes or other external features, and the other five at genus or tribe level.
Bipalium kewense (Fig. 4) has been identified by the characteristic shape of the anterior end and the pattern of stripes along the dorsal and ventral body surfaces. One specimen preserved in 70% ethanol from Bordils locality (Loc code V in Table 1) has been deposited at the Natural History Museum of the United Kingdom (NHMUK) with voucher number NHMUK 2014.5.13.6.
For Caenoplana bicolor (Graff, 1899) there is no published description of a sexually mature specimen, hence the identification of the only specimen obtained, also an immature individual, relied exclusively on its external appearance (Fig. 5). This specimen is deposited in the tissue collection of the Genetics Department (Universitat de Barcelona). Among the specimens with an external morphology initially ascribable to the Caenoplana coerulea phenotype, we have found two morphotypes basing on their color pattern. Morphotype Ca1 (Fig. 6) presents a dorsal coloration in dark blue with a yellow middle-dorsal stripe, and a ventral light blue region (characteristic pattern of Caenoplana coerulea). The histological study of one specimen from El Prat de Llobregat locality (Loc code O in Table 1) (NHMUK 2014.5.13.14) reveals that it may belong to the Caenoplana coerulea species. Morphotype Ca2 (Fig. 7) presents a light brown dorsal region with a pale yellow middle-dorsal stripe, and a ventral light blue-greenish region. The histological study of one specimen from Bordils locality (Loc code V in Table 1) (NHMUK 2014.5.13.12) has revealed that its copulatory apparatus characters do not fit any of the described Caenoplana species.
Kontikia ventrolineata (Dendy, 1892) ( Fig. 9) was externally identified, following Great Britain Non-Native Species Secretariat (2013). We assigned the specific name following Jones, Johns & Winsor (1998), who considered Parakontikia Winsor, 1991 as a junior synonym of Kontikia Froehlich, 1955. Three specimens from Granollers locality (Loc code M in Table 1) are deposited at the NHMUK (NHMUK 2014.5.13.3-5). We found one morphotype externally ascribable to the genus Rhynchodemus, but not to a known species (Fig. 10). Rhynchodemus morph Rs1 has a dark brown pigmented body with two black longitudinal stripes, and two large eyes situated a little distant from the anterior tip. One specimen from Vall d'en Bas locality (Loc code P in Table 1) (NHMUK 2014.5.13.9) was histologically studied but, unluckily, presented a copulatory apparatus not well developed, preventing us from determining whether it could belong to Rhynchodemus sylvaticus (Leidy, 1851) to which it was extremely externally similar.
A morphotype externally ascribable to the tribe Rhynchodemini was found in Benamargosa locality (Loc-code G), but its morphological features did not allow assigning it to any genus. Rhynchodemini morph Ri1 presents a dark brown pigmented body with one dorsal black line (no image available).

Phylogenetic results
We inferred ML trees to check the diagnosis of the introduced specimens and to determine their level of relatedness to the ones from the original areas of distribution. For this reason, the datasets included, when possible, sequences belonging to morphologicallý diagnosed specimens from the original area of distribution of the putative introduced species (obtained for this study or coming from GenBank; Table 2).
We obtained 28S sequences for 15 individuals. One or two sequences from each morphotype were aligned together with 19 GenBank ingroup sequences and 3 outgroup sequences belonging to the Dugesia genus (terrestrial planarians sister group; Carranza et al., 1998;Álvarez-Presas, Baguñà & Riutort, 2008). Cox1 sequences were obtained for all individuals included in the study (Table 2). To obtain a more detailed picture of the situation within the main clades, including introduced planarians found on the concatenated analysis, we split the Cox1 sequences into four new datasets, one for each subfamily, tribe or genus: Caenoplanini (56 ingroup + 4 outgroup), Geoplaninae (26 ingroup + 2 outgroup), Bipaliinae (9 ingroup + 3 outgroup) and Rhynchodemini (29 ingroup + 4 outgroup). For each clade, its sister group was selected as the outgroup as shown on the concatenated analysis and/or previous studies (Álvarez-Presas, Baguñà & Riutort, 2008). The best-fit model of sequence evolution for the 28S was GTR + G and for Cox1 was GTR + I + G. We inferred a ML tree with partitions from a concatenated dataset including 37 individuals for which both 28S and Cox1 sequences had been obtained (Fig. 12). The ML trees obtained from the Cox1 datasets are shown in Figs. 13-16.
For the concatenated dataset, the ML tree showed most introduced specimens constitute monophyletic groups together with representatives of their species coming from the original distribution area or other introduced localities. We have found introduced planarians in the IP for all non-autochthonous terrestrial planarians subfamilies; in the case of the Rhynchodeminae there are even representatives from two tribes (Rhynchodemini and Caenoplanini).
Within the Bipaliinae, Bipalium specimens found in the IP constitute a monophyletic group together with Bipalium sequences from other species, B. adventitium being the closest relative in the Cox1 tree (Figs. 12 and 13). The genetic diversity among the four B. kewense sequences, coming from the IP and Açores Islands, was very small.
In the Geoplaninae clade (Figs. 12 and 14) the introduced specimens found in the IP constitute a monophyletic group with a still not-described species from Brazil (Obama sp. 6 after Carbayo et al., 2013, Fig. 14). In the Cox1 tree, specimens coming from the IP, United Kingdom (both introduced) and Brazil (original area) constitute a highly-supported monophyletic group. Within this group, the introduced individuals aré  Fig. 14), also distinctly separated from the Brazilian individuals. All the UK individuals fall within the clade Obama sp.A.
The Caenoplanini clade (Figs. 12 and 15) includes a high number of introduced individuals and the broadest diversity of sequences. Even Caenoplana coerulea sequences, either coming from GenBank, or from the individuals sent by our collaborator in Australia, are found in very distinct genetic clades pointing to the existence of more than one species (see Discussion). For this reason, we use the name Caenoplana coerulea s.l. to refer to all those specimens. In the concatenated tree, the representative of Caenoplana morphotype Ca1 is closely related to Caenoplana coerulea s.l. from Australia, while Caenoplana morphotype Ca2 is the sister group of a clade constituted by C. coerulea s.l. and C. bicolor. The divergence among these three lineages can be appreciated when compared to the otheŕ   subfamilies present in the tree. In the Cox1 tree ( Fig. 15) genus, Caenoplana again shows high levels of genetic diversity, evidenced by the long branches separating its subclades. Most Caenoplana morphotype Ca1 from the IP constitute a low diversity clade including C. coerulea s.l. from its original area (Australia) and also from UK and Menorca (also introduced). This clade is sister to another group including C. coerulea s.l. originally from different localities in Australia (Sunnucks et al., 2006); however, the differentiation among these two clades is extremely high. The other two Caenoplana morphotype Ca1 individuals, coming from Townsville (Australia), constitute a highly differentiated clade that also includes a GenBank sequence identified only to the genus level and one of the introduced individuals. Finally, there is a clade including only introduced animals, one of them identified as C. bicolor and the rest as morphotype Ca2. The genetic differentiation between the two lineages within this clade is nonetheless extremely high.
In the Rhyncodemini clade (Figs. 12 and 16) we find representatives of three genera in the IP. Dolichoplana striata sequences form a monophyletic clade in the Cox1 tree, including three introduced animals in the IP and one coming from Brazil. The individualś  assigned to Rhynchodemini morphotype Ri1 collected in Málaga (Spain, Loc code G in Table 1; Vila-Farré et al., 2011) cannot be assigned to any species, although they probably belong to Dolichoplana given the relationships they show in the Cox1 tree. The four K. ventrolineata specimens analyzed constitute a monophyletic group with a low variability, the French representative being the more divergent. The genus Rhynchodemus is represented by at least three species in the Cox1 tree. Rhynchodemus sylvaticus (considered an European autochthonous species), Rhynchodemus morphotype Rs1, and a clade including two individuals from Panamá that we had ascribed to the Rhynchodemini bý  their external appearance, and they appear likely to belong to the genus Rhynchodemus. It should be noted that the specific identification of all R. sylvaticus specimens found in the IP (Boix & Sala, 2001;Mateos et al., 2009;Vila-Farré et al., 2008;Vila-Farré et al., 2011) have been made based exclusively on external morphology (for this reason all these specimens have been considered Rhynchodemus cf. sylvaticus). Rhynchodemus cf. sylvaticuś clade, including representatives from Spain, Portugal, UK and France, together with a specimen identified in a previous study (Rhynchodemus cf. sylvaticus (Canyamars)) is a sister group of a clade constituted by Rhynchodemus morphotype Rs1 and one specimen of R. cf. sylvaticus (specimen 219). Figure 3 and Tables 1 and 2 show the sampling localities of the animals analyzed in this study. In all the plant nurseries, only one species of terrestrial planarian was found (Bipalium kewense, Rhynchodemini Ri1 or Obama sp.), except in Bordils where six species were found (Table 1, Loc-code V). The rest of the localities also contained a single species, with the exception of Treto (a garden, Loc-code U) with two species, and the two "semi natural areas" situated in Vall d'en Bas (Loc-code P) and in Granollers (Loc-code M) also with two species each. Obama sp. was the species most frequently found in plant nurseries, while B. kewense predominated in private gardens. In the semi natural areas only the species K. ventrolineata, C. coerulea s.l., and Rhynchodemus Rs1 (not found anywhere else) have been found.

Potential species distribution modelling
The result of projecting models for the potential distribution of C. coerulea s.l. in the IP presents mean values of AUC beyond 0.9 (0.974) and significance for all tests of omission, which indicates good performance of the models. Furthermore, predictions were significantly different from random because binomial omission test thresholds were significant (p < 0.01) in all 1,000 runs. A composite map showing the potentiaĺ Figure 17 Potential distribution of Caenoplana coerulea species across the Iberian Peninsula. The color gradient indicates the predicted likelihood that the environmental conditions suitable for the species based on the MaxEnt average output. Letters indicate localities where C. coerulea has been found, locality codes correspond to those in Table 1. distribution models for C. coerulea s.l. species projected on current climate layers is provided in Fig. 17.
The results of the potential distribution of the species in the IP, based on data from its current distribution in their region of origin (Australia), show that the species can find extremely suitable areas for its survival and expansion is the northern region, where the appropriate temperature and humidity conditions occur.

Species identification, or, how many species are out there?
External morphology , analysis of histological sections, and phylogenetic inference from molecular data  have revealed the presence of five clearly identifiable species of introduced exotic land planarians in the IP: Bipalium kewense (Bipaliinae), Caenoplana bicolor, Caenoplana coerulea s.l. (Ca1), Dolichoplana striata (Rhynchodeminae, Rhynchodemini), and Kontikia ventrolineata (Rhynchodeminae, Caenoplanini). However, the phylogenetic trees obtained and the analysis of the external appearance of the specimens indicate that probably at least five more species were present, including Rhynchodemini morph Ri1, Rhynchodemus morph Rs1, Obama sp. and two more species within Caenoplana:Caenoplana morph Ca2 and probably some individuals of Caenoplana morph Ca1 (see below).
The assignation of Bipalium kewense is based on its characteristic external morphology (see Hill & Merickel, 2011;Jones, 1998). There are no published Cox1 gene sequences for this species in Genbank, so those presented in this paper are the first available. Phylogenetic analysis of these sequences point to an introduction from the same lineage. Surprisingly, all sequences belonging to Kontikia ventrolineata (coming from Spain, France and UK) are situated within the Rhynchodemini clade with high support in both trees. This situation contradicts the taxonomy proposed by Sluys et al. (2009) where the genus Kontikia belongs to tribe Caenoplanini.
The genetic differentiation observed within the group constituted by the genus Caenoplana, monophyletic in the trees, leads us to predict that it includes more than one species. In the Cox1 tree (Fig. 15), at least three monophyletic groups seem to be clearly defined and probably represent different species. In fact, C. coerulea is considered by a specialist in this group (L Winsor, pers. comm., 2013) as a complex of species, on the basis of internal anatomical characters and stripe morphology. According to Winsor, there are at least three species that are distinguishable morphologically; but there are probably more than three species in the area of origin. One of the problems with the group is that the type of the species is non-sexually mature. Hence, to clarify the situation and number of the species in this group, a broad sampling in its original area of distribution is required, followed by a thorough morphological and molecular study. Nonetheless, for the purpose of the present paper, the evidence is clear that at least three different genetic lineages from Australia have been introduced in the IP, probably independently.
In the case of Rhynchodemus Rs1, we cannot be sure if this is a distinct species or simply a differentiated lineage of R. cf. sylvaticus. The latter has been generally regarded as an introduced species in Europe from USA (Jones, 1988), but it is also considered as probable species native to Europe (Jones, 1998;Jones, 2005) and introduced in the USA from Europe (Ogren & Kawakatsu, 1998). The type locality of R. sylvaticus is Philadelphia, Pennsylvania, USA (1851). This species has a wide distribution in the IP and two of the locations are plant nurseries, one in Barcelona (Vila-Farré et al., 2008) and one in Málaga (Vila-Farré et al., 2011), while the other localities can be considered natural habitats. In our molecular analysis there was no separation of specimens according to their locality type (natural or artificial). Two distinct clades of European Rhynchodemus were obtained (Fig. 16), suggesting the existence of two different species with similar external morphology.
In the case of Rhyncodemini Ri1, this species probably belongs to the genus Dolichoplana; however, we were only able to obtain three specimens and none of them were sexually mature.
When specimens of Obama sp. were first found in the UK and the IP, they were identified as O. marmorata (Schultze & Müller, 1857) due to their external appearance; however, molecular data (M Riutort, unpublished data, 2014) showed that the European specimens did not constitute a monophyletic group with that species, indicating that they belonged to an unknown, still undescribed, Geoplaninae. Sampling performed in Brazil since then has found another species (Obama sp.6), which is also externally very similar. Molecular data show that it is closely related to the individuals found in Europe (M Riutort, unpublished data, 2014). As in the previous case, a morphological and molecular study should be undertaken to clearly delimit and describe the new species. The two clades found in our Cox1 tree (Fig. 14), that may represent two different species, suggest that there have been two independent introductions into the IP from different native sites in Brazil.
Overall, we have shown that at least ten introductions have occurred in the IP. These introductions include species from all the non-European terrestrial planarian subfamilies from native localities as far as South America and Australia. Since most of these species have previously been reported to have been introduced in other countries, the introductions into the IP have probably not been directly from the source countries, but were more likely to be indirect, following plant trade routes. In most cases, all the individuals from the same species found in the different localities are nearly identical, even when compared between Spain and the UK, which can be interpreted as the result of a single introduction (or a single exportation from the place of origin). In others, as in the case of Caenoplana, the observed diversity clearly indicates that the introductions were from different lineages within this group and is likely to be the consequence of more than one export from the native area.

What makes terrestrial planarians so successful as introduced species?
Temperature, humidity and food availability are the three basic factors determining the geographical distribution of terrestrial planarians (Boag, Yeates & Johns, 1998). The feeding habits of the introduced species in the IP indicate that all of them feed on invertebrate soil fauna (Table 3). In plant nurseries and greenhouses microclimatic conditions are maintained artificially (high humidity and stable temperature) and are likely to favor the presence of stable populations of many species of terrestrial invertebrates. In nurseries we visited, especially under flowerpots, we have observed the presence of numerous specimens of snails, slugs, earthworms, millipedes, isopods, beetles and various groups of microarthropods, including springtails. Therefore, in this very suitable artificial microhabitat, there is likely to be a greater number of species of terrestrial flatworms (as is the case of Bordils, Loc code V in Table 1, where six species were detected in the same greenhouse).
Land planarians and their cocoons are very often associated with the soil of plants in pots and certain types of fresh vegetables (Ogren, 1985;Mather & Christensen, 1992;Hogan & Dunne, 1996). The transport of these pots and materials (which can occur over international and intercontinental distances) may permit the transport of associated planarians and/or cocoons, which is the primary means of introduction of exotic species of terrestrial planarians into different contaminated countries (Winsor, Johns & Barker, 2004). The suitable conditions in the plant nurseries and garden centers may explain their introduction success. In recent decades, the adoption of free market policies and trade agreements have reduced barriers to plant trade among different countries (Dehnen-Schmutz et al., 2010), but there has been insufficient attention given to how such structural change in international trade can affect the risk of spread of invasive species (Drew, Anderson & Andow, 2010). Depending on the intricate network of commercial interactions among European countries (see Dehnen-Schmutz et al., 2010), we expect a huge European dispersal of exotic animal species associated with this trade.

Will planarians become invasive in the Iberian Peninsula as has occurred in other areas?
Exotic species present in an area could be categorized as introduced (detected in the area but with unknown status), adventives or not established (they reproduce occasionally in the area not constituting stable populations), naturalized or established (they form stable reproductive populations in the area) and invasive (established and well spread in the area) (Richardson et al., 2000;Simberloff et al., 2013). The "tens rule" (Williamson & Fitter, 1996;Williamson & Brown, 1986;Williamson, 1996) predicts that just one of hundreds of introduced species becomes invasive (about 10% of the introduced species are established, and that 10% of those become invasive). Based on the premise of the "tens rule", some researchers minimize the potential impact of exotic species (National Research Council, 2002;Campbell & Gibson, 2001), while others warn that this risk minimization is dangerous and, with respect to the possible impact of introduced species, the adoption of the precautionary principle is crucial (Jarić & Cvijanović, 2012), but unlikely! The problem with this sort of assumption or calculation is that, in most cases, we simply have no knowledge of the unsuccessful introductions.
In the case of terrestrial planarians, some species are very tolerant of habitat modification (Cannon et al., 1999;Carbayo, Leal-Zanchet & Vieira, 2002), facilitating theiŕ survival in humanized environments. Many introduced species of terrestrial planarians are found confined to these types of habitats (parks, private gardens, plant nurseries), but it is not known whether this distribution is so restricted due to environmental constraints (planarians, coming from tropical habitats cannot live outside these artificial habitats in the European environment) or to a low velocity of dispersion to natural habitats (Ducey & Noce, 1998). In our case, most specimens occurred in confined areas (gardens and nurseries). However, Rhynchodemus Rs1, C. coerulea s.l. and K. ventrolineata have been also found in recently restored areas that were more or less connected to natural and agricultural environments, which increases the danger of their becoming naturalized or even invasive.
In the particular case of C. coerulea s.l., we performed a potential distribution study to check whether the area around its present introduced localities in the IP may be suitable for its expansion. The results show that the potential distribution of the species (Fig. 17) indeed includes the countryside that was nearby to the localities of the IP where it is already present. The most suitable area is the northern IP. This is not surprising when we consider that in this northern region, the climatic conditions (temperature and humidity) are also more optimal for the presence of native land planarians (Mateos et al., 2009;Alvarez-Presas et al., 2012). Thus, we show that by having suitable climatic databases, it is possible to model the potential distribution of introduced species, and thus predict their risk of becoming invasive. If we add to this information the knowledge of some biological features of the terrestrial planarian species, such as their prey preferences, we may be able to make an even more precise image of the sites where it is more likely for the species to become invasive and thus concentrate prevention efforts in those areas.
Our results show that C. coerulea s.l. is apparently the most successful colonizer, since it is the only species present in all three unconfined (semi natural) areas sampled. This may be because it feeds on several groups of arthropods that are abundant in areas where this species has been detected (isopods, beetles, diplopods). The three species (Rhynchodemus Rs1, C. coerulea s.l. and K. ventrolineata) we find in unconfined environments feed on arthropods, whereas the other species (found only in confined environments) do not feed on arthropods, but instead on other invertebrates that require extremely wet habitats. Hence, land planarian species that feed on arthropods have their food "secured" in environments with a Mediterranean climate and, as a consequence, have a higher likelihood of being successful and becoming established or even invasive.

What consequences might the introduction of flatworms have on human economies and biodiversity?
Another important question is: what are the negative effects of the spread of these species? In literature, the primary problems reported are related to economic consequences for agricultural activities (Boag & Neilson, 2006;Boag, Neilson & Jones, 2010). As predators of earthworms, planarians can cause soil drainage and fertility to be severely compromised. The ecological consequences of the presence of these predators depends on theiŕ propagation speed and efficiency, but could have significant effects on processes mediated by earthworms in both agroecosystems and forests (Lilleskov et al., 2010). Although there is still no direct impact study of the presence of invasive planarians on agricultural production (Boag, Neilson & Jones, 2010), data from farmers with infested farmland and from the scientific literature have suggested that it could reduce grass yields significantly (Boag & Neilson, 2006).
No reference has been made to the effect of these species on autochthonous populations of terrestrial planarians, probably because the knowledge of the autochthonous fauna is very scarce. In the IP we have already performed some studies on the autochthonous terrestrial planarian fauna and found that it is very diverse, including at least 15 species, of which some contain a great deal of genetic diversity (Mateos et al., 2009;Álvarez-Presas et al., 2012). The potential arrival of some of these introduced species in natural habitats, where the autochthonous ones are localized (as predicted by the potential distribution studies), would have very negative consequences. Since exotic planarians are, in general, bigger in size, more voracious, have more aggressive behavior, and sometimes appear to have a generalist diet (pers. obs.), they may be more resistant to extreme conditions than the native species.

A cautionary tale: plant trading and landscape restoration
An important question raised by all these observations is whether governments in Europe should be asked to propose new, more restrictive rules on the trade of plants coming from outside, or alternatively, to establish better controls or protocols to avoid the introduction of unwanted organisms together with the plants. However, it is probably now too late to have an impact on the transport of species around the world. Nonetheless, we are still in time to avoid invasions of terrestrial planarians. The restoration of degraded areas involves planting native plant species. These plants are available from nurseries and transported to the restoration areas accompanied by a certain amount of soil on the roots. If this land is not subject to any preventive treatment, it may be contaminated with organisms that are also introduced in the area that is being restored. Among these organisms may be unwanted species that, if given the right conditions, can become invasive. It is important to warn agencies conducting such restorations of these dangers and ask stakeholders to include in the protocols of landscape restoration the necessary steps to avoid these unwanted introductions.
Some simple, easy-to-perform sanitizing procedures, such as heating the soil (EPPO, 2000a;EPPO, 2000b;SEERAD, 2000;Sugiura, 2008) before transplanting the nursery plants to the natural environment, may be sufficiently effective and reliable to ensure that there is no concomitant dispersal of flatworms. Such procedures, together with a periodic analysis of the introduced species present in garden centers and nurseries, and a study of the potential areas of flatworm distribution, would also help avoid the introduction of terrestrial planarians into areas where they are more likely to become invasive (DEFRA, 2005;DOVE, 2012).