# Hammerhead worms everywhere? Modelling the invasion of bipaliin flatworms in a changing climate

Yoan Fourcade, Leigh Winsor, Jean-Lou Justine

► To cite this version:

Yoan Fourcade, Leigh Winsor, Jean-Lou Justine. Hammerhead worms everywhere? Modelling the invasion of bipaliin flatworms in a changing climate. Diversity and Distributions, 2022, 28 (4), pp.844-858. 10.1111/ddi.13489. hal-03607897

HAL Id: hal-03607897

https://hal.u-pec.fr/hal-03607897

Submitted on 14 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Hammerhead worms everywhere? Modelling the invasion of bipaliin flatworms in a changing climate

Yoan Fourcade<fnref n="1" /> ID | Leigh Winsor<fnref n="2" /> ID | Jean-Lou Justine<fnref n="3" /> ID

**Correspondence**

Yoan Fourcade, CNRS, IRD, INRAE, Institut d'écologie et des sciences de l'environnement, IEES, Univ Paris Est Creteil, Sorbonne Université, F-94010 Créteil, France.

Email: yoan.fourcade@u-pec.fr

**Editor:** José Brito

**Abstract**

**Aim:** Alien species introduced into new ecosystems occasionally predate upon or out-compete native species. Land planarians (Geoplanidae) are a family of carnivorous Platyhelminthes among which several species are found outside their native range. Specifically, hammerhead flatworms originate from Asia and Madagascar but have now reached many new locations worldwide through the transport of exotic plants. Because they are predators of earthworms and snails, they are considered a potential threat to native ecosystems. In this context, to anticipate their potential impacts and to inform early preventative actions, it is necessary to know where these species could spread to in future, or where they might already be present but undetected.

**Location:** Worldwide.

**Methods:** Here, we used occurrence records from online databases and climatic and soil variables to model the potential distribution of five hammerhead flatworm species (*Bipalium adventitium*, *B. kewense*, *B. pennsylvanicum*, *B. vagum* and *Diversibipalium multilineatum*) that are known to occur outside their native range.

**Results:** We demonstrate that precipitation is an important factor determining their distribution, which is in accordance with their known affinity for humidity. We show that some areas have the potential to be invaded by all five species, including regions that are relatively spared so far. This includes the River Plate basin in South America, which already harbours a diverse fauna of native terrestrial flatworms and which appears to also be a potential hotspot for the establishment of alien hammerhead flatworms. According to scenarios of future climate change, two species (*B. kewense* and *B. vagum*) that currently have the largest observed global range are predicted to further increase their potential distribution.

**Main conclusions:** The results we report can be used to provide guidance for monitoring the potential sources of introduction of alien hammerhead flatworms in regions that are suitable, but which are not yet colonized.

**KEYWORDS**

*Bipalium*, *Diversibipalium*, *Geoplanidae*, invasion, Platyhelminthes, species distribution modelling

<footnote n="1">¹Univ Paris Est Creteil, Sorbonne Université, CNRS, IRD, INRAE, Institut d'écologie et des sciences de l'environnement IEES, Créteil, France</footnote>

<footnote n="2">²College of Science and Engineering, James Cook University, Townsville, Queensland, Australia</footnote>

<footnote n="3">³ISYEB Institut de Systématique Évolution Biodiversité, UMR7205 MNHN, CNRS, EPHE, UPMC, Université des Antilles, Muséum National d'Histoire Naturelle, Paris, France</footnote>
# 1 | INTRODUCTION

The last few decades have seen an enormous increase in the amount and frequency of movements of people and goods on a global scale (Ascensão & Capinha, 2017). Trade and transport networks that connect distant regions of Earth contribute to the dispersal of living organisms from their native ranges to new areas (Banks et al., 2015), where they sometimes establish viable populations when the local environmental conditions resemble the conditions to which they were adapted in their native range (Richardson et al., 2000). When non-native species spread at the expense of the local fauna and flora (i.e. when they become invasive), native communities can be disturbed to the point that outcompeted or predated species may critically decline or even be locally extirpated (Doherty et al., 2016; Maron & Marler, 2008). Invasive species may therefore pose a severe threat to biodiversity and ecosystem processes (Tylianakis et al., 2008), such that monitoring their current distribution and impact, as well as predicting their potential spread, is essential to implement adequate mitigation strategies (Gallien et al., 2010; Simberloff et al., 2005).

Land planarians (Geoplanidae) are a family of Platyhelminthes that are entirely free-living carnivorous flatworms, and are apex predators in their microhabitat (Sluys, 2016). There are about 900 described species of terrestrial flatworms worldwide, and many more species await description. Most terrestrial flatworms occur naturally in the tropics and in the temperate Southern Hemisphere, in a large range of habitat types, including from mesophile to xenophile habitats, from alpine herb fields to sandy semi-desert and from subantarctic rata forest to tropical rain forest (Winsor et al., 1998). However, some species are found outside their native range in large parts of the world such as Europe (Cannon et al., 1999; Čapka & Čejka, 2021; Carbayo et al., 2016; Jones, 2019; Jones et al., 2020; Jones & Sluys, 2016; Justine et al., 2014, 2020, 2022; Justine, Winsor, et al., 2018; Mori et al., 2022), the Americas (Justine et al., 2015, 2019, 2021), Asia (Chaisiri et al., 2019; Hu et al., 2019) and Polynesia (Justine, Lemarcis, et al., 2018; Justine, Winsor, et al., 2018), and some of them occasionally became invasive (Sluys, 2016).

Among land planarians, hammerhead flatworms, that is species belonging to the subfamily Bipaliinae and mainly to the genera *Bipalium* and *Diversibipalium*, originate from Asia and Madagascar, but have now reached many new locations worldwide (Ducey et al., 2007; see, e.g., Ducey & Noce, 1998; Justine, Winsor, et al., 2018; Morffe et al., 2016). *Bipalium kewense*, for example, which has a natural range in Southeast Asia, was already present on all continents except Antarctica as early as 40 years ago (Winsor, 1983). These species have mainly been introduced by the transport of potted exotic ornamental plants in soil and as a consequence are mostly found, in their non-native range, in hothouses, gardens and urban parks. Since there is evidence that hammerhead flatworms are able to survive and reproduce in non-native regions, there are concerns that they may expand and establish in natural ecosystems, or even that they may already be present but undetected so far.

Although the invasive nature of hammerhead flatworms is not always immediately obvious, expansion of their range is worrying because they are predators of earthworms or molluscs (Fiore et al., 2004; Winsor et al., 2004). Earthworms play a crucial role in ecosystems, especially for the regulation of many ecological and physical processes in the soil (Liu et al., 2019). Terrestrial molluscs, on the contrary, are among the most threatened species groups, totalling more than one third of all extinctions recorded to date (Lydeard et al., 2004). Therefore, the global expansion of hammerhead flatworms raises concerns about their potential effect on natural ecosystems. This is especially worrying since other closely related species have documented impacts on native communities. One of them is *Arthurdendyus triangulatus*, the New Zealand flatworm, which is a predator of lumbricid earthworms, invasive in the United Kingdom, where it has been shown to decrease populations of anecic earthworm species (Jones et al., 2001; Murchie & Gordon, 2013). Similarly, *Platydemus manokwari* is a flatworm species that has been introduced since the 1950s in several tropical islands in the hope that it would control an invasive predatory snail. Since then, however, *P. manokwari* has caused the decline or even the extinction of several native snail species (Gerlach, 2019; Gerlach et al., 2021; Sugiura, 2009; Sugiura & Yamaura, 2009).

To anticipate the possible impacts of hammerhead flatworms on native communities of earthworms and molluscs, and in order to inform early prevention actions, it is necessary to know where species could spread to in future, or where they might already be present but remain undetected. Without doubt, human transport is the primary driver of current flatworms' distribution in non-native areas (Justine, Winsor, et al., 2018). However, other environmental factors such as climate or soil properties may play a role in the probability of establishment at a given location, following transport and release into nature. In this regard, correlative species distribution models (SDMs), or ecological niche models (ENMs), are statistical tools that can link species occurrences to environmental variables in order to model ecological niches and project their suitability in space and time (Elith & Leathwick, 2009). The applicability of SDM to invasion biology lies in the hypothesis that species conserve their climatic niche (Barbet-Massin et al., 2018). Under this assumption, species occupy the same range of climatic conditions in their invasive range as in their native range, in such a way that the distribution of climatically suitable areas globally, currently and in future is predictable from the distribution of the species as presently observed. However, even if this assumption is not met—there is contradictory evidence regarding the prevalence of niche conservatism in invasive species (Atwater et al., 2018; Liu et al., 2020; Petitpierre et al., 2012)—it remains possible to produce useful predictions of the potential distribution of invasive species as long as enough occurrences are sampled throughout the entire occupied range, and providing that models are properly parametrized (Elith et al., 2010; Fourcade, 2021).

The aim of the present study was to predict the global distribution of suitable areas of five species of hammerhead flatworms, which are considered potentially invasive because they have been introduced in non-native areas and have successfully established
populations. Species distribution models, trained from georefer-
enced species occurrences from online databases and climate and
soil variables, will be used to answer the following questions: (1)
What environmental factors determine the global niche of terrestrial
flatworms, and are there differences between species? (2) Where
are currently distributed the suitable habitats of these species, that
is what regions could each species potentially invade if they were in-
troduced by human transport? (3) What regions are the most suscep-
tible to invasion of terrestrial flatworms? (4) Do existing scenarios of
climate change predict spatial shifts in the distribution of species'
environmental suitability in future that would modify the areas of
potential invasion? The suitability maps that our models produce
provide a spatially explicit assessment of the potential spread of
alien hammerhead flatworms that can be used as an early warning of
their possible expansion in natural environments.

2 | METHODS

2.1 | Bipaliin land planarian study species

We focus in this study on five bipaliin species that are known to have
spread outside their native range. *Bipalium adventitium* Hyman, 1943
was first identified in gardens in 1943, in Berkley and Pasadena,
California, United States. Its most likely countries of origin are Japan
and/or South Korea where similar native species exist. This species
predates exclusively on earthworms. *Bipalium kewense* Moseley,
1878 was described from a specimen found in 1876-1877 in one
of the hot houses of the Royal Botanic Gardens, Kew, England,
United Kingdom. It originates from the central highlands of Vietnam
and is known to predate on earthworms. So far, all specimens col-
lected from several continents were genetically identical for the
cox1 gene widely used for barcoding, suggesting that a single animal
and its clones are responsible for all documented invasions (Justine,

Winsor, et al., 2018). *Bipalium pennsylvanicum* Ogren, 1987 was first observed in 1978 on the St David's campus of Eastern College, Delaware, Pennsylvania, United States (Ogren, 1986). It is generally thought that its country of origin is Japan, where an identical species was observed near Nagoya in 2018. It is a predator of earthworms. *Bipalium vagum* Jones and Sterrer, 2005 was first found in 2003 in the Aberfeldy Nursery, Paget Parish, Bermuda, and most likely originates from Southeast Asia, possibly Thailand, where similar species occur. Contrary to the other study species, *B. vagum* predates on land snails (Ducey et al., 2007). Finally, *Diversibipalium multilineatum* (Makino and Shirasawa, 1983) was first found in an urban garden in Hino, Tokyo, Japan. It is possibly native to Japan, though it is known there only from non-sexual specimens suggesting that its country of origin may be elsewhere, such as China where there are similar species. This species is known to predate on earthworms. Brief descriptions and illustrations of these species are available in Wallace and Winsor (2020), and their known distribution from observation records is described in Figure 1 (see Section 2.2).

2.2 | Occurrence data

As baseline data to train species distribution models, we used geo-
referenced occurrences of each of the five species, obtained from
two public databases of biodiversity data (Figure 1 and Table 1):
the Global Biodiversity Information Facility (GBIF, www.gbif.org)
and iNaturalist (www.inaturalist.org). The GBIF database compiles
observation records from a variety of sources (Edwards, 2004), in-
cluding museum specimens, scientific surveys, opportunistic records
or citizen science programmes (see, e.g., the citizen science project
aimed at collecting observations of land planarians in France: http://
bit.ly/Plathelminthe). iNaturalist is a citizen science project built as a
social network where users can easily record species' observations
through a web portal or a mobile application. One interesting feature

<figure><img src="image_2.png" /><figcaption>FIGURE 1 Maps of species occurrences obtained from GBIF and iNaturalist databases, after spatial thinning (see Table 1 for the detail of the number of occurrences)</figcaption></figure>
of iNaturalist is that observations are accompanied by photographs, which helps the taxonomic identification of the observed species via both an AI and an expert review by experienced users. The quality of data is subsequently classified from “casual” to “research grade” depending on the outcomes of the community identification process.

As the first step of this project, two of us (LW and JLJ), acting as members of iNaturalist, examined photographs and edited more than 6000 geoplanid records worldwide, removing inconsistencies and misidentified species, to clean the database before downloading records for further analysis. We downloaded occurrences from GBIF.org (2021) and iNaturalist through R (data accessed on 1 October 2021), using the “spocc” package (Chamberlain, 2020). We removed occurrences that had obvious georeferencing errors and discarded duplicated coordinates. Additionally, we subjected each occurrence to a careful screening to discard suspicious data such as species recorded in regions where they are very unlikely to occur. To reduce the effect of sampling bias that occurs when all regions are not surveyed with the same effort (Fourcade et al., 2014; Kramer-Schadt et al., 2013), we performed a spatial thinning of occurrences, keeping only one coordinate per grid cell in a raster of 50 × 50 km resolution. This ensures that the final dataset does not contain clusters of high density of occurrences that could bias model predictions.

<table><caption>TABLE 1 Number of occurrences obtained from GBIF and iNaturalist for each study species, before and after spatial thinning at a 50-km scale</caption><thead><tr><th rowspan="2">Species</th><th colspan="2">No. of unique occurrences in raw data</th><th colspan="2">No. of unique occurrences in spatially thinned data</th></tr><tr><th>iNaturalist</th><th>GBIF</th><th>iNaturalist</th><th>GBIF</th></tr></thead><tbody><tr><td>Bipalium adventitium</td><td>761</td><td>392</td><td>220</td><td>0</td></tr><tr><td>Bipalium kewense</td><td>3038</td><td>1515</td><td>655</td><td>42</td></tr><tr><td>Bipalium pennsylvanicum</td><td>161</td><td>79</td><td>43</td><td>0</td></tr><tr><td>Bipalium vagum</td><td>663</td><td>339</td><td>240</td><td>12</td></tr><tr><td>Diversibipalium multilineatum</td><td>93</td><td>66</td><td>48</td><td>11</td></tr></tbody></table>

Note: The latter were used to train distribution models, and their location is shown in Figure 1.

Here, since we carefully checked and edited all iNaturalist records, we prioritized these data in the spatial thinning process such that in the end, the final dataset essentially consists of iNaturalist records complemented by GBIF data (Table 1).

## 2.3 | Environmental variables

Although we did not know beforehand what factors drove the distribution of the five study species, there is some evidence that land planarians are sensitive to variation in temperature (Daly & Matthews, 1982; Negrete et al., 2020) and soil moisture (Boag et al., 1998). Therefore, species distribution models were trained using as predictors six environmental variables that represent both climatic conditions and soil properties (Table 2). Among climatic variables, we initially chose six bioclimatic variables from the CHELSA climate dataset (chelsa-climate.org, Karger et al., 2017), which accounted for the mean and extreme temperature and precipitation: mean annual temperature (BIO1), maximum temperature of the warmest month (BIO5), minimum temperature of the coldest month (BIO6), annual precipitation (BIO12), precipitation of the wettest month (BIO13) and precipitation of the driest month (BIO14). Because there was evidence of high correlation between these variables, we discarded BIO1 and BIO12 so that the variance inflation factor among all remaining variables was <5. As soil variables, we chose soil pH in water (Hengl, 2018) and soil organic carbon (Hengl & Wheeler, 2018) at the surface (0 cm), both estimated through machine learning in the OpenLandMap project (openlandmap.org). All environmental layers were obtained as raster files with a resolution of 30 arc second and aggregated to 5 arc minutes for faster computation and because we aimed to conduct broad analyses at the global scale, and not necessarily to provide small-scale predictions that could be dependent on other, non-accounted factors.

One objective of the present study was to quantify changes in the amount and spatial distribution of suitable hammerhead flatworm habitats that may occur in future as a result of climate change. For this purpose, we projected our distribution models into future conditions, according to various scenarios of climate change. In this approach, we assumed that soil properties will remain the same, both because there is no estimate of future soil conditions and because

<table><caption>TABLE 2 Summary of the environmental variables that were included as predictors in distribution models</caption><thead><tr><th>Abbreviation</th><th>Meaning</th><th>Unit</th><th>Source</th><th>VIF</th></tr></thead><tbody><tr><td>Bio 5</td><td>Maximum temperature of the warmest month</td><td>°C</td><td>CHELSA (http://chelsa-climate.org/)</td><td>2.61</td></tr><tr><td>Bio 6</td><td>Minimum temperature of the coldest month</td><td>°C</td><td>CHELSA (http://chelsa-climate.org/)</td><td>3.29</td></tr><tr><td>Bio 13</td><td>Precipitation of the wettest month</td><td>kg/m<sup>2</sup></td><td>CHELSA (http://chelsa-climate.org/)</td><td>2.19</td></tr><tr><td>Bio 14</td><td>Precipitation of the driest month</td><td>kg/m<sup>2</sup></td><td>CHELSA (http://chelsa-climate.org/)</td><td>1.55</td></tr><tr><td>pH</td><td>Soil pH in H<sub>2</sub>O at 0 cm</td><td>-</td><td>OpenLandMap (https://zenodo.org/recor<br>d/2525664)</td><td>1.76</td></tr><tr><td>SOC</td><td>Soil organic carbon</td><td>%</td><td>OpenLandMap (https://zenodo.org/recor<br>d/2525664)</td><td>1.14</td></tr></tbody></table>

Note: The variance inflation factor (VIF) indicates the level of correlation between variables.
it turned out that soil variables contributed little to species distribu-
tions (see Section 3). We obtained future predictions of the same four
bioclimatic variables as used during modelling for the years 2041–
2060 (referred to as 2050 in the figures) and 2061–2080 (referred
to as 2070 in the figures) according to two scenarios of greenhouse
gases emissions (RCP 4.5 and RCP 8.5) and five global circulation
models (GCMs): ACCESS1.0, CESM1-BGC, CanESM2, HadGEM2-AO
and MIROC-ESM. The RCP (Representative Concentration Pathway)
8.5 is the most pessimistic climate change scenario, which assumes
that greenhouse gases emissions will continue to rise throughout the
21st century, to the point that it overestimates the amount of fossil
fuels that can be realistically exploited (Burgess et al., 2020). It is
thus viewed as a worst-case scenario that is unlikely to happen (but
see Schwalm et al., 2020). On the contrary, RCP 4.5 is an optimistic
scenario that assumes temperature rise will level off prior to the end
of the 21st century thanks to the global adoption of effective ac-
tions to mitigate greenhouse gases emissions (Thomson et al., 2011).
Projecting changes in habitat suitability under both these scenarios
is therefore a way to explore a range of possible trends in climate
change when the actual fate of the global climate remains uncertain.

## 2.4 | Species distribution modelling

To predict the potential distribution of hammerhead flatworms,
we adopted an ensemble modelling framework, which consists in
producing species distribution models using different approaches
and combining all outputs into an averaged model. This strategy is
expected to provide more accurate results than individual models
when the relative performance of the different modelling algorithms
available is highly context-dependent and cannot be assessed be-
forehand (Hao et al., 2020). We trained species distribution models
and combined them using the ensemble modelling strategy imple-
mented in the "biomod2" R package (Thuiller et al., 2009).

Here, we considered four different distribution model algo-
rithms that are routinely used in similar studies: generalized linear
model (GLM), boosted regression trees (BRT, Elith et al., 2008),
maximum entropy (MaxEnt, Phillips et al., 2006) and random for-
est (Breiman, 2001). Automatic variable selection in GLMs was
performed based on the AIC criterion, and we allowed quadratic
relationships and one-way interactions between variables. In BRT,
we set the number of trees to 2500 and the interaction depth to
7. In MaxEnt models, we allowed only linear, quadratic and hinge
features to obtain relatively smooth response curves, and we set
the regularization parameter to three instead of the default value
of one to prevent overfitting. Finally, random forest models were
trained with 5000 trees.

When only presence data are available, as is the case here, the
environment at occurrence points must be compared with the en-
vironmental values in a set of background points representing the
available conditions, instead of true absence points. Ideally, the
background must be sampled within the areas that could have po-
tentially been reached by the species (Barve et al., 2011). When this

area is not known, and especially for species that are not at equilib-
rium with their environment, a common strategy is to restrict the
selection of background points in space so that they remain at rela-
tively close distance to species' records. Here, we sampled for each
species 10000 background points within 500-km buffers around
the occurrence points.

During the modelling process, each species occurrence dataset
was divided into four geographical subsets corresponding to the
four corners of occurrences' coordinates separated by latitude and
longitude mid-points. Each subset was used in turn to evaluate the
predictive performance of the model, while the other three were
used as training data, resulting in four different models per algo-
rithm per species. This ensures that models were evaluated using a
spatially independent dataset (Bahn & McGill, 2013; Fourcade et al.,
2018). Model performance was estimated by the area under the re-
ceiver operating characteristic curve (AUC, Hanley & McNeil, 1982).
The applicability of AUC for SDM studies has been increasingly
questioned in recent years (Fourcade et al., 2018; Jiménez-Valverde,
2012; Lobo et al., 2008), but remains one of the most popular met-
rics to quantify the performance of SDMs, mostly because other
available metrics appeared to suffer from the same issues. Owing
to these criticisms, we did not interpret the AUC in terms of the
probability that our models produce accurate predictions of species
distributions. However, we relied on it to rank the different models
in the ensemble approach and to eliminate those that were clearly
unable to fit our data (see below). Evaluation metrics are reported in
Appendix 1, Table S1.

For each species, we produced ensemble models as the mean
of all algorithms and data subsets weighted by their AUC, dis-
carding models with AUC <0.5. The resulting average model rep-
resents a global map of predicted habitat suitability on a scale of
0 (unsuitable) to 1 (fully suitable), at the same resolution as the
input environmental variables (5 arc-min). All models were pro-
jected into global current and future conditions; the latter were
averaged across all circulation models, resulting in 4 future predic-
tions per species: 2041–2060 and 2061–2080 according to RCP4.5
or RCP8.5 scenarios.

In addition, we report for each model the relative importance
of each variable evaluated using permutation, as well as response
curves, to represent the relationships between species' presence
and environmental variables that underlie the models. Species dis-
tribution modelling is summarized in a ODMAP protocol (Appendix
S1), as recommended by Zurell et al. (2020).

## 2.5 | Post-modelling analyses

When models are projected in space and time, it is possible that
model predictions are extrapolated in areas where the combina-
tion of environmental variables does not exist in the training region.
To quantify where and to what extent models extrapolate in novel
environments, we computed multivariate environmental similarity
surface (MESS, Elith et al., 2010) maps using the "dismo" R package
(Hijmans et al., 2020). This analysis compares the environment en-
countered in the training data (presence + background) with the en-
tire environmental rasters (since models are projected at the global
scale), separately for current and each future condition, to produce
maps of environmental similarity, highlighting regions of extrapola-
tion (shown in Appendix S4).

We converted each continuous suitability map into binary pres-
ence/absence maps, using the suitability threshold that maximizes
the sum of sensitivity and specificity, as recommended by Liu et al.
(2013). For future climatic conditions, since we obtained predic-
tions for 5 different GCM, we considered as suitable the pixels that
were above the threshold for at least three models. These binary
predictions were summed to produce predicted species richness
maps in current and future, allowing us to estimate the number of
hammerhead flatworm species that could potentially colonize the
same areas. Note that this is directly derived from merging species-
specific SDMs and does not consider species interactions between
co-occurring hammerhead flatworm species. Moreover, using binary
maps, we also calculated the global potential range area (in km²) pred-
icted for each species in current and future climatic conditions. We
also estimated the mean latitudinal shift (in absolute values since
species are—potentially—distributed across both hemispheres) be-
tween the currently predicted range and model projections in future
climate.

To estimate the degree of overlap between the distributions
of each species, we used Schoener's *D* index (Schoener, 1968), be-
cause it has been shown to be the best descriptor of niche overlap
in the geographical space (Rödder & Engler, 2011). We computed
Schoener's *D* between each pair of species in current and future
climatic conditions so that we could quantify changes in overlap
between the potential distributions of the different species of ham-
merhead flatworms. Calculation was performed by comparing pairs
of ensemble suitability maps using functions in the "dismo" R pack-
age (Hijmans et al., 2020). As such, this approach is complementary
to the mapping of predicted species richness that shows potential
overlap in space but relies on binarized maps instead of the contin-
uous model outputs.

3 | RESULTS

Although there was no single variable that appeared to unambigu-
ously constrain the distribution of all the five species modelled, we
observed that the relative importance of at least one climate variable
was always higher than the soil variables (Figure 2a). Precipitation of
the driest month (BIO 14) and minimum temperature of the coldest
month (BIO 6) were for all species the two most important variables
in the models (with the exception of *B. adventitium* for which BIO 5

was also equally important), both having on average an importance
of ca. 0.38. In terms of precipitation, we observed an increase in
habitat suitability with increasing BIO 14, up to a plateau reached
at various values depending on the species (e.g. ca. 250 mm for
*B. pennsylvanicum* and ca. 400 mm for *D. multilineatum*). Response
curves also showed that precipitation of the wettest month (BIO
13) was positively related to suitability for *B. kewense* and *B. vagum*
(Figure 2b). Species' dependency on temperature took various
forms, but was always characterized by a minimum environmental
suitability for small values of BIO 6, showing that flatworm species
were systematically absent from regions exhibiting cold winter tem-
peratures. Some species, especially *B. adventitium* and *B. kewense*,
also showed a sudden increase in suitability after 25°C of maximum
temperature of the warmest month (BIO 5) (Figure 2b).

Overall, all five species had higher suitability in the current cli-
mate in the same regions: Southeast Asia, New Zealand and the
eastern coast of Australia, south-eastern South America (northern
Argentina, Uruguay, Paraguay and southern Brazil) and north-west
of South America, eastern coast of the USA, western Europe and
central Africa (Figure 3). Differences between species exist, though,
as evidenced by the fact that only small areas are predicted to
host (potentially) all five species together (only a few grid cells in
Uruguay; see species richness map in Figure 3). Converting contin-
uous suitability predictions into binary maps revealed considerable
differences between species in terms of their potential range size
(Appendix S1, Table S2 and Appendix S3). Only *B. adventitium* was
predicted to have >10 × 10⁶ km² of suitable area worldwide, while
on the contrary, the potential range area of *D. multilineatum* and *B.
pennsylvanicum* was estimated at 1.10 × 10⁶ km² and 0.16 × 10⁶ km²
only, respectively. Both *B. kewense* and *B. vagum* had a predicted
range of intermediate size (7.68 × 10⁶ km² and 7.47 × 10⁶ km², re-
spectively). Measures of overlap revealed that *D. multilineatum* and
*B. pennsylvanicum* had the most similar potential distributions in cur-
rent conditions, while *D. multilineatum* and *B. kewense* were the most
dissimilar (Appendix S1, Figure S1).

At the spatial scale of our models, there was no strong shift in
the distribution of suitability for the five species. As such, the poten-
tial hotspots of highest species richness of hammerhead flatworms
remained located in the same regions. However, it was evident that
there was a tendency for a shift towards higher latitudes, especially
in the 2061-2080 period, and for the scenario of maximum climate
change, RCP 8.5 (Figure 4, Appendix S2 and Appendix S1, Figure S3),
except for *B. vagum*. Similarly, overlaps between potential ranges of
the five species did not change much in future climate compared with
current conditions (Appendix S1, Figure S2). However, binary predic-
tions unveiled substantial changes in the size of suitable ranges in
future for most species (Figure 5). Out of the five study species, two
of them (*B. adventitium*, *D. multilineatum*) showed a clear decrease

FIGURE 2 Species response to environmental variables. (a) Variable importance estimated by permutation in each species distribution model. Boxplots represent the five quartiles of importance, while model-specific values are shown as dots (16 values per species per variable since four different algorithms were trained, with four replicates each). (b) Ensemble response curves estimated in each species distribution model by computing predicted suitability across a range of values for each variable while the others were kept at their mean
<figure><img src="image_3.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /><figcaption>FIGURE 3 Predicted suitability in current conditions obtained from ensemble species distribution models for each of the five hammerhead flatworm species. The bottom-right map shows the predicted potential species richness—from 0 (unsuitable for all species) to 5 (suitable for all species)—as the sum of binary predictions obtained after converting suitability maps according to the threshold that maximizes sensitivity and specificity</figcaption></figure>

in their potential range in all scenarios of climate change by ca. 30–45%. On the contrary, it is predicted that the size of suitable area will increase for *B. vagum* (up to ca. +50% in 2061–2080 for the RCP 8.5 scenario) and *B. kewense* (up to ca. +17% in 2061–2080 for the RCP 8.5 scenario). There were large discrepancies in the potential change in range size predicted for *B. pennsylvanicum* in future scenarios: in 2041–2070, models showed on average an increase in the area of suitable range by ca. 25%–40%, while the mean predicted change in 2061–2080 was negative. However, it must be noted that this was associated with large uncertainty depending on the global circulation model; some models predicted a large decrease, while others predicted a substantial increase in suitable area.

4 | DISCUSSION

Invasive species are undoubtedly one of the most severe threats that biodiversity faces in the 21st century (Doherty et al., 2016). Often, the impact of invasive species is discovered only once they have already significantly affected native populations. By then, their removal may be extremely challenging or even unattainable, especially for fast-reproducing species, such that the functioning of invaded ecosystems may be durably altered. Therefore, an important strategy to mitigate the impact of alien species worldwide is to predict in advance the potential spread of the species that may imperil native species if they were allowed to expand in new areas (Gallien et al., 2010). Here, we employed state-of-the-art modelling techniques to predict the global distribution of suitable areas for a group of flatworm species that are already present outside their native range and that deserve particular attention as predators of local fauna.

Aside from suitability maps, the first lesson learnt from our
models resides in the species response to environmental variables

that revealed similarities, but also some differences, between the
five study species. All of them, especially *B. adventitium*, *B. pennsyl-
vanicum* and *D. multilineatum* for which this pattern was particularly
strong, showed an affinity for higher precipitations during the driest
month of the year, demonstrating that dry regions are unsuitable
for them. This is in accordance with observational and experimen-
tal evidence that proves humidity to be a key factor for the spread
and establishment of land planarians. For example, for some spe-
cies of terrestrial planarians in Victoria, Australia, there is a strong
correlation between their coarse-scale distribution (10 x 10 minute
grid cells) and the Bassian (Kosciuskan) and Eyrean zoogeographic
subregions, which can be regarded as moisture zones (Winsor, 1977,
1998, 2003). Ecological and physiological observations of terrestrial
flatworms also confirm their affinity for moist habitats. Terrestrial
flatworms lack any special physiological or anatomical mechanism
for water conservation (Kawaguti, 1932) and therefore rely strongly
on water conserving behaviours for their survival. They require a
humid microhabitat found beneath leaf litter and rotting fallen logs,
from which they emerge when humidity is high (early observations
summarized by Graff, 1889; Schirch, 1929; Steel, 1897). When con-
ditions are too hot and dry, or possibly too cold for them, flatworms
are also known to retreat several centimetres into the soil or deep
into rotten timber (Boag et al., 1998; Cardale, 1940; Darwin, 1844;
Fletcher & Hamilton, 1887; Winsor et al., 2004). Overall, previous
evidence suggests that the global distribution of the study species
is limited to humid regions of the world, while they are absent from
areas that have a dry season that would make them face a risk of des-
iccation, which appears to be confirmed by our modelling approach.

The distribution of our study species was also clearly associated
with increasing temperature of the coldest month (BIO 6), with the
exception of *B. pennsylvanicum* for which this was more ambigu-
ous. This demonstrates that they are unable to establish in regions
2050 (RCP 8.5)

<figure><img src="image_6.png" /><figcaption>FIGURE 4 Predicted suitability in future conditions (top: 2041–2060, median 2050; bottom: 2061–2080, median 2070) under the RCP 8.5 scenario, obtained from ensemble species distribution models for each of the five hammerhead flatworm species. The bottom-left maps show the predicted potential species richness—from 0 (unsuitable for all species) to 1 (suitable for all species)—as the sum of binary predictions obtained after converting suitability maps according to the threshold that maximizes sensitivity and specificity. Predictions for a milder scenario of climate change (RCP 4.5) are presented in Appendix S2</figcaption></figure>

characterized by cold winters. In particular, three species (*B. adventitium*, *B. kewense* and *B. vagum*) appeared to be extremely influenced by temperature variables, while their response to precipitation is less evident. Interestingly, these species also have, according to our models, a much larger potential range than the two species that responded the most strongly to precipitation (*D. multilineatum* and *B. pennsylvanicum*). It is thus possible that, by being more tolerant to drought, *B. adventitium*, *B. kewense* and *B. vagum* could colonize a larger global range and thus be limited in their expansion only by cold winter temperatures. Comparing predicted occurrence with actual occurrence records (see Figure 1), *B. kewense* and *B. vagum* are effectively present in a large part of the globe, especially *B. kewense*

that has been found on all continents but Antarctica. Of course, our interpretation relies on the fact that the occurrences we obtained compare reasonably well with the current distribution of the species; in case all species are in fact present worldwide, observed differences in response to temperature and precipitation by the species could be an artefact of biases in the data (Yackulic et al., 2013).

The fact that we did not detect a strong effect of soil variables on species distributions may be due to the spatial resolution of our analysis; we cannot exclude that at a finer scale, soil characteristics could be an important driver of local distribution within areas that are otherwise climatically suitable (Thuiller, 2013). In this regard, Steel (1897) observed that land planarians were most plentiful
<figure><img src="image_7.png" /><figcaption>FIGURE 5 Change in predicted range area in future, according to different scenarios of climate change, compared with potential range area in current conditions. Range size was obtained after converting continuous suitability predictions into binary presence-absence maps using the threshold that maximizes sensitivity and specificity</figcaption></figure>

in areas where the soil was rich, and he postulated that rich fertile land that retained moisture best supported the land planarians and a better supply of their potential prey, especially in dry weather, than was present in less favourable districts. Moreover, we know that, in the absence of surface cover, terrestrial flatworms can live entirely in soil as permanent members of the soil fauna (Winsor et al., 1998), which suggests that soil properties are likely to play a role in their ecological niche.

A remarkable feature of our distribution models is that they all predict suitable habitats in the same region of South America, corresponding to the River Plate basin (covering Uruguay, north-eastern Argentina, south-eastern Paraguay and southern Brazil) and partly expanding through parts of the Atlantic Forest—even though only *B. kewense*, and *B. vagum* to a lesser extent, was recorded there. This is the only place on Earth where there are areas predicted as suitable for all five species together (predicted potential richness of five species; see Figure 3). These regions have already been shown to host a great diversity of land planarians (Álvarez-Presas et al., 2011, 2014; Carbayo et al., 2002; Iturralde et al., 2021; Negrete, Do Amaral, et al., 2020; Negrete et al., 2019; Sluys, 1999). Moreover, two recent distribution modelling studies (Fourcade, 2021; Negrete, Francavilla, et al., 2020) identified the same area as highly suitable for another flatworm species (*Obama nungara* Carbayo et al., 2016), a predator of earthworms and snails whose expansion in non-native areas (including Europe, Justine et al., 2020) is also a source of concern. This region of South America thus has both the current highest diversity of endemic land flatworms and the higher potential for being colonized by alien species. In the event that these species are effectively introduced in this area, there is a high risk that they may alter its original biodiversity. On the contrary, one of the reasons for the success of alien species in invading new areas, including land flatworms, is the absence of predators and competition (Duncan, 1997; Keane & Crawley, 2002). In this regard, invading South America might be more challenging for hammerhead flatworms than in other areas in the world. Indeed, competition has been described for several species in this area (Boll & Leal-Zanchet, 2016; Boll et al., 2015), and it is thus likely that competition with the numerous native terrestrial

flatworms may also play a role in limiting the invasibility of South America.

In some parts of the globe, the distribution of environmental suitability for the study flatworm species can be directly related to the distribution of our input data: all five species have suitable habitats in western Europe and North America, where most of the occurrence records come from. Interestingly, this is the case even for *B. adventitium* and to a lesser extent for *B. vagum* and *B. pennsylvanicum*, which have not yet been recorded in Europe. Looking more precisely at Europe, we observed that our models predicted suitable habitats for at least three species in the eastern part of France, in the Balkans, and around the Black Sea, which are areas where none of the five study species have been reported so far. Another alien flatworm species (*O. nungara*) has been recorded in eastern France (Justine et al., 2020), Italy (Mori et al., 2022) and Slovakia (Čapka & Čejka, 2021), so far in low numbers, showing that these regions may be suitable for alien terrestrial flatworms. Similarly, New Zealand and eastern Australia appeared to be among the most suitable areas for almost all species, despite the fact that the occurrence dataset we compiled included records in this region for *B. kewense* only. In this regard, it is noticeable that New Zealand is the native range of *Arthurdendyus triangulatus*, a flatworm species that has invaded the Faroe Islands and the British Isles, but which has not established in continental Europe or elsewhere (Murchie & Gordon, 2013). Altogether, this shows that several species of flatworms are adapted to a temperate climate, in such a way that the regions we describe above may potentially have a high invasion risk. This calls for careful monitoring of these areas, where hammerhead flatworms might already occur without our knowledge, or may establish viable populations if they were to be transported there.

Thanks to the availability of climate variables projected under scenarios of climate change, we were able to use our species distribution models trained in current climate to predict what might become the distribution of habitat suitability for hammerhead flatworms in future. At the global scale, the areas predicted to be suitable did not change dramatically; the potential hotspots remained located in western Europe, North America, the River
Plate basin, Southeast Asia, eastern Australia and New Zealand. However, there were some substantial differences worthy of investigation. For example, while there were areas in South America that had a potential richness of the five study species in current conditions, this was no longer the case in most of the climate change scenarios. A careful examination revealed that *D. multilineatum* was predicted to lose all its suitable area in this region in 2061–2070 under the 8.5 RCP (the RCP 4.5 scenario predicted a strong but not complete range retraction). Similarly, our models predicted a large reduction in habitat suitability in South America for *B. adventitium* in all future scenarios. In the light of these results, we may expect that invasion risk will decrease in this area, which appeared to be a potential hotspot for invasive hammerhead flatworm in the current climate. However, since the same region is home to many native flatworm species (Álvarez-Presas et al., 2011, 2014; Carbayo et al., 2002; Iturralde et al., 2021; Sluys, 1999), this also raises concerns for future of the local fauna.

A consistent prediction from our models is a poleward shift of environmental suitability in future, all the more so for scenarios that predict large changes in climatic conditions. A latitudinal shift in species distributions is a classical outcome of temperature warming and has already been observed empirically in a wide range of organisms (Chen et al., 2011). Here, a consequence is that both species that have the largest observed global distribution (*B. vagum* and *B. kewense*) may benefit from climate change to expand their (potential) range even more. Regarding *B. vagum*, models predict a relatively stable suitable area with respect to average latitude but point to an important increase in the size of this suitable area. These species have already been introduced in Europe (Justine, Winsor, et al., 2018; Mori et al., 2022) and North America, where populations have been maintained in the long term (Ducey & Noce, 1998; Justine, Winsor, et al., 2018; Winsor, 1983). Whether they will spread even more because of climate change will thus depend on their intrinsic dispersal ability and on human transport within the invaded range. However, ideally, special attention such as monitoring should be given to areas that may appear unsuitable currently, but which may become suitable in the near future, so as to ensure that no new introductions occur in regions that have been spared until now. Moreover, contrary to the present-day introduced individuals of *B. kewense* that are genetically similar and possibly clonal (Justine, Winsor, et al., 2018), more genetic variability may exist in *B. kewense*'s native area, such that additional variants not yet identified in non-native areas, and thus not fully taken into account here, might in fact have a much higher potential area for invasion.

There were two species for which models unambiguously predicted a reduction in potential range size in future. This concerns both the species that has the largest potential distribution in current conditions (*B. adventitium*) and the species that has the most restricted range (*D. pennsylvanicum*). Both *B. adventitium* and *B. pennsylvanicum* exhibit a unimodal response to a temperature variable (BIO 5 and BIO 6, respectively), with a decrease in suitability for the highest temperatures. The observed reduction in potential range in

a warming climate may thus reflect increased mortality during warm summers or winters—provided that the effect is real and direct, since this variable could also be a proxy for the effect of another factor, non-accounted for. In a management perspective, these findings may indicate that these species could become extirpated from part of their invaded range without further intervention. However, given the multiplicity of factors that could not be considered here, such as human transport, microclimate or biotic interactions, these predictions should not be used to justify a relaxation of measures aimed at preventing the spread or the introduction of these species in non-native areas.

The spread of land planarians outside their native range has received increasing attention recently, including from the public who has been made aware of this issue by popular science publications and citizen science programmes (Justine, Winsor, et al., 2018; Justine et al., 2020). The ability of some species to be rapidly dispersed over all continents has made some authors say that they are “trying to take over the world” (Negrete, Francavilla, et al., 2020). Although the primary factor that favours this expansion is the transport of adult individuals, regenerating fragments, or possibly by cocoons by humans, it is possible to identify the regions that are the most at risk of invasion by modelling species environmental niches and projecting suitability in geographical space (Barbet-Massin et al., 2018; Gallien et al., 2010). The results we report constitute the first study to use modern modelling approaches to predict the potential range of five hammerhead flatworm species that are already present in several non-native regions. Although large-scale models such as the ones we present here must always be taken with a grain of salt given the many uncertainties that exist with these methods, they can definitely prove useful in order to prioritize management efforts. In this regard, we have identified some potential invasion hotspots that are particularly important to consider, especially those that the species we modelled here have not reached yet (e.g. River Plate basin). These modelling approaches could not be implemented without the existence of massive databases of species’ records; the contribution of enthusiastic naturalists and highly motivated citizens is thus instrumental in management efforts and modelling exercises targeted at non-native species. More precise predictions of the impact of hammerhead flatworms on native ecosystems will now require a deeper evaluation of their interactions with native species, including their prey and potential competitors.

## ACKNOWLEDGEMENTS

We would like to thank the many contributors to the GBIF and iN-aturalist databases, without whom the invasion of terrestrial flat-worms could not be monitored, and who provided the baseline knowledge necessary for conducting this study. We also thank Rebekah D. Wallace, Centre for Invasive Species and Ecosystem Health, University of Georgia, for her collaboration in reviewing iN-aturalist records of the five study species in the USA.

## CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.
## DATA AVAILABILITY STATEMENT

All data used in this study are publicly available in online databases. Occurrence records are available at www.gbif.org (https://doi.org/10.15468/dl.awpn2e) and www.inaturalist.org. Climate variables can be found at https://chelsa-climate.org, and soil variables are available at https://zenodo.org/record/2525664 and https://zenodo.org/record/2525553.

## ORCID

Yoan Fourcade ID https://orcid.org/0000-0003-3820-946X

Leigh Winsor ID https://orcid.org/0000-0002-6679-470X

Jean-Lou Justine ID https://orcid.org/0000-0002-7155-4540

## REFERENCES

Álvarez-Presas, M., Carbayo, F., Rozas, J., & Riutort, M. (2011). Land planarians (Platyhelminthes) as a model organism for fine-scale phylo-geographic studies: understanding patterns of biodiversity in the Brazilian Atlantic Forest hotspot. *Journal of Evolutionary Biology*, *24*, 887-896. https://doi.org/10.1111/j.1420-9101.2010.02220.x

Álvarez-Presas, M., Sánchez-Gracia, A., Carbayo, F., Rozas, J., & Riutort, M. (2014). Insights into the origin and distribution of biodiversity in the Brazilian Atlantic forest hot spot: a statistical phylogeographic study using a low-dispersal organism. *Heredity*, *112*, 656-665. https://doi.org/10.1038/hdy.2014.3

Ascensão, F., & Capinha, C. (2017). Aliens on the Move: Transportation Networks and Non-native Species. In L. Borda-de-Água, R. Barrientos, P. Beja, & H. M. Pereira (Eds.), *Railway ecology* (pp. 65–80). Springer International Publishing.

Atwater, D. Z., Ervine, C., & Barney, J. N. (2018). Climatic niche shifts are common in introduced plants. *Nature Ecology & Evolution*, *2*, 34–43. https://doi.org/10.1038/s41559-017-0396-z

Bahn, V., & McGill, B. J. (2013). Testing the predictive performance of distribution models. *Oikos*, *122*, 321–331. https://doi.org/10.1111/j.1600-0706.2012.00299.x

Banks, N. C., Paini, D. R., Bayliss, K. L., & Hodda, M. (2015). The role of global trade and transport network topology in the human-mediated dispersal of alien species. *Ecology Letters*, *18*, 188–199. https://doi.org/10.1111/ele.12397

Barbet-Massin, M., Rome, Q., Villemant, C., & Courchamp, F. (2018). Can species distribution models really predict the expansion of invasive species? *PLoS One*, *13*, e0193085. https://doi.org/10.1371/journal.al.pone.0193085

Barve, N., Barve, V., Jimenez-Valverde, A., Lira-Noriega, A., Maher, S. P., Peterson, A. T., Soberón, J., & Villalobos, F. (2011). The crucial role of the accessible area in ecological niche modeling and species distribution modeling. *Ecological Modelling*, *222*, 1810–1819. https://doi.org/10.1016/j.ecolmodel.2011.02.011

Boag, B., Yeates, G. W., & Johns, P. M. (1998). Limitations to the distribution and spread of terrestrial flatworms with special reference to the New Zealand flatworm (*Artioposthia triangulata*). *Pedobiologia*, *42*, 495–503.

Boll, P. K., & Leal-Zanchet, A. M. (2016). Preference for different prey allows the coexistence of several land planarians in areas of the Atlantic Forest. *Zoology*, *119*, 162–168. https://doi.org/10.1016/j.zool.2016.04.002

Boll, P. K., Rossi, I., Amaral, S. V., & Leal-Zanchet, A. (2015). A taste for exotic food: Neotropical land planarians feeding on an invasive flatworm. *PeerJ*, *3*, e1307. https://doi.org/10.7717/peerj.1307

Breiman, L. (2001). Random forests. *Machine Learning*, *45*, 5–32. https://doi.org/10.1023/A:1010933404324

Burgess, M. G., Ritchie, J., Shapland, J., & Pielke, R. (2020). IPCC baseline scenarios have over-projected CO₂ emissions and economic

growth. *Environmental Research Letters*, *16*, 14016. https://doi.org/10.1088/1748-9326/abcdd2

Cannon, R. J. C., Baker, R. H. A., Taylor, M. C., & Moore, J. P. (1999). A review of the status of the New Zealand flatworm in the UK. *Annals of Applied Biology*, *135*, 597–614. https://doi.org/10.1111/j.1744-7348.1999.tb00892.x

Čapka, J., & Čejka, T. (2021). First record of Obama nungara in Slovakia (Platyhelminthes: Geoplanidae). *Biodiversity & Environment*, *13*, 41–44.

Carbayo, F., Álvarez-Presas, M., Jones, H. D., & Riutort, M. (2016). The true identity of Obama (Platyhelminthes: Geoplanidae) flatworm spreading across Europe. *Zoological Journal of the Linnean Society*, *177*, 5–28. https://doi.org/10.1111/zoj.12358

Carbayo, F., Leal-Zanchet, A. M., & Vieira, E. M. (2002). Terrestrial flatworm (Platyhelminthes: Tricladida: Terricola) diversity versus man-induced disturbance in an ombrophilous forest in southern Brazil. *Biodiversity and Conservation*, *11*, 1091–1104. https://doi.org/10.1023/A:1015865005604

Cardale, E. (1940). A new species of New Zealand Land-Planarian, *Artioposthia civis*. Notes on anatomy, histology and physiology with ecological data. *Transactions of the Royal Society of New Zealand*, *70*, 277–281.

Chaisiri, K., Dusitsittipon, S., Panitvong, N., Ketboonlue, T., Nuamtanong, S., Thaenkham, U., Morand, S., & Dekumyoy, P. (2019). Distribution of the newly invasive New Guinea flatworm *Platydemus manokwari* (Platyhelminthes: Geoplanidae) in Thailand and its potential role as a paratenic host carrying Angiostrongylus malaysiensis larvae. *Journal of Helminthology*, *93*, 711–719. https://doi.org/10.1017/S0022149X18000834

Chamberlain, S. (2020). spocc: Interface to species occurrence data sources. R package version 1.1.0. https://CRAN.R-project.org/package-spocc

Chen, I. C., Hill, J. K., Ohlemuller, R., Roy, D. B., & Thomas, C. D. (2011). Rapid range shifts of species associated with high levels of climate warming. *Science*, *333*, 1024–1026. https://doi.org/10.1126/science.1206432

Daly, J. J., & Matthews, H. M. (1982). Effect of weight and temperature upon oxygen consumption of the land Planarian Bipalium kewense. *Physiological Zoology*, *55*, 148–154. https://doi.org/10.1086/physzool.55.2.30155850

Darwin, C. (1844). Brief descriptions of several terrestrial Planariæ, and of some remarkable marine species, with an account of their habits. *The Annals and Magazine of Natural History*, *14*, 241–251. https://doi.org/10.1080/037454809495165

Doherty, T. S., Glen, A. S., Nimmo, D. G., Ritchie, E. G., & Dickman, C. R. (2016). Invasive predators and global biodiversity loss. *Proceedings of the National Academy of Sciences of the United States of America*, *113*, 11261–11265. https://doi.org/10.1073/pnas.1602480113

Ducey, P. K., McCormick, M., & Davidson, E. (2007). Natural History observations on Bipalium cf. vagum Jones and Sterrer (Platyhelminthes: Tricladida), a terrestrial broadhead planarian new to North America. *Southeastern Naturalist*, *6*, 449–460. https://doi.org/10.1656/1528-7092(2007)6[449:NHOOBC]2.0.CO;2

Ducey, P. K., & Noce, S. (1998). Successful invasion of New York State by the Terrestrial Flatworm, Bipalium adventitium. *Northeastern Naturalist*, *5*, 199–206. https://doi.org/10.2307/3858619

Duncan, R. P. (1997). The role of competition and introduction effort in the success of Passeriform birds introduced to New Zealand. *The American Naturalist*, *149*, 903–915. https://doi.org/10.1086/286029

Edwards, J. L. (2004). Research and societal benefits of the Global Biodiversity Information Facility. *BioScience*, *54*, 485–486. https://doi.org/10.1641/0006-3568(2004)054[0486:RASBOT]2.0.CO;2

Elith, J., Kearney, M., & Phillips, S. J. (2010). The art of modelling range-shifting species. *Methods in Ecology and Evolution*, *1*, 330–342. https://doi.org/10.1111/j.2041-210X.2010.00036.x
Elith, J., & Leathwick, J. R. (2009). Species distribution models: ecological explanation and prediction across space and time. *Annual Review of Ecology Evolution and Systematics*, *40*, 677-697. https://doi.org/10.1146/annurev.ecolsys.110308.120159

Elith, J., Leathwick, J. R., & Hastie, T. (2008). A working guide to boosted regression trees. *Journal of Animal Ecology*, *77*, 802-813. https://doi.org/10.1111/j.1365-2656.2008.01390.x

Fiore, C., Tull, J. L., Zehner, S., & Ducey, P. K. (2004). Tracking and prediction on earthworms by the invasive terrestrial planarian Bipalium adventitium (Tricladida, Platyhelminthes). *Behavioural Processes*, *67*, 327-334. https://doi.org/10.1016/j.beproc.2004.06.001

Fletcher, J. J., & Hamilton, A. G. (1887). Notes on Australian land planarians, with descriptions of some new species. Part I. *Proceedings of the Linnean Society of New South Wales*, *2*, 349-374. https://doi.org/10.5962/bhl.part.29186

Fourcade, Y. (2021). Fine-tuning niche models matters in invasion ecology. A lesson from the land planarian Obama nungara. *Ecological Modelling*, *457*, 109686. https://doi.org/10.1016/j.ecolomdel.2021.109686

Fourcade, Y., Besnard, A. G., & Secondi, J. (2018). Paintings predict the distribution of species, or the challenge of selecting environmental predictors and evaluation statistics. *Global Ecology & Biogeography*, *27*, 245-256. https://doi.org/10.1111/geb.12684

Fourcade, Y., Engler, J. O., Rödder, D., & Secondi, J. (2014). Mapping species distributions with MAXENT using a geographically biased sample of presence data: a performance assessment of methods for correcting sampling bias. *PLoS One*, *9*, e97122. https://doi.org/10.1371/journal.pone.0097122

Gallien, L., Münkemüller, T., Albert, C. H., Boulangeat, I., & Thuiller, W. (2010). Predicting potential distributions of invasive species: Where to go from here? *Diversity and Distributions*, *16*, 331-342. https://doi.org/10.1111/j.1472-4642.2010.00652.x

GBIF.org (2021). GBIF Occurrence Download. https://doi.org/10.15468/dl.awpn2e

Gerlach, J. (2019). Predation by invasive Platydemus manokwari flatworms: A laboratory study. *Biological Letters*, *54*, 47-60. https://doi.org/10.2478/biolet-2019-0005

Gerlach, J., Barker, G. M., Bick, C. S., Bouchet, P., Brodie, G., Christensen, C. C., Collins, T., Coote, T., Cowie, R. H., Fiedler, G. C., Griffiths, O. L., Florens, F. B. V., Hayes, K. A., Kim, J., Meyer, J.-Y., Meyer, W. M., Richling, I., Slapcinsky, J. D., Winsor, L., & Yeung, N. W. (2021). Negative impacts of invasive predators used as biological control agents against the pest snail Lissachatina fulica: the snail Euglandina 'rosea' and the flatworm Platydemus manokwari. *Biological Invasions*, *23*, 997-1031. https://doi.org/10.1007/s10530-020-02436-w

Hanley, A., & McNeil, J. (1982). The meaning and use of the area under a receiver operating characteristic (ROC) curve. *Radiology*, *143*, 29-36. https://doi.org/10.1148/radiology.143.1.7063747

Hao, T., Elith, J., Lahoz-Monfort, J. J., & Guillera-Arroita, G. (2020). Testing whether ensemble modelling is advantageous for maximising predictive performance of species distribution models. *Ecography*, *43*, 549-558. https://doi.org/10.1111/ecog.04890

Hengl, T. (2018). Soil pH in H2O at 6 standard depths (0, 10, 30, 60, 100 and 200 cm) at 250 m resolution. https://doi.org/10.5281/zenodo.2525664

Hengl, T., & Wheeler, I. (2018). Soil organic carbon content in x 5 g/kg at 6 standard depths (0, 10, 30, 60, 100 and 200 cm) at 250 m resolution. https://doi.org/10.5281/zenodo.2525553

Hijmans, R. J., Phillips, S., Leathwick, J., & Elith, J. (2020). *dismo: Species Distribution Modeling*.

Hu, J., Yang, M., Ye, E. R., Ye, Y., & Niu, Y. (2019). First record of the New Guinea flatworm *Platydemus manokwari* (Platyhelminthes, Geoplanidae) as an alien species in Hong Kong Island, China. *ZooKeys*, *873*, 1-7. https://doi.org/10.3897/zookeys.873.36458

Iturralde, G. G., Allgayer, H., Valiati, V. H., & Leal-Zanchet, A. M. (2021). Seeing the true colours: Three new species of Obama (Platyhelminthes: Continenticola) from remnants of Atlantic Forest in Southern Brazil based on an Integrative Approach. *Invertebrate Systematics*, *35*, 312-331. https://doi.org/10.1071/IS20043

Jiménez-Valverde, A. (2012). Insights into the area under the receiver operating characteristic curve (AUC) as a discrimination measure in species distribution modelling. *Global Ecology and Biogeography*, *21*, 498-507. https://doi.org/10.1111/j.1466-8238.2011.00683.x

Jones, H. D. (2019). Another alien terrestrial planarian in the United Kingdom: Australopacifica atrata (Steel, 1897) (Platyhelminthes: Tricladida: Continenticola). *Zootaxa*, *4604*, 575-587. https://doi.org/10.11646/zootaxa.4604.3.12

Jones, H. D., Mateos, E., Riutort, M., & Álvarez-Presas, M. (2020). The identity of the invasive yellow-striped terrestrial planarian found recently in Europe: Caenoplana variegata (Fletcher & Hamilton, 1888) or Caenoplana bicolor (Graff, 1899)? *Zootaxa*, *4731*, 193-222. https://doi.org/10.11646/zootaxa.4731.2.2

Jones, H. D., Santoro, G., Boag, B., & Neilson, R. (2001). The diversity of earthworms in 200 Scottish fields and the possible effect of New Zealand land flatworms (*Arthurdendyus triangulatus*) on earthworm populations. *Annals of Applied Biology*, *139*, 75-92. https://doi.org/10.1111/j.1744-7348.2001.tb00132.x

Jones, H. D., & Sluys, R. (2016). A new terrestrial planarian species of the genus Marionfyea (Platyhelminthes: Tricladida) found in Europe. *Journal of Natural History*, *50*, 2673-2690. https://doi.org/10.1080/00222933.2016.1208907

Justine, J.-L., Gastineau, R., Gros, P., Ruzzier, E., Charles, L., & Winsor, L. (2022). Hammerhead flatworms (Platyhelminthes, Geoplanidae, Bipaliinae): mitochondrial genomes and description of two new species from France, Italy, and Mayotte. *PeerJ*, *10*, e12725. https://doi.org/10.7717/peerj.12725

Justine, J.-L., Gey, D., Vasseur, J., Thévenot, J., Coulis, M., & Winsor, L. (2021). Presence of the invasive land flatworm Platydemus manokwari (Platyhelminthes, Geoplanidae) in Guadeloupe, Martinique and Saint Martin (French West Indies). *Zootaxa*, *4951*, 381-390. https://doi.org/10.11646/zootaxa.4951.2.11

Justine, J.-L., Lemarcis, T., Gerlach, J., & Winsor, L. (2018). First report of the land planarian Endeavouria septemlineata (Hyman, 1939) (Platyhelminthes, Tricladida, Continenticola, Geoplanidae) in French Polynesia. *Zootaxa*, *4450*, 297-300. https://doi.org/10.11646/zootaxa.4450.2.10

Justine, J.-L., Théry, T., Gey, D., & Winsor, L. (2019). First record of the invasive land flatworm Bipalium adventitium (Platyhelminthes, Geoplanidae) in Canada. *Zootaxa*, *4656*, 591-595. https://doi.org/10.11646/zootaxa.4656.3.13

Justine, J.-L., Winsor, L., Barrière, P., Fanai, C., Gey, D., Han, A. W. K., Quay-Velázquez, G. L., Lee, B.-P.-Y.-H., Lefevre, J.-M., Meyer, J.-Y., Philippart, D., Robinson, D. G., Thévenot, J., & Tsatsia, F. (2015). The invasive land planarian Platydemus manokwari (Platyhelminthes, Geoplanidae): records from six new localities, including the first in the USA. *PeerJ*, *3*, e1037. https://doi.org/10.7717/peerj.1037

Justine, J.-L., Winsor, L., Gey, D., Gros, P., & Thévenot, J. (2014). The invasive New Guinea flatworm Platydemus manokwari in France, the first record for Europe: Time for action is now. *PeerJ*, *2*, e297. https://doi.org/10.7717/peerj.297

Justine, J.-L., Winsor, L., Gey, D., Gros, P., & Thévenot, J. (2018). Giant worms chez moi! Hammerhead flatworms (Platyhelminthes, Geoplanidae, Bipalium spp., Diversibipalium spp.) in metropolitan France and overseas French territories. *PeerJ*, *6*, e4672. https://doi.org/10.7717/peerj.4672

Justine, J.-L., Winsor, L., Gey, D., Gros, P., & Thévenot, J. (2020). Obama chez moi! The invasion of metropolitan France by the land planarian Obama nungara (Platyhelminthes, Geoplanidae). *PeerJ*, *8*, e8385. https://doi.org/10.7717/peerj.8385
Karger, D. N., Conrad, O., Bohner, J., Kawohl, T., Kreft, H., Soria-Azuza, R. W., Zimmermann, N. E., Linder, H. P., & Kessler, M. (2017). Climatologies at high resolution for the earth's land surface areas. *Scientific Data*, *4*, 170122. https://doi.org/10.1038/sdata.2017.122

Kawaguti, S. (1932). On the physiology of land planarians. III. The problems of desiccation. *Memoirs of the Faculty of Science and Agriculture, Taihoku Imperial University, Formosa, Japan*, *7*, 39-55.

Keane, R. M., & Crawley, M. J. (2002). Exotic plant invasions and the enemy release hypothesis. *Trends in Ecology & Evolution*, *17*, 164-170. https://doi.org/10.1016/S0169-5347(02)02499-0

Kramer-Schadt, S., Niedballa, J., Pilgrim, J. D., Schröder, B., Lindenborn, J., Reinfelder, V., Stillfried, M., Heckmann, I., Scharf, A. K., Augeri, D. M., Cheyne, S. M., Hearn, A. J., Ross, J., Macdonald, D. W., Mathai, J., Eaton, J., Marshall, A. J., Semiadi, G., Rustam, R., ... Wilting, A. (2013). The importance of correcting for sampling bias in MaxEnt species distribution models. *Diversity and Distributions*, *19*, 1366-1379. https://doi.org/10.1111/ddi.12096

Liu, C., White, M., & Newell, G. (2013). Selecting thresholds for the prediction of species occurrence with presence-only data. *Journal of Biogeography*, *40*, 778-789. https://doi.org/10.1111/jbi.12058

Liu, C., Wolter, C., Xian, W., & Jeschke, J. M. (2020). Most invasive species largely conserve their climatic niche. *Proceedings of the National Academy of Sciences of the United States of America*, *117*, 23643-23651. https://doi.org/10.1073/pnas.2004289117

Liu, T., Chen, X., Gong, X., Lubbers, I. M., Jiang, Y., Feng, W., Li, X., Whalen, J. K., Bonkowski, M., Griffiths, B. S., Hu, F., & Liu, M. (2019). Earthworms coordinate soil biota to improve multiple ecosystem functions. *Current Biology*, *29*, 3420-3429.e5. https://doi.org/10.1016/j.cub.2019.08.045

Lobo, J. M., Jimenez-Valverde, A., & Real, R. (2008). AUC: a misleading measure of the performance of predictive distribution models. *Global Ecology and Biogeography*, *17*, 145-151. https://doi.org/10.1111/j.1466-8238.2007.00358.x

Lydeard, C., Cowie, R. H., Ponder, W. F., Bogan, A. E., Bouchet, P., Clark, S. A., Cummings, K. S., Frest, T. J., Gargominy, O., Herbert, D. G., Hershler, R., Perez, K. E., Roth, B., Seddon, M., Strong, E. E., & Thompson, F. G. (2004). The Global Decline of Nonmarine Mollusks. *BioScience*, *54*, 321-330. https://doi.org/10.1641/0006-3568(2004)054[0321:TGDONM]2.0.CO;2

Maron, J.L., & Marler, M. (2008). Field-based competitive impacts between invaders and natives at varying resource supply. *Journal of Ecology*, *96*, 1187-1197. https://doi.org/10.1111/j.1365-2745.2008.01440.x

Morffe, J., García, N., Adams, B., & Hasegawa, K. (2016). First record of the land planarian Bipalium kewense Moseley, 1878 (Tricladida: Geoplanidae: Bipaliinae) from Cuba. *Biolnvasions Records*, *5*, 127-132. https://doi.org/10.3391/bir.2016.5.3.01

Mori, E., Magoga, G., Panella, M., Montagna, M., Winsor, L., Justine, J.-L., Menchetti, M., Schifani, E., Melone, B., & Mazza, G. (2022). Discovering the Pandora's box: The invasion of alien flatworms in Italy. *Biological Invasions*, *24*, 205-216. https://doi.org/10.1007/s10530-021-02638-w

Murchie, A. K., & Gordon, A. W. (2013). The impact of the 'New Zealand flatworm', *Arthurdendyus triangulatus*, on earthworm populations in the field. *Biological Invasions*, *15*, 569-586. https://doi.org/10.1007/s10530-012-0309-7

Negrete, L., Do Amaral, S. V., Ribeiro, G. G., Wolmann Gonçalves, J., Valiati, V. H., Damborenea, C., Brusa, F., & Leal-Zanchet, A. M. (2020). Far away, so close! Integrative taxonomy reveals a new genus and species of land flatworm (Platyhelminthes: Geoplanidae) from southern South America. *Zoological Journal of the Linnean Society*, *189*, 722-744. https://doi.org/10.1093/zoolinnean/zlz131

Negrete, L., Francavilla, M. L., Damborenea, C., & Brusa, F. (2020). Trying to take over the world: Potential distribution of Obama nungara (Platyhelminthes: Geoplanidae), the Neotropical land planarian that has reached Europe. *Global Change Biology*, *26*, 4907-4918. https://doi.org/10.1111/gcb.15208

Negrete, L., Gira, R. D., & Brusa, F. (2019). Two new species of land planarians (Platyhelminthes, Tricladida, Geoplanidae) from protected areas in the southern extreme of the Paranaense Rainforest, Argentina. *Zoologischer Anzeiger*, *279*, 38-51. https://doi.org/10.1016/j.jcz.2019.01.002

Ogren, R. E. (1986). The human factor in the spread of an exotic land planarian in Pennsylvania. *Proceedings of the Pennsylvania Academy of Science*, *59*, 117-118.

Petitpierre, B., Kueffer, C., Broennimann, O., Randin, C., Daehler, C., & Guisan, A. (2012). Climatic niche shifts are rare among terrestrial plant invaders. *Science*, *335*, 1344-1348. https://doi.org/10.1126/science.1215933

Phillips, S. J., Anderson, R. P., & Schapire, R. E. (2006). Maximum entropy modeling of species geographic distributions. *Ecological Modelling*, *190*, 231-259. https://doi.org/10.1016/j.ecolmodel.2005.03.026

Richardson, D. M., Pysek, P., Rejmanek, M., Barbour, M. G., Panetta, F. D., & West, C. J. (2000). Naturalization and invasion of alien plants: concepts and definitions. *Diversity and Distributions*, *6*, 93-107. https://doi.org/10.1046/j.1472-4642.2000.00083.x

Rödder, D., & Engler, J. O. (2011). Quantitative metrics of overlaps in Grinnellian niches: Advances and possible drawbacks. *Global Ecology and Biogeography*, *20*, 915-927. https://doi.org/10.1111/j.1466-8238.2011.00659.x

Schirch, P. (1929). Sobre as planarias terrestres do Brasil. *Boletim do Museu Nacional Rio de Janeiro*, *5*, 27-38.

Schoener, T. W. (1968). The Anolis lizards of Bimini: resource partitioning in a complex fauna. *Ecology*, *49*, 704-726. https://doi.org/10.2307/1935534

Schwalm, C. R., Glendon, S., & Duffy, P. B. (2020). RCP8.5 tracks cumulative CO2 emissions. *Proceedings of the National Academy of Sciences of the United States of America*, *117*, 19656-19657. https://doi.org/10.1073/pnas.2007117117

Simberloff, D., Parker, I. M., & Windle, P. N. (2005). Introduced species policy, management, and future research needs. *Frontiers in Ecology and the Environment*, *3*, 12-20. https://doi.org/10.1890/1540-9295(2005)003[0012:ISPMAF]2.0.CO;2

Sluys, R. (1999). Global diversity of land planarians (Platyhelminthes, Tricladida, Terricola): A new indicator-taxon in biodiversity and conservation studies. *Biodiversity and Conservation*, *8*, 1663-1681. https://doi.org/10.1023/A:1008994925673

Sluys, R. (2016). Invasion of the flatworms. *American Scientist*, *104*, 288. https://doi.org/10.1511/2016.122.288

Steel, T. (1897). Australian land planarians: Descriptions of new species and notes on collecting and preserving. *Proceedings of the Linnean Society of New South Wales*, *22*, 104-119. https://doi.org/10.5962/bhl.part.12706

Sugiura, S. (2009). Seasonal fluctuation of invasive flatworm predation pressure on land snails: Implications for the range expansion and impacts of invasive species. *Biological Conservation*, *142*, 3013-3019. https://doi.org/10.1016/j.biocon.2009.07.032

Sugiura, S., & Yamaura, Y. (2009). Potential impacts of the invasive flatworm Platydemus manokwari on arboreal snails. *Biological Invasions*, *11*, 737-742. https://doi.org/10.1007/s10530-008-9287-1

Thomson, A. M., Calvin, K. V., Smith, S. J., Kyle, G. P., Volke, A., Patel, P., Delgado-Arias, S., Bond-Lamberty, B., Wise, M. A., Clarke, L. E., & Edmonds, J. A. (2011). RCP4.5: a pathway for stabilization of radiative forcing by 2100. *Climatic Change*, *109*, 77. https://doi.org/10.1007/s10584-011-0151-4

Thuiller, W. (2013). On the importance of edaphic variables to predict plant species distributions - limits and prospects. *Journal of Vegetation Science*, *24*, 591-592. https://doi.org/10.1111/jvs.12076

Thuiller, W., Lafourcade, B., Engler, R., & Araújo, M. B. (2009). BIOMOD - a platform for ensemble forecasting of species distributions. *Ecography*, *32*, 369-373. https://doi.org/10.1111/j.1600-0587.2008.05742.x

Tylianakis, J. M., Didham, R. K., Bascompte, J., & Wardle, D. A. (2008). Global change and species interactions in
terrestrial ecosystems. *Ecology Letters*, **11**, 1351-1363. https://doi.org/10.1111/j.1461-0248.2008.01250.x

von Graff, L. (1889). Monographie der Turbellarien. II. Tricladida Terricola (Landplanarien). In: *Atlas von Achtundfunfzig Tafeln zur Monographie der Turbellarien. II. Tricladida Terricola (Landplanarien)* (pp. ixv+1-574, Pls. I-LVIII). Verlag von Wilhelm Engelmann.

Wallace, R. D., & Winsor, L. (2020). *Hammerhead worms Fact Sheet*. Centre for Invasive Species and Ecosystem Health University of Georgia. https://resources.ipmcenters.org/view/resource.cfm?rid=23878

Winsor, L. (1977). Terrestrial planarians and nemerteans of the Otway Region. *Proceedings of the Royal Society of Victoria*, **89**, 137-146.

Winsor, L. (1983). A revision of the cosmopolitan land planarian Bipalium kewense Moseley, 1878 (Turbellaria: Tricladida: Terricola). *Zoological Journal of the Linnean Society*, **79**, 61-100. https://doi.org/10.1111/j.1096-3642.1983.tb01161.x

Winsor, L. (1998). The Australian terrestrial flatworm fauna (Tricladida: Terricola). *Pedobiologia*, **42**, 457-463.

Winsor, L. (2003). *Studies on the systematics and biogeography of terrestrial flatworms (Platyhelminthes: Tricladida: Terricola) of the Australian region (PhD)*. James Cook University.

Winsor, L., Johns, P. M., & Barker, G. M. (2004). Terrestrial planarians (Platyhelminthes: Tricladida: Terricola) predaceous on terrestrial gastropods. In G. M. Barker (Ed.), *Natural enemies of terrestrial molluscs* (pp. 227-278). CABI.

Winsor, L., Johns, P., & Yeates, G. (1998). Introduction, and ecological and systematic background, to the Terricola (Tricladida). *Pedobiologia*, **42**, 389-404.

Yackulic, C. B., Chandler, R., Zipkin, E. F., Royle, J. A., Nichols, J. D., Campbell Grant, E. H., & Veran, S. (2013). Presence-only modelling using MAXENT: When can we trust the inferences? *Methods in Ecology and Evolution*, **4**, 236-243. https://doi.org/10.1111/2041-210x.12004

Zurell, D., Franklin, J., König, C., Bouchet, P. J., Dormann, C. F., Elith, J., Fandos, G., Feng, X., Guillera-Arroita, G., Guisan, A., Lahoz-Monfort, J. J., Leitão, P. J., Park, D. S., Peterson, A. T., Rapacciuolo, G., Schmatz, D. R., Schröder, B., Serra-Diaz, J. M., Thuiller, W., ... Merow, C. (2020). A standard protocol for reporting species distribution models. *Ecography*, **43**, 1261-1277. https://doi.org/10.1111/ecog.04960

## BIOSKETCHES

Yoan Fourcade is Associate Professor at Université Paris Est Créteil and in the Institute of Ecology and Environmental Sciences of Paris. He is interested in conservation biology and macroecology, especially with regard to the impact of climate change on biodiversity. Leigh Winsor is Adjunct Zoologist at James Cook University, Queensland, Australia, specializing in the taxonomy, systematics and histology of land flatworms. Jean-Lou Justine is Professor in the National Museum of Natural History in Paris. His interests include parasitology and zoology, especially invasive land flatworms.

Author contribution: Yoan Fourcade: Conceptualization, Formal Analysis, Methodology, Visualization, Writing – Original Draft Preparation, Writing – Review & Editing; Leigh Winsor: Data Curation, Writing – Review & Editing; Jean-Lou Justine: Conceptualization, Data Curation, Writing – Review & Editing

## SUPPORTING INFORMATION

Additional supporting information may be found in the online version of the article at the publisher's website.

**How to cite this article:** Fourcade, Y., Winsor, L., & Justine, J.-L. (2022). Hammerhead worms everywhere? Modelling the invasion of bipalini flatworms in a changing climate. *Diversity and Distributions*, **00**, 1-15. https://doi.org/10.1111/ddi.13489
# Appendix S1

**Methods S1:** ODMAP protocol for reporting details about the SDM procedure

## OVERVIEW

**Title of study:** Hammerhead worms everywhere? Modelling the invasion of bipaliin flatworms in a changing climate

**Authorship:** **Yoan Fourcade**
Univ Paris Est Creteil, Sorbonne Université, CNRS, IRD, INRAE, Institut d'écologie et des sciences de l'environnement, IEES, F-94010 Créteil, France

**Leigh Winsor**
College of Science and Engineering, James Cook University, Townsville, QLD, Australia

**Jean-Lou Justine**
ISYEB, Institut de Systématique Évolution Biodiversité, UMR7205 MNHN, CNRS, EPHE, UPMC, Université des Antilles, Muséum National d'Histoire Naturelle, 43 Rue Cuvier, 75005 Paris, France

**Contact:** yoan.fourcade@u-pec.fr

## MODEL OBJECTIVE

**Model objective:** Forecast and transfer

**Target output:** Areas of potential invasion for invasive hammerhead flatworms

**Focal Taxon:** *Bipalium adventitium*, *B. kewense*, *Diversibipalium multilineatum*, *B. pennsylvanicum* and *B. vagum*

## LOCATION

**Location:** Global

## SCALE OF ANALYSIS

**Spatial extent:** -180, 180, -90, 90 (xmin, xmax, ymin, ymax)

**Spatial resolution:** 2.5 arc-min

**Temporal extent:** present (1979-2013), future (2041-2060 and 2061-2080)

**Boundary:** None

## BIODIVERSITY DATA

**Observation type:** Opportunistic records

**Response data type:** Presence vs. pseudo-absence / background

## PREDICTORS

**Predictor types:** Climatic, edaphic

## ASSUMPTIONS

**Model assumptions:** We assumed that at the scale of the study, hammerhead flatworms' distributions are mostly controlled by climate, although we also include soil predictors. We did not account for dispersal in the models because we aimed to
predict global suitability, i.e. regions of potential invasion in the case of accidental human transport.

## ALGORITHMS

**Modelling techniques:** Generalised linear model (GLM), boosted regression trees (BRT), maximum entropy (MaxEnt) and random forest (RF)

**Model averaging:** Mean of all algorithms and data subsets weighted by their AUC, discarding models with AUC < 0.5

## WORKFLOW

**Model workflow:**
(1) Occurrence data were obtained from GBIF and iNaturalist databases, (2) data were checked and spatially filtered to remove errors and sampling bias, (3) ensemble models were trained using four algorithms, (4) models were projected in the present and in two future scenarios of climate change and two time periods, merging predictions from five global circulation models

## SOFTWARE

**Software:** biomod2 package in R

**Code availability:** Can be provided upon request

**Data availability:**
Climate variables: https://chelsa-climate.org/
Soil pH: https://zenodo.org/record/2525664
Soil organic carbon content: https://zenodo.org/record/2525553
Earthworm data: iNaturalist and GBIF

# DATA

## BIODIVERSITY DATA

**Taxon names:**
5 species of hammerhead flatworms:
- *Bipalium adventitium*
- *Bipalium kewense*
- *Diversibipalium multilineatum*
- *Bipalium pennsylvanicum*
- *Bipalium vagum*

**Ecological level:** Species

**Data sources:** Data from iNaturalist and GBIF

**Sampling design:** No sampling design; opportunistic records

**Sample size:** 159 to 4553 occurrences per species in raw data; 43 to 697 after spatial thinning (mean = 254)

**Clipping:** None

**Cleaning:**
More than 6,000 photographs were examined and edited, removing inconsistencies and misidentified species.
Spatial thinning of occurrences, keeping only one coordinate per grid cell in a raster of 50 × 50 km resolution, prioritizing iNaturalist data.

**Absence data:** None

**Background data:** 10000 background points within 500-km buffers around occurrences

## Data partitioning

**Training data:** Occurrences divided into four geographical subsets corresponding to the four corners of occurrences' coordinates separated by latitude and longitude midpoints (i.e. spatial blocks)
**Validation data:** Occurrences divided into four geographical subsets corresponding to the four corners of occurrences' coordinates separated by latitude and longitude midpoints (i.e. spatial blocks)

## Predictor variables

**Predictor variables:** Maximum temperature of the warmest month (BIO5), minimum temperature of the coldest month (BIO6), precipitation of the wettest month (BIO13), and precipitation of the driest month (BIO14). We removed annual precipitation (BIO12) and mean annual temperature (BIO1) because of high correlation with other predictors.

Soil pH in water and soil organic carbon at the surface (0 cm).

**Data sources:** CHELSA (climate): https://chelsa-climate.org/,

OpenLandMap (soil): https://zenodo.org/record/2525664
https://zenodo.org/record/2525553

**Spatial extent:** -180, 180, -90, 90 (xmin, xmax, ymin, ymax)

**Spatial resolution:** 2.5 arc-min

**Coordinate reference system:** Initially EPSG 4326, model predictions were projected to EPSG 54009

**Temporal extent:** Climate: 1979-2013
Soil: 1950-2017

**Data processing:** Environmental layers were obtained as raster files with a resolution of 30 arc-sec, then were aggregated to 5 arc-min.

## Transfer data

**Data sources:** CHELSA (climate): https://chelsa-climate.org/

**Spatial extent:** -180, 180, -90, 90 (xmin, xmax, ymin, ymax)

**Spatial resolution:** 2.5 arc-min

**Temporal extent:** Climate: 2041-2060, 2061-2080

**Models and scenarios:**
Global circulation models: ACCESS1.0, CESM1-BGC, CanESM2, HadGEM2-AO, MIROC-ESM
Climate change scenarios: RCP 4.5 and RCP 8.5

**Data processing:** Environmental layers were obtained as raster files with a resolution of 30 arc-sec, then were aggregated to 5 arc-min.

**Quantification of Novelty:** Multivariate environmental similarity surfaces (MESS)

# MODEL

## Multicollinearity

**Multicollinearity:** We removed annual precipitation (BIO12) and mean annual temperature (BIO1) from the initial set of bioclimatic variables because of high correlation with other predictors.

## Model settings

**GLM:** Variable selection based on AIC; quadratic relationships and one-way interactions

**MaxEnt:** Features: linear, quadratic and hinge; regularization parameter: 3
BRT: Number of trees: 2500; interaction depth: 7

RF: Number of trees: 5000

## Threshold selection

Threshold selection: Based on the threshold that maximizes the sum of sensitivity and specificity

## ASSESSMENT

### Performance statistics

Performance on validation data: Area under the ROC curve (AUC) calculated from four geographical subsets corresponding to the four corners of occurrences' coordinates separated by latitude and longitude midpoints

### Plausibility check

Response shapes: Extracted for each variable as the predicted suitability (according to the ensemble model) across the range of available values, keeping the other predictors at their mean values

Expert judgement: Visualization of output binary maps

## PREDICTION

### Prediction output

Prediction unit: Continuous suitability maps
Binary presence/absence maps

Post-processing: We stacked presence/absence predictions to map potential species richness of invasive hammerhead flatworms. Global potential range area was calculated in current and future climatic conditions. We estimated mean latitudinal shift between current and future predicted ranges. Overlap between species' distributions was assessed using Schoener's *D* index.

### Uncertainty quantification

Scenario uncertainty: Projection into two scenarios and two time periods; using five global circulation models

Novel environments: Multivariate environmental similarity surfaces (MESS)
<table><caption>**Table S1:** Area under the ROC curve (AUC) obtained for each species and modelling algorithm. Since all models were conducted four times, corresponding with four spatially independent validation datasets, we show the mean AUC as well as the minimum and maximum AUC across these replicates.</caption><thead><tr><th>Species</th><th>Boosted reg. Trees</th><th>GLM</th><th>MaxEnt</th><th>Random forest</th><th>Mean by species</th></tr></thead><tbody><tr><td>Bipalium adventitium</td><td>0.81 [0.70-0.87]</td><td>0.76 [0.71-0.81]</td><td>0.78 [0.70-0.87]</td><td>0.73 [0.68-0.78]</td><td>0.77</td></tr><tr><td>Bipalium kewense</td><td>0.79 [0.71-0.88]</td><td>0.79 [0.72-0.88]</td><td>0.75 [0.65-0.86]</td><td>0.77 [0.73-0.84]</td><td>0.78</td></tr><tr><td>Diversibipalium multilineatum</td><td>0.66 [0.55-0.77]</td><td>0.59 [0.58-0.60]</td><td>0.60 [0.55-0.68]</td><td>0.51 [0.50-0.54]</td><td>0.59</td></tr><tr><td>Bipalium pennsylvanicum</td><td>0.78 [0.55-0.96]</td><td>0.73 [0.60-0.84]</td><td>0.78 [0.64-0.85]</td><td>0.53 [0.20-0.79]</td><td>0.70</td></tr><tr><td>Bipalium vagum</td><td>0.74 [0.57-0.90]</td><td>0.73 [0.61-0.86]</td><td>0.74 [0.56-0.89]</td><td>0.69 [0.54-0.83]</td><td>0.72</td></tr><tr><td>Mean by algorithm</td><td>0.76</td><td>0.72</td><td>0.73</td><td>0.64</td><td></td></tr></tbody></table>
<table><caption>**Table S2:** Potential species range size obtained after converting continuous suitability predictions into binary presence/absence maps using the threshold that maximises sensitivity and specificity, displayed in million km². Blue bars represent range areas relative to the minimum and maximum values across all models and species.</caption><thead><tr><th rowspan="3">Species</th><th rowspan="3">Current</th><th colspan="5">Future</th></tr><tr><th rowspan="2">GCM</th><th colspan="2">RCP 4.5</th><th colspan="2">RCP 8.5</th></tr><tr><th>2041 - 2060</th><th>2061 - 2080</th><th>2041 - 2060</th><th>2061 - 2080</th></tr></thead><tbody><tr><td rowspan="5">Bipalium adventitium</td><td>13.22</td><td>ACCESS1</td><td>7.05</td><td>7.71</td><td>6.40</td><td>7.14</td></tr><tr><td></td><td>CanESM2</td><td>8.48</td><td>8.17</td><td>8.96</td><td>9.33</td></tr><tr><td></td><td>CESM1</td><td>8.38</td><td>7.29</td><td>8.15</td><td>9.16</td></tr><tr><td></td><td>HadGEM</td><td>6.99</td><td>6.42</td><td>6.15</td><td>7.06</td></tr><tr><td></td><td>MIROC</td><td>9.27</td><td>8.49</td><td>8.78</td><td>11.13</td></tr><tr><td rowspan="5">Bipalium kewense</td><td>7.68</td><td>ACCESS1</td><td>7.68</td><td>7.84</td><td>8.16</td><td>8.95</td></tr><tr><td></td><td>CanESM2</td><td>7.33</td><td>7.45</td><td>7.67</td><td>9.61</td></tr><tr><td></td><td>CESM1</td><td>7.52</td><td>7.45</td><td>7.54</td><td>7.66</td></tr><tr><td></td><td>HadGEM</td><td>7.82</td><td>8.13</td><td>10.10</td><td>9.05</td></tr><tr><td></td><td>MIROC</td><td>7.93</td><td>10.28</td><td>8.22</td><td>9.53</td></tr><tr><td rowspan="5">Diversibipalium multilineatum</td><td>1.10</td><td>ACCESS1</td><td>0.58</td><td>0.66</td><td>0.60</td><td>0.62</td></tr><tr><td></td><td>CanESM2</td><td>0.65</td><td>0.63</td><td>0.58</td><td>0.63</td></tr><tr><td></td><td>CESM1</td><td>0.82</td><td>0.83</td><td>0.88</td><td>0.49</td></tr><tr><td></td><td>HadGEM</td><td>0.68</td><td>0.51</td><td>0.73</td><td>0.56</td></tr><tr><td></td><td>MIROC</td><td>0.76</td><td>0.77</td><td>0.85</td><td>0.76</td></tr><tr><td rowspan="5">Bipalium pennsylvanicum</td><td>0.16</td><td>ACCESS1</td><td>0.19</td><td>0.18</td><td>0.14</td><td>0.11</td></tr><tr><td></td><td>CanESM2</td><td>0.13</td><td>0.09</td><td>0.26</td><td>0.05</td></tr><tr><td></td><td>CESM1</td><td>0.23</td><td>0.30</td><td>0.31</td><td>0.23</td></tr><tr><td></td><td>HadGEM</td><td>0.25</td><td>0.12</td><td>0.20</td><td>0.13</td></tr><tr><td></td><td>MIROC</td><td>0.19</td><td>0.11</td><td>0.21</td><td>0.12</td></tr><tr><td rowspan="5">Bipalium vagum</td><td>7.47</td><td>ACCESS1</td><td>8.87</td><td>10.04</td><td>10.25</td><td>11.78</td></tr><tr><td></td><td>CanESM2</td><td>9.02</td><td>9.28</td><td>10.03</td><td>11.73</td></tr><tr><td></td><td>CESM1</td><td>9.30</td><td>9.00</td><td>9.22</td><td>9.52</td></tr><tr><td></td><td>HadGEM</td><td>9.06</td><td>9.70</td><td>9.43</td><td>11.28</td></tr><tr><td></td><td>MIROC</td><td>9.34</td><td>9.73</td><td>9.58</td><td>11.62</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Figure S1:** Overlap between modelled distributions in current climate, estimated using Schoener's *D*.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure S2:** Overlap between modelled distributions in future climate. Overlap is estimated using the mean Schoener's *D* value across five global circulation models.</figcaption></figure>
**Figure S3:** Latitudinal range shifts predicted for each species in each climate change scenario. Values shown are the mean (+/- 95% confidence interval) differences in the absolute latitude of suitable grid cells (from binary maps) between each future projection and the predicted current distribution. Negative values reflect shifts towards the poles.

<figure><img src="image_10.png" /></figure>
2050 (RCP 4.5)

<figure><img src="image_11.png" /><figcaption>Bipalium kewense</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Diversibipalium multilineatum</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Bipalium pennsylvanicum</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Bipalium vagum</figcaption></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

2070 (RCP 4.5)

<figure><img src="image_18.png" /><figcaption>Bipalium kewense</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Diversibipalium multilineatum</figcaption></figure>

<figure><img src="image_20.png" /></figure>

Species
richness

<figure><img src="image_21.png" /><figcaption>Bipalium pennsylvanicum</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Bipalium vagum</figcaption></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>
<figure><img src="image_25.png" /></figure>

**Appendix S3:** Binary predictions of global environmental suitability, in current and future climate. Continuous suitability maps were converted into binary predictions using the threshold that maximises sensitivity and specificity. For future predictions, we considered areas suitable when at least three GCM agreed. The sum of these binary maps is presented in Figures 3, 4 and Appendix S2 as an estimate of potential species richness.
<figure><img src="image_26.png" /></figure>

**Appendix S4:** Multivariate Environmental Similarity Surfaces (MESS) summarised into similar environments (positive values, in grey) and novel environments (negative values, in dark red), showing where models extrapolated into environmental conditions that were not present in the training data. Future MESS maps are based on the mean MESS values across all GCMs.