Introduction

The order Tricladida represents one of the best-known groups of free-living flatworms, with representatives in all biogeographical regions of the world (Schockaert et al., 2008), of which the suborder Continenticola (including both terrestrial and freshwater representatives) houses the largest number of species. Recent studies have reported the introduction of many species of tropical terrestrial planarians all over the world (Álvarez-Presas et al., 2014; Justine et al., 2014; Jones & Sluys, 2016; Mazza et al., 2016; Mateos et al., 2020; Negrete et al., 2020), with some species representing a problem for agriculture due to their annelid-based diet depleting the soils from this important aeration fauna (Murchie & Gordon, 2013). These studies not only have unravelled the country of origin of these introduced species but have pointed out also the putative causes of these introductions, as well as potential areas for future range extension (Álvarez-Presas et al., 2014; Negrete et al., 2020), which represents vital knowledge for the development of control measures. In contrast, introductions of invasive freshwater planarians have received much less attention, of which the genus Girardia Ball, 1974 is of special importance.

The genus Girardia includes about 60 valid species, the natural distribution of which extends from the southern regions of Argentina to Southern Canada. Furthermore, species of Girardia have been introduced into many other regions of the world. The species G. tigrina (Girard, 1850) was first reported in Germany in the 1920s, after which it was recorded also from other parts of mainland Western Europe, as well as the Balearic Islands and the Azores, its spread most likely effectuated through the international trade in aquatic plants and the activity of aquarists (Stocchino et al., 2019 and references therein). Representatives of the genus have been reported also from Australia, Japan, and Hawaii (Kawakatsu et al., 1984; Sluys et al., 1995a, b; Sluys et al., 2005, 2010). In addition, records of Girardia in China, Serbia, and Eastern Europe have been published in the past few years (Chen et al., 2015; Ilić et al., 2018; Kanana & Riutort, 2019), and recently the genus has been detected for the first time in Africa (Mabrouki et al., 2023; Benítez-Álvarez et al., 2022).

For a long time, it has been assumed that only G. tigrina was introduced in Europe, although some studies have pointed to the possibility that other species belonging to the genus Girardia may be present, in view of the diversity in morphology, karyology, and reproductive biology of the populations analysed from the Iberian Peninsula (Ribas et al., 1989) and Italy (Benazzi, 1993; Stocchino et al., 2019). However, due to the fact that, generally, the introduced populations of Girardia show only asexual reproduction and do not possess a copulatory apparatus, their identification to species level is greatly hampered, since the main diagnostic characters reside in the reproductive structures (Sluys & Riutort, 2018). This problem of identification of asexual strains is solved using molecular markers, a procedure successfully applied in various studies on planarians, both for identification of species and for resolution of their phylogenetic relationships (see Álvarez-Presas & Riutort, 2014). Based on this methodology, a recent molecular study discovered that three species of Girardia belonging to North American lineages have dispersed to Europe and other regions of the world (Benítez-Álvarez et al., 2023a).

An aspect of considerable importance and concern regarding introduced alien species involves their potential invasiveness and detrimental effects on autochthonous ecosystems, an issue that is particularly important for freshwater habitats, since these communities are highly sensitive to compositional changes due to the introduction of invasive species (Havel et al., 2015). Unfortunately, only few studies have tried to evaluate the effects of introduced Girardia on local ecosystems, mainly in the United Kingdom (Pickavance, 1971; Velde, 1975; Wright, 1987; Gee & Young, 1993; Gee et al., 1998). Therefore, at present, we lack adequate knowledge on the potential danger of introduced Girardia to autochthonous species and possible cascading negative effects on the freshwater environment.

Ecological niche modelling offers the possibility to evaluate the potential distribution of species and to estimate, at the same time, the environmental parameters underlying the distributional spread. By analysing and comparing the environmental requirements of introduced, as well as native species, an indirect evaluation is provided of the potential influence exerted by the introduced species on autochthonous ones (Warren & Seifert, 2011; Radosavljevic & Anderson, 2014; Puchałka et al., 2020).

The main objectives of the present study were (a) determining the origin and current diversity of introduced European and worldwide populations of Girardia, (b) identifying their potential area of distribution and environmental requirements compared to some native species, and (c) predicting future tendencies in their potential distribution at European and global scales. To this end, we constructed a molecular phylogeny of the genus Girardia, including both native and introduced representatives from all over the world. In addition, a distribution modelling analysis was performed on the introduced species, to determine their potential area of distribution, with a special focus on the Iberian Peninsula, where the distribution of autochthonous species is well established and equivalent studies have already been performed. Finally, we ran a niche modelling analysis under two different climatic scenarios to predict future potentially suitable areas for the three introduced species.

Materials and methods

Data acquisition

Sequences of Cytochrome Oxidase I (COI) and Elongation factor 1 alpha (EF1α) genes from Girardia representatives from Asia, Australia, Hawaii, the Americas, and Europe were downloaded from GenBank. We analysed a total of 98 Girardia representatives, originating from 76 localities that are distributed all over the native range of the genus as well as from several regions in which species were introduced; Europe constitutes the most densely represented region (Fig. 1; Online Resources Table S1).

Fig. 1
figure 1

Maps showing the sampling localities for all individuals analysed in the present study. A The Americas, B Western Europe, C Hawaii, and D Asia and Oceania. For a finer resolution, see the interactive map at https://www.ub.edu/planarian-maps/

Full size image

Sequence alignment

Sequences of COI and EF1α were aligned independently with ClustalW on the BioEdit Sequence Alignment Editor using default parameters (Hall, 1999). Each gene was translated into amino acids with the corresponding genetic code to check for the absence of stop codons and to produce the alignment and, thereafter, converted again to nucleotides. A concatenated dataset of both genes was obtained in Mesquite v3.04 (Maddison & Maddison, 2015). Missing data, due to the short length of some downloaded sequences or the lack of one of the two genes for some individuals, were coded by Ns. PartitionFinder v2.1.1 (Lanfear et al., 2012) was used to find the best partition scheme for each gene, considering the score for the Bayesian Information Criterion (BIC). The results of the PartitionFinder program validated a partition scheme considering each codon site independently, both for COI and EF1α. For each partition the best model was General Time Reversible + Gamma Distribution + Invariable Sites (GTR + Г + I).

Phylogenetic inference

Bayesian Inference (BI) method was used to infer a phylogenetic tree from the concatenated dataset. In a previous phylogenetic study of the genus Girardia (Benítez-Álvarez et al., 2023a) based on the same genes, maximum likelihood resulted in identical topologies to Bayesian inference, albeit with lower supports. However, we found that bootstrap supports were sensitive to the high proportion of missing data in some sequences, due to the resampling strategy of the method, while it has been shown that posterior probability increases sensitivity to the phylogenetic signal, as compared with bootstrap values (Alfaro et al., 2003). Additionally, empirical and simulated data point to the advantage of including missing data over the radical elimination of genes or taxa (Wiens, 2003, 2006; Fulton & Strobeck, 2006; Jiang et al., 2014). Furthermore, other studies have also pointed out the better performance of the BI method when analysing alignments with high missing data content (Wiens & Morrill, 2011 and references therein). For these reasons, we decided to perform exclusively a Bayesian analysis in the present study. Bayesian analyses were run in MrBayes v3.2.2 (Ronquist et al., 2012) with 10 million generations, sampling every 1,000 generations, and 25% burn-in (default setting) to obtain the tree with the highest probability and the posterior probability values. The best evolutionary model and partition scheme obtained in PartitionFinder was applied with parameters unlinked. Convergence of the parameter values and topologies was examined by checking that the average standard deviation of split frequencies was below 0.01. Estimated sample size values (ESS) of each run were inspected in Tracer v1.5 (Rambaut et al., 2018) and only runs with values for all parameters above 200 were accepted.

Species distribution modelling

Based on the molecular analysis, we identified the lineages of Girardia that had been introduced into different regions of the world, corresponding to three species, viz., G. tigrina, G. sinensis Chen & Wang, 2015, and G. dorotocephala (Woodworth, 1897) (see Results). The software MaxEnt 3.4.0 (Phillips & Dudík, 2008; Phillips et al., 2017) was used to infer the potential areas of distribution of these three species, both at a global scale and for the Iberian Peninsula in particular. In this analysis, worldwide distributional presence data was included for each of the three introduced Girardia species, including localities in both native and non-native regions. Six hydro-environmental variables were downloaded from the RiverATLAS database (Linke et al., 2019), at a spatial resolution of either 15 or 30 arc-seconds: Temperature, Natural Discharge, Precipitation, Land Cover Classes, Lithological Classes, and Slope. These variables were specifically chosen because they had already been shown to be informative (and non-correlated to each other) to predict the potential distribution of freshwater planarians (Leria et al., 2022). Future potential distribution of these three Girardia species was inferred under two different scenarios of climate change. The bioclimatic variables Annual Mean Temperature (Bio1) and Annual Precipitation (Bio12) were downloaded from the WorldClim2 database (Fick & Hijmans, 2017) at a spatial resolution of 2.5 arc-min for the time periods corresponding to Current, 2021–2040, and 2081–2100. For each of the future time periods, two scenarios of greenhouse gas emissions were considered, viz., (1) highest decrease of emissions (SSP1–2.6), and (2) highest increase of emissions (SSP5–8.5) (see https://www.carbonbrief.org/cmip6-the-next-generation-of-climate-models-explained for the detailed characteristics of each scenario; last date visited 20 September 2022). All MaxEnt analyses were performed by using 10 cross-validated replicates and setting the output to cloglog format. In each cross-validated replicate all data was used for validation, thus making a better use of small datasets than setting a single training/test split. MaxEnt was allowed to automatically determine 10,000 background points as pseudo-absences and to calculate and apply the regularisation coefficients for each variable. For each analysis, performance of the model was evaluated by checking the Area Under the Curve (AUC) of the Receiver Operating Characteristic (ROC) (Fielding & Bell, 1997). In total, the following number of localities could be used for each species: 16 (G. tigrina), 17 (G. sinensis), 18 (G. dorotocephala) (Table S1). We are aware that these values are well suited for a narrow-ranged species, but may be somewhat scarce for a widespread species, for which 20 to 30 records are recommendable (van Proosdij et al., 2016). As a result, our AUC values are expected to be slightly inflated. However, as we included all presently available data, we considered it important and worthwhile to perform the niche modelling analysis to make a preliminary prediction of the potential ranges of distribution of the three invasive species of Girardia.

Results

Sequences, alignment, phylogenetic analysis, and distribution

A total of 98 Girardia specimens were used for the present work. The concatenated dataset included 94 sequences of COI (837 base pairs, bp) and 78 of EF1α (879 bp) (Table S1). To avoid rooting with outgroup taxa that might be too distantly related to the ingroup, the tree was rooted on G. schubarti (Marcus, 1946) + Girardia sp. 1 (below indicated as the G. schubarti group). A recent study on the phylogenetic relationships within the genus Girardia supports this rooting (Benítez-Álvarez et al.,  2023a).

Our phylogenetic tree (Fig. 2) recovered monophyletic groups for the individuals of each of the identified species and also delimited seven clades for the unidentified individuals originating from Mexico (Girardia sp. 2 and sp.7), Brazil, and Chile (Girardia sp. 1, and Girardia sp. 3 to sp. 6) all included in a South American clade, collapsed in the present work for the sake of clarity (Fig. 2; Table S1; see Benítez-Álvarez et al., 2023a for a further discussion on these groups). All of these clades were fully supported (> 0.99 PP), except for G. sanchezi (Hyman, 1959) (0.67 PP).

Fig. 2
figure 2

Bayesian Inference tree inferred from the concatenated dataset (COI + EF1α). Most clades have been collapsed for the sake of clarity, excepting those of the three invasive species of Girardia. The text boxes at the right column indicate for certain clades the name of the species included and their areas of distribution. Circles at nodes indicate posterior probability. Scale bar number of substitutions per nucleotide position

Full size image

In the phylogenetic tree (Fig. 2), all non-American specimens grouped into three clades, together with American individuals corresponding to the species G. tigrina, G. sinensis, and G. dorotocephala. The topology shows that these three species are closely related, being connected by relatively short branches. The former two species form a clade, which shares a sister–group relationship with the last-mentioned species. Together, the clade of these three species is sister to an unidentified species from Mexico (Girardia sp. 7, 1.00 PP).

The G. tigrina clade (1.00 PP) comprises individuals from Michigan, USA (code 679sp_USADL), and Nova Scotia, Canada (code 659sp_CanNS). All other specimens in this clade come from various European countries, viz., France, Spain, Germany, and Italy (Fig. 2; Online Resources S1). It is noteworthy that the G. dorotocephala clade (1.00 PP) comprises individuals from widely separated areas, such as USA, Mexico, Japan, Hawaii, Brazil, Portugal, Spain, and France (Fig. 2 and S1). In similar vein, it is noteworthy that the G. sinensis clade (1.00 PP) comprises individuals from Australia, China, Cuba, Spain (mainland, as well as the island of Menorca), Italy (Sardinia Island), France, Germany, and The Netherlands (Fig. 2 and S1).

Niche modelling

The distribution models inferred with MaxEnt for G. tigrina, G. sinensis and G. dorotocephala gave AUC values close to 1.00 (Table 1), indicating that all models had high predictive power for the potential distribution of these species. Two of the environmental variables that contributed the most to the current models were Temperature (T) and Land Cover Classes (LCC), which in all cases showed a combined percentage of contribution > 60%, while the variables that contributed less were Slope (S) and Natural Discharge (ND) of rivers (Table 1).

Table 1 Contribution, in percentages, of the various variables used in the niche modelling for predicting the potential distribution of Girardia tigrina, G. sinensis, and G. dorotocephala
Full size table

In the following we have considered as being of “low suitability” values equal or below 0.5, “moderately suitable” values above 0.5 but below 0.7, and “highly suitable” being represented by values equal to or above 0.7. Regarding LCC, all three clades showed very high suitability for artificial surfaces and associated areas (e.g., canals and dams), and high suitability for different types of tree cover (with differences in their preferences, depending on the species), and for cultivated and managed areas (Fig. 3). Furthermore, G. dorotocephala also showed high suitability for mosaic regions of croplands and tree cover (Fig. 3). The most suitable Temperature (T) for G. sinensis and G. dorotocephala was ~ 15°C, while for G. tigrina it was slightly lower (~ 11°C), in all cases with a broad range of moderately suitable temperatures (Fig. 3). Regarding the Lithological Classes (LC), the only terrain that was determined as being suitable for all three species was the one composed of carbonated sedimentary rocks, while the rest of lithological preferences were different for each species (or shared between two of the three species) (Fig. 3).

Fig. 3
figure 3

Relative effects of the six hydro-environmental variables used in the niche modelling for predicting their effect on the potential distribution of Girardia tigrina, G. sinensis, and G. dorotocephala. Y-axis: habitat suitability ranging from 0 to 1 (low suitability < 0.5, moderate suitability ≥ 0.5 < 0.7, high suitability ≥ 0.7); X-axis: for continuous variables (Temperature, Precipitation, Slope, and Water discharge), presenting range of variation of each variable. The graph line for the continuous variables represents the mean effect of different replicates. Land Cover and Lithological Classes include only classes with high suitability

Full size image

All three clades showed the highest suitability for regions with moderate precipitation regimes with an annual average of approximately 1,000 mm (1,000 l/m2). However, the three species show moderate suitability for a wide range of rainfall regimes, going from scarcely any rain during the year (close to 0 mm/year) to regimes reaching accumulated precipitations of 3,000 mm/year (G. dorotocephala) and 4,000 mm/year (G. sinensis) (Fig. 3), with G. tigrina showing the steepest decay in suitability over 1,000 mm/year. No clear preference was detected for Water Discharge (WD) of rivers in any of the three species. Girardia sinensis and G. dorotocephala were found to be highly suited to regions with a moderate terrain slope (between 0° and 15°), thus differing from G. tigrina, which turned out to be highly suited to a wide range of terrain slopes (between 5° and 40°).

Predicted current potential distribution of the three species revealed a predominant suitability for the Northern Hemisphere, especially for regions such as Eastern USA, Europe, and Eastern Asia (Fig. S1). Additionally, these three species also showed some potentially suitable regions in the southern parts of South America, Africa, and Australia, with G. sinensis presenting the broader continuous potential distribution in these southern areas.

On the Iberian Peninsula, G. tigrina and G. sinensis showed extensive regions of high suitability, with the former having its most extensive area of high suitability being located in the northern region of the Peninsula, particularly in the area ranging from the Basque Country to Asturias (Fig. 4). In contrast, maximally suitable conditions for G. sinensis were found in the main rivers in the western region of the Iberian Peninsula, especially in the hydrological basin of the River Tajo. However, it should be borne in mind that the rest of regions of the Iberian Peninsula showed a high suitability for these two species. In contrast to G. tigrina and G. sinensis, MaxEnt did not predict any extensive regions of high suitability in the Iberian Peninsula for G. dorotocephala. Nonetheless, there were some areas, such as the surroundings of Lisbon or some scattered spots in Catalonia, Andalusia, and the northern part of the Iberian Peninsula, where potentially highly suitable areas were predicted (between 0.8 and 1.0) for this species.

Fig. 4
figure 4

Potential geographic distribution of Girardia tigrina, G. sinensis, and G. dorotocephala in the Iberian Peninsula, estimated through niche modelling by using six hydro-environmental variables. Colour scale from blue to red indicates habitat suitability, ranging from 0 to 1 (low suitability < 0.5, moderate suitability ≥ 0.5 < 0.7, high suitability ≥ 0.7). Black pixels: no hydro-environmental data

Full size image

The predicted future potential distribution of the three introduced species revealed a gradual northwards expansion into suitable areas in the Northern Hemisphere and a gradual reduction of their suitable areas in the Southern Hemisphere (Online Resources Fig. S2). Although this pattern was observed in both of our two scenarios on greenhouse gas emissions, it was much more extreme in the scenario with the highest increase in emissions (SSP5–8.5) than in the scenario with the highest decrease in emissions (SSP1–2.6). Figure 5 represents a summary of these results focused on our main area of interest, Europe and the Mediterranean region. In this area, G. tigrina, currently having the more extensive potentially suitable area of the three species, will be able to extend its distribution to Scandinavia in the best scenario (SSP1–2.6), reaching the northern coast of Norway in the worst-case scenario (SSP5–8.5), while it will have no suitable areas in the North of Africa. For the other two species, the variations are less pronounced in the best scenario (SSP1–2.6), while in the worst-case scenario they will be able to expand also over nearly all of Europe (except high mountains, such as the Alps) and reduce, but do not completely lose, their suitable regions in Northern Africa.

Fig. 5
figure 5

Predicted future geographic distribution of Girardia tigrina, G. sinensis, and G. dorotocephala in Europe, estimated through niche modelling under two different scenarios of climate change (SSP1–2.6: low greenhouse gas emissions; SSP5–8.5: high greenhouse gas emissions) by using the bioclimatic variables Mean Annual Temperature and Annual Precipitation. Colours correspond to three different time periods (green: Current; blue: 2021–2040; red: 2081–2100). Dotted lines denote the boundaries of the overlapping distributions. Detailed figures for each species and scenario are found in online resources Fig. S2

Full size image

Discussion

Three North American species of Girardia colonise the world

Our phylogenetic analysis showed that all non-American Girardia samples belonged to either G. tigrina, G. sinensis or G. dorotocephala. The close phylogenetic relationship between G. tigrina and G. dorotocephala and the fact that they are of North American origin, suggest that they diversified only relatively recently on the North American continent. However, G. sinensis, which is phylogenetically closely related to G. tigrina and G. dorotocephala (Fig. 2), was described on the basis of specimens from China, but the phylogenetic tree shows the presence of this species also on the Island of Cuba, which lies on the North American Plate (Fig. 2; Fig. S1). Therefore, and on the basis of the phylogenetic relationships and the geographic distribution of these three species, G. sinensis is actually a North American species that was introduced into China, but thus far is undetected in North America (Benítez-Álvarez et al., 2023a).

Thus, our results corroborate (1) the earlier hypothesis that G. tigrina was introduced from North America, and (2) the introduction of G. dorotocephala in Japan and Hawaii. In addition, our results revealed the introduction of G. dorotocephala in Europe and Brazil and point to a third introduced species, G. sinensis, that has been introduced from North America into China, Europe, and Australia (Figs. 1, 2).

Our present results indicate multiple introductions of Girardia into Europe, corroborating previous results of Ribas et al. (1989) and Stocchino et al. (2019), who had already described the presence of four morphological biotypes (classes A–D) of G. tigrina. Ribas et al. (1989) detected genetic differentiation between animals of class A and those of classes B + C and suggested that this may signal the presence of at least two separate taxa, either at the subspecies or species level. Moreover, both Ribas et al. (1989) and Stocchino et al. (2019) suggested that presence of the four morphological groups in Europe may be the result of multiple independent introductions. Our present results corroborate this notion of multiple introductions, albeit that now we can establish molecularly that two species are involved in the areas where Ribas et al. (1989) and Stocchino et al. (2019) analysed the presence of G. tigrina, viz., G. tigrina and G. sinensis.

With respect to their external appearance and internal morphology, G. sinensis and G. tigrina are highly similar. This paucity of morphological differences between both species may be the reason why G. sinensis has gone undetected in its original area of distribution, North America. In this sense, it is noteworthy that the copulatory apparatus of the sexually reproducing animals on Menorca analysed by Ribas et al. (1989) conforms to that of G. tigrina, and their habitus and pharynx pigmentation pattern to class C, i.e., spotted on a greyish background with two longitudinal stripes delimiting a lighter midline area and a pharynx very lightly pigmented (Ribas et al., 1989, Fig. 2). On the other hand, the specimens from Menorca analysed in our study, and which came from a locality close to the B-type locality of Ribas et al. (1989) inhabited by fissiparous individuals with a pigmentation pattern similar to that of class C, belong to G. sinensis. All of this provides substantial evidence that G. tigrina and G. sinensis have colonised the same areas, where they may even occur at exactly the same site, exemplified also by the sampling locality #138 in mainland Spain (Table S1).

The morphological similarity of G. tigrina and G. sinensis and their co-occurrence in Western Europe, implies that records of presumed G. tigrina simply based on external appearance can no longer be considered to be reliable. In such cases where the animals only show asexual reproduction through the process of fission, as frequently being the case with introduced populations, discrimination between specimens of G. tigrina and G. sinensis can be achieved only through molecular markers. In that respect, it is noteworthy that earlier reports mentioned presumed G. tigrina for Australia (see Sluys et al., 1995a, b and references therein), while the individuals that we analysed from Queensland and Tasmania unequivocally ranked as representatives of G. sinensis.

Planarians are poor dispersers because they lack protection against water loss (Ball, 1974), and in Dugesiidae the cocoons (capsules containing fertilised eggs and nurse cells) are mostly stalked, cemented and usually left attached by an endplate to the under surfaces of stones, fallen leaves or other objects (Sluys & Riutort, 2018), rendering their biochore dispersal difficult, if not impossible (Ball, 1974). This fact, in combination with the widely separated geographic localities where all three species were introduced, strongly suggest that there were multiple introductions and that these were due to human-induced dispersal events. Most likely, the latter were effectuated through the international trade in aquatic plants and the activity of aquarists or through the importation and subsequent culturing of North American strains of G. tigrina and G. dorotocephala as model organisms in education and scientific studies all around the world. The fact that G. tigrina is very sticky, clearly enhances its attachment to all kinds of surfaces and thereby facilitates its passive dispersal (Stocchino et al., 2019). Wright (1987) suggested that multiple occasions of disposal of the contents of aquaria and of small ornamental garden ponds into local rivers might underlie the pattern of distribution of G. tigrina in Great Britain. In addition, disposal of such contents into waste-water pipes might facilitate entrance of flatworms into a river via surface water drains (Wright, 1987). The present data does not allow us to infer whether introductions into different countries came directly from North America or may have occurred in a stepping-stone model from one introduction site to the next. Most probably, both types of events have contributed to the spread of alien species of Girardia. In this process, Girardia alien species may even have been reintroduced to America. Future analyses of more variable markers, as for instance RADseq data and population structure analyses that allow fine detection of haplotype frequencies and potential directions of gene flow, may be necessary for determining the various routes that may have been followed by the three North American species of Girardia in their conquest of the world.

What makes Girardia species successful colonisers?

Our niche modelling indicates that G. tigrina, G. dorotocephala, and G. sinensis do not only have the ecological potential to successfully survive at the localities where they have been introduced but are also capable of expanding their distributional range. Interestingly, in the MaxEnt model the responses of the three species to different environmental variables indicated a key feature explaining their high capacity for colonization, viz., shared tolerance for anthropogenic habitats, such as artificial freshwater habitats (like canals) or freshwater habitats in cultivated and managed areas (Fig. 3). The importance of this variable has been reported also for many other introduced plant and animal species (Rickart et al., 2011; Salomidi et al., 2013; Simkanin et al., 2013; Johnston et al., 2017; González-Ortegón & Moreno-Andrés, 2021), including terrestrial planarians (Álvarez-Presas et al., 2014) and is also one of the main factors explaining the distribution of another introduced freshwater planarian, Dugesia sicula (Leria et al., 2022). In that respect, it is noteworthy that the factor “tolerance for anthropogenic habitats” played no role in the explanation of the autochthonous distribution of D. subtentaculata in the Iberian Peninsula (Leria et al., 2022). Anthropogenic habitats generally are characterised by low autochthonous species diversity, thus offering empty niches for invasive species that can use these as a springboard to natural habitats (Dietz & Edwards, 2006).

Apart from association with humans and anthropogenic habitats, other traits of successful colonisers that may become invaders, are, for example, high abundance in the native range, short generation time, high genetic variability, wide physiological tolerance, phenotypic plasticity, and asexual reproduction (Kolar & Lodge, 2001). These are general features that may vary or combine in different ways, depending on the intrinsic characteristics of each group of organisms and of the receiver area. Our MaxEnt analysis allowed us to determine that the introduced species of Girardia present some of these traits.

A case of wide physiological tolerance concerns temperature, with the three Girardia species presenting a range of highly suitable temperatures that is much broader than in other freshwater planarian species (Vila-Farré & Rink, 2018; Leria et al., 2022). For example, Polycelis felina exhibits an optimality maximum at 4°C, with the optimality rapidly decreasing at temperatures above 10°C, while D. subtentaculata shows a narrow temperature range between 12 and 15°C with a high suitability (Leria et al., 2022). In contrast, in Girardia high suitability was found for temperatures ranging from around 4°C to roughly 20°C. This minimum value of about 4°C found for Girardia in our niche modelling, coincides with ecological studies that found temperature to be a limiting factor in the expansion of G. tigrina in the United Kingdom and, curiously, catalogued it as a warm-water species, with a lower tolerance limit of 6°C for feeding and 10°C for asexual reproduction (Wright, 1987 and references therein). On the other hand, in the present study, G. dorotocephala has been found in Mexico and G. sinensis in Cuba in waters with temperatures ranging around 19.5–20.5°C and 24°C, respectively.

Our niche modelling also revealed that the three species are suitable to a wide range of lithological classes and water discharge regimes, as well as precipitation regimes (Fig. 3). However, for slope, not all of the introduced species presented a flexible suitability, as only G. tigrina would be highly suited to an extremely wide range of terrain slopes. This is surprising, since it has been reported that G. tigrina generally occurs in the lower reaches of lowland rivers with low slopes (Wright, 1987). All of these results suggest that each of the three introduced Girardia species exhibits high tolerances to most of these key environmental variables.

As remarked above, asexual reproduction is considered to be a trait enhancing successful colonization and invasion. For G. tigrina, and G. dorotocephala, it has been established that they show both sexual and asexual reproduction in their native areas, with populations being either strictly sexual or asexual or showing both reproductive modalities (Kenk, 1937, 1972; Puccinelli & Deri, 1991; Knakievicz et al., 2006). For G. sinensis, in the laboratory it was established that the species can reproduce both sexually and asexually by fission (Chen et al., 2015). However, most of the introduced populations are constituted by fissiparous individuals, as evidenced by the presence of individuals with an anterior or posterior blastema, and the absence of cocoons at the various localities. Thus far, only five immigrant sexual populations are known from Western Europe, located in Great Britain, Spain (Menorca Island), Italy, and France (Stocchino et al., 2019 and references therein), while all other known populations reproduce by fission. This prevalence of asexual populations might be advantageous for dispersion and population growth (Wright, 1987; Gee et al., 1998) and has already been proposed as a major factor explaining the successful recent colonisation of D. sicula in the Mediterranean region (Lázaro & Riutort, 2013).

Another important characteristic of planarians possibly contributing to their colonisation capability concerns their top position in the trophic web. Although there is a large number of species known to predate on freshwater planarians, the impact of predation on natural populations of flatworms generally is rather low (Vila-Farré & Rink, 2018). Evidently, this may promote rapid population growth, which, in turn, fulfils a basic requirement for high rates of spread.

Potential ecological impact of introduced Girardia

Since freshwater communities are highly sensitive to composition changes through the introduction of invasive species (Havel et al., 2015; Gallardo et al., 2016), prediction and risk assessment of potential invader species are vitally important for the protection of these habitats. However, risk assessment is distinct from prediction whether a certain species will become an invader, and even more difficult to achieve than the prediction. With respect to freshwater planarians, ecological studies on the specific impact of introduced species on their receiver freshwater communities are (1) limited, (2) focused on only small geographic regions, or (3) are even contradictory. For example, while some studies concluded that G. tigrina may be a competitive invader (Reynoldson, 1985; Wright, 1987; Gee & Young, 1993), its ecological impact on native communities was downplayed in a recent paper (Ilić et al., 2018).

For making predictions whether certain non-native species will become invasive, it is unfortunate that there is no suite of characteristics that can unequivocally signal successful invaders (Kolar & Lodge, 2001). Nevertheless, a non-native species must go through five phases to become an invader: (1) transport to a new region, (2) release or escapement to the wild, (3) establishment, (4) dispersal or spread, and (5) integration or impact (García-Berthou, 2007). From that perspective, all three species of introduced Girardia have already reached phases four and five, albeit that their ecological impact has not been properly analysed.

In the present work, we have run a preliminary study on the potential impact of introduced species of Girardia on autochthonous freshwater planarians that was restricted to the Iberian Peninsula, for which there is a good knowledge about its autochthonous species diversity. Most native species in this region occur at high altitudes Crenobia alpina (Dana, 1766), P. felina (Dalyell, 1814), and various species of Phagocata Leidy, 1847, while others, such as Phagocata ullala Sluys, 1995, are endemic to very restricted areas in which Girardia has never been detected (Baguñà et al., 1980, 1981; Sluys et al., 1995a, b). Nonetheless, D. subtentaculata (Draparnaud, 1801), the most widely distributed autochthonous species, has been found together with Girardia representatives. At locality #247, D. subtentaculata coexists with G. sinensis, while in the north of the Iberian Peninsula, it co-occurs with unidentified Girardia individuals (M. Vila-Farré, pers. comm.) and with G. dorotocephala (recently identified by DNA barcoding in our laboratory). Moreover, optimal environmental conditions for D. subtentaculata were recently estimated in a niche modelling analysis (Leria et al., 2022). For this species, optimal conditions are present along the northern and western coasts of the Iberian Peninsula, a region with principally deciduous tree cover, an annual average temperature around 14°C, and relatively high precipitation regimes. These environmental conditions partially overlap with the optimal ecological niche for G. sinensis, as identified in the present study (Fig. 4). Moreover, niche modelling showed that the Iberian Peninsula presents extensive unsuitable areas for D. subtentaculata (zero suitability values) (Leria et al. 2022), while in our modelling analyses of Girardia, almost all regions of the entire Peninsula received high suitability values (> 0.7), especially for G. tigrina and G. sinensis. In that light, it has been shown that introduced species of Girardia may exert a negative impact on native planarian species, for instance by competing for the same food. Gee & Young (1993) showed that introduced G. tigrina had considerable dietary overlap with the native planarian species Polycelis tenuis Ijima, 1884 and P. nigra (Müller, 1774).

Our data also evidences that the genetic diversity within and among the three introduced species is quite low, as reflected by the short branch lengths in these lineages and, thus, suggests recent speciation, probably after the closure of the Isthmus of Panama (for a discussion, see Benítez-Álvarez et al., 2023a). On the other hand, in the Western Mediterranean more than 10 species of Dugesia have been characterised, representing a broad diversification that evolved during the last 25 My (Dols-Serrate et al., 2020; Leria et al., 2020, 2022; Benítez-Àlvarez et al., 2023b). Therefore, genetic diversity of introduced populations of Girardia is much lower than in the autochthonous species of the genus Dugesia in Western Europe. If changes in ecosystems induced by human activities favour alien species over autochthonous ones, we run the risk that a highly diverse and locally adapted fauna of freshwater planarians in Western Europe will be overtaken by more generalist and genetically low-diverse invaders. Evidently, the ecological impact of introduced Girardia may not be restricted to other planarians but may concern also other freshwater organisms that they may predate or compete with.

Future expansion trends of introduced Girardia

Our niche modelling indicated that at present the three introduced Girardia species have the ecological potential to expand their distributional range (Online Resources Fig. S1), while modelling of future situations, regardless of the particular climatic scenario, pointed to an extension of their potential areas of suitability (Fig. 5 and S2). Their potential for range expansion was mostly restricted to suitable areas in the Northern Hemisphere, in particular to regions throughout Europe, but without diminishing high suitability of regions in the south. In this light, and also considering the generalist characteristics of the species and their demonstrated suitability for human-modified habitats, it is only to be expected that in our current increasingly anthropogenic world they may have an advantage over native species.

From this perspective, it is especially worrisome that G. sinensis has recently been detected for the first time in Africa (Benítez-Álvarez et al., 2022). It is noteworthy that this concerned a region for which our niche modelling already indicated its suitability and that harbours a high diversity of autochthonous species of Dugesia (Leria et al., 2022; Benítez-Álvarez et al., 2023b). Therefore, it will be important to avoid future introduction of G. dorotocephala, an invasive species that has an even wider predicted suitable region than G. sinensis, covering the entire Atlas Mountains, even in future suitability predictions (Fig. 5).

Measures to avoid future introductions may involve creating awareness among aquarists and the staff of educational and research institutions using planarians as model organisms that their culturing waters should not be disposed of by simply pouring these into the sink, or, even worse, directly into the environment. Such water should first be fully cleansed from its inhabitants. Moreover, the implementation of regular surveys, for instance using eDNA techniques, now feasible taking profit of the nuclear markers used here, especially on those areas shown to present more suitability for introduced species, could be a way to detect the new introductions and avoid rapid expansions of the alien Girardia.

Final remarks

In the present work, we do a preliminary dissection of a long-held invasion. We corroborate that human activities have effectuated great changes in the geographic distribution of three invasive species of Girardia in a relatively short period of time (Benítez-Àlvarez et al., 2023a). We characterise the features that may explain such successful colonisation, especially as compared to some autochthonous species, concluding that there exists an overlap between autochthonous and Girardia suitability areas, with an advantage in area extension and physiological plasticity for the introduced species. Moreover, we demonstrate the usefulness of two molecular markers to rapidly and cheaply identify the introduced species at any locality, which may facilitate a continuous survey of freshwater habitats for new introductions and expansions. Even taking into account the limitations in our sampling that may have slightly affected our model accuracy, niche modelling demonstrates its value to predict potential areas with a higher risk of being colonised, which may help focus future projects to monitor the expansion and potential detrimental effects of the three species.