Predicting terrestrial dispersal corridors of the invasive African clawed frog Xenopus laevis in Portugal

Invasive species, such as the mainly aquatic African clawed frog Xenopus laevis, are a main threat to global biodiversity. The identification of dispersal corridors is necessary to restrict further expansion of these species and help to elaborate management plans for their control and eradication. Here we use remote sensing derived resistance surfaces, based on the normalised difference vegetation index (NDVI) and the normalised difference water index (NDWI) accounting for behavioural and physiological dispersal limitations of the species, in combination with elevation layers, to determine fine scale dispersal patterns of invasive populations of X. laevis in Portugal, where the frog had established populations in two rivers. We reconstruct past dispersal routes between these two invaded rivers and highlight high risk areas for future expansion. Our models suggest terrestrial dispersal corridors that connect both invaded rivers and identify artificial water bodies as stepping stones for overland movement of X. laevis. Additionally, we found several potential stepping stones into novel areas and provide concrete information for invasive

introduction Worldwide, invasive species are a main threat to biodiversity (e.g. Fritts and Rodda 1998;Rosenzweig 2001;Chornesky and Randall 2003;Davis 2003;Doherty et al. 2016), having also played a role in global amphibian decline (e.g. Gibbons et al. 2000;Chornesky and Randall 2003;Falaschi et al. 2020). Some invasive amphibians are known to have disastrous impacts on native ecological communities and to alter sensitive ecological relationships through competition for resources, predation or spread of infectious diseases (Kraus 2009;Kraus 2015 and references therein; Measey et al. 2016 and references therein; Kumschick et al. 2017).
The African clawed frog (Xenopus laevis) is native to southern Africa and has been moved worldwide as a model organism for laboratory research (Measey et al. 2012;van Sittert and Measey 2016). This species has unique physiological and demographical traits, including a tolerance to saltwater and eutrophic conditions and behavioural adaptations, such as terrestrial migration or the ability to burrow into substrate to persist in drought and extreme temperature events. This trait combination confers it with an enormous invasive potential (for a thorough review, see Measey et al. 2012;Sousa et al. 2018;Scalera et al. 2019). Recently, it was ranked second amongst all invasive amphibian species considering its environmental and socio-economic impacts .
To date, X. laevis has established invasive populations in numerous countries across four continents due to both deliberate and accidental release of laboratory animals and the pet trade (Measey et al. 2012). In Europe, it successfully invaded lotic and lentic freshwater habitats (Moreira et al. 2017) and established populations are known in the U.K., France, Sicily (Italy) and Portugal (see Measey et al. 2012 and references therein). Comparisons of mitochondrial DNA suggest that all the Portuguese frogs resulted from presumably a single introduction event and descend from animals exported from the South-western Cape clade, in Mediterranean from South Africa . Correlative SDM approaches revealed the Mediterranean region of Portugal as climatically highly suitable for X. laevis (Measey et al. 2012;Ihlow et al. 2016). However, while the populations in Sicily and France are spreading fast (Faraone et al. 2008;Louppe et al. 2017), the expansion of the Portuguese populations was comparatively slow -approximately 30 years after the introduction event, the species was still confined to a 30 km 2 region (Sousa et al. 2018).
Dispersal is essential for successful spread of an invasive species (cf. Blackburn et al. 2011). Animal dispersal occurs at different life stages and is triggered to evade competition, to acquire resources, to reduce mortality or for reproduction (Bowler and Benton 2005 and references therein;Van Dyck and Baguette 2005). Amphibians are often referred to as poor dispersers, but some species disperse over considerable distances of more than 10 km, for example, from the terrestrial habitats to spatially disjunct breeding sites for migration (e.g. Avise 2000; Smith and Green 2005 and references therein). Traditionally, X. laevis was thought to be strictly aquatic with all life stages inhabiting the same aquatic environment, constraining the invasive potential of the species to connected aquatic habitats. However, there is now sufficient evidence that the species migrates overland amongst various types of lentic and lotic water bodies, which significantly expanded our view of the species' dispersal and, therefore, invasive potential (see Measey et al. 2016 for recent review). Terrestrial movement seems to be constrained to a fraction of the population (21-36% [Measey and Tinsley 1998]), involves adults of both sexes (De Villiers and Measey 2017), occurs mostly nocturnally, involves Euclidean distances up to 2.36 km (De Villiers 2015) and a maximum velocity of 0.2 km/h (Measey and Tinsley 1998). It is known that a drying habitat or the reduction of resources due to high numbers of conspecifics can lead to mass migration events, although other potential factors that trigger terrestrial dispersal remain unknown .
To predict dispersal pathways and, therefore, be able to block or hamper further expansion of an invasive species, it is possible to build resistance surfaces that reflect different costs for a species to move through the landscape using vegetation cover, elevation, slope or other landscape features (Landguth et al. 2012). As a multitude of paths may exist between two points, either multiple low-cost paths or smoothed output paths using a probability-density function can be considered (Cushman et al. 2009;Pinto and Keitt 2009). The latter approach allows a variety of smoothing functions (Gaussian, Epanechnikov, uniform, triangle, biweight, triweight and cosine function) referred to as kernel density estimations (Li and Racine 2007).
In the present study, we used occurrence records of X. laevis in Portugal and fine scale remote sensing data to build landscape resistance kernels that predict the influence of landscape structure on the dispersal dynamics of this invasive frog. Landscape resistance is subsequently used to identify past dispersal routes and to highlight areas at risk of future invasion by X. laevis. This study provides insights into the role of landscape configuration on dispersal patterns and provides a tool for future management of this species, as well as of others with similar dispersal patterns.

Study area
West Portugal is characterised by a Mediterranean-type climate (Rubel and Kottek 2010). The first record of X. laevis in western Portugal occurred in Laje River, which runs through a densely-urbanised part of Oeiras County ca. 20 km west of Lisbon, in 2006 (Sousa et al. 2018) (Fig. 1A). However, this introduction likely occurred much earlier and is probably the result of accidental escape from nearby research laboratories after the strong winter floods of 1979/80 (Rebelo et al. 2010;Sousa et al. 2018). Due to its cryptic lifestyle combined with a lack of interest in the wildlife of urban rivers, the species established and spread along the river, undetected for more than 25 years (Sousa et al. 2018). The frog then spread from Laje River into a close parallel-flowing second river (Barcarena), where it was found in 2008 (Rebelo et al. 2010). The maximum invaded area by the frog on the river sections occurred along 5.86 km in the main Figure 1. Study area A overview of study area: Elevation and important locations (i.e. the two invaded rivers and other localities that could be threatened by invasion). Features as presence and absence points (streams and ponds), water bodies and the site of introduction are highlighted B areas of low and high risk of invasion: Landscape resistance and the connectivity (including all presences as starting points) of the study area. Features as presence and absence points (streams and ponds) and water bodies are highlighted. Two areas of low but possible risk of invasion are surrounded by black circles C reconstruction of past dispersal from Laje into Barcarena River: Landscape resistance and the connectivity (including only presences from Laje River as starting points) of the study area. Features as presence and absence points (streams and ponds), and water bodies are highlighted. The golf-course ponds, which were used as stepping stones into Barcarena River, are surrounded by a red circle. stream of Laje plus one of its tributaries and 6.39 km in the main stream of Barcarena plus three tributaries (Moreira et al. 2017). The invaded river basins are roughly 3 km apart, but the headwaters of two small tributaries of Laje River are nearby (ca.1 km) to the headwaters of a tributary of Barcarena River (Fig. 1A). Both rivers are permanent, approximately 10 m wide in most stretches and 2 m deep in summer in the deepest stream pools (Moreira et al. 2017). The tributaries inhabited by X. laevis are also permanent, ca. 1 m wide and up to 1.5 m deep (Moreira et al. 2017). Both rivers flow into the Tagus estuary and, although X. laevis can tolerate a moderate salinity, it seems unlikely that it has used this path to cross between both river basins (Sousa et al. 2018). The area has just a few still and artificial water bodies, most of them 20-110 m apart from the streams (Moreira et al. 2017). Three large golf-course ponds are located between the headwaters of two tributaries flowing in opposite directions, one to each river. The invasive populations seem constrained to the two rivers, some of their tributaries, the three golf-course ponds and some man-made ponds, covering a total area of approximately 30 km 2 .

Landscape resistance
We calculated fine scale resistance kernels to determine connectivity and predict potential overland dispersal for the invasive Portuguese population. We used literature-based GPS data of confirmed presences (Rebelo et al. 2010), updated them to a total of 201 locations and added 19 confirmed absences along the streams, as well as in isolated ponds, according to our own field research. To define the study area, we chose a circular buffer of 5 km around the GPS-points, which is about twice the maximum documented terrestrial dispersal distance during a dispersal event of the frog (Measey 2016).

Remote sensing derived resistance surfaces
We obtained high resolution multispectral satellite imagery containing the invaded Portuguese distribution range (625 km² × 4 title IDs = 2500 km 2 ) as A3 products of the RapidEye satellite (Blackbridge 2014). We used the satellite images with title-IDs 2956913, 2956914, 2957013 and 2957014 and card IDs 26196070, 26196539, 26196860, 26195477, 26196831, 26196078, 26196867, 26196894 and 26196746. The dataset contains four orthorectified raster images, each with five remote sensing channels (blue: 440-510 nm, green: 520-590 nm, red: 630-685 nm, red edge: 690-730 nm, NIR: 760-850 nm wavelengths). Each raster image covers an area of 625 km². The spatial resolution is 5 m grid cell size and spatial accuracy is 10 m (Blackbridge 2014). The subsequent corrections were applied to the five raw remote sensing channels: 'top of atmosphere correction' (TOA), 'cloud cover correction' and 'histogram correction' using the packages LANDSAT (Goslee 2011), RASTER (Hijmans 2015) and LMODEL2 (Legendre 2014) for R according to the product specifications. Using RS TOOLBOX (Leutner and Horning 2016) and the above-mentioned packages for R, the raster images for each study area were mosaicked before the 'Normalised Difference Vegetation Index' (NDVI) as a measure for vegetation cover and the 'Normalised Difference Water Index' (NDWI), showing differences in the water content of vegetation, were computed.
Based on remote sensing variables using a threshold-based water detection method, the larger still and flowing freshwater bodies within the study area were detected (Klemenjak et al. 2012;Tetteh and Schönert 2015). To improve the algorithms' capacity to detect water, the precise locations of all verified water bodies were identified manually in Google Earth and used to train a bioclim model using the DISMO package for R (Hijmans, Phillips et al. 2013;R Core Team 2017). The output of the bioclim model was again verified by hand. Small water bodies that had not been detected by this measure were georeferenced by hand using Google Earth. We used only one randomly selected occurrence record per 50 × 50 m grid cell because computation time increases exponentially with the number of species occurrences and the size of the study area. We collected seven presences comprising two invaded river sections (Laje 5.6 km and Barcarena 6.3 km) and identified 612,536 locations potentially adequate for invasion.
We calculated resistance surfaces by combining NDVI and NDWI scores giving higher priority to vegetation cover, but acknowledging that humid areas may be preferred by the frogs (i.e. NDVI + NDWI / 10). We applied an inverse monomolecular transformation using relevant functions of the RESISTANCEGA package for R (Peterman 2014, Peterman et al. 2014) to account for the higher expected permeability of areas covered by humid vegetation. The equation of the inverse monomolecular transformation is y = -r(1 -exp -bx ), with r = resistance surface, which is controlled by shape (x) and magnitude (b) parameters that are varied during optimisation (see Peterman 2018 for more details). This transformation appreciates that the permeability of vegetated areas may not shift too much with declining vegetation cover. However, resistance increases exponentially in more open landscapes, especially in the Mediterranean climate, characterised by hot and dry summers. As our resistance surface is based on vegetation and humidity indices, we can differentiate various states of vegetation, but less so for bare soil. Bare soil usually falls within an NDVI range between 0.1-0.2 and plants will always have positive values between 0.2 and 1. Therefore, an exponential function counterbalances the differences in the discrimination ability of our indices. The resistance surfaces were scaled to range from 0-1. Subsequently, a threshold was determined to detect man-made surfaces, such as roads and buildings by comparing all scores between 0.70-0.85 in steps of 0.01 to areal pictures of the same area and reclassifying all scores above the best-matching threshold to a resistance of 50. This threshold does not differentiate between roads, which are semi-permeable and buildings which represent absolute barriers. A score of 50 allows the frogs to cross a maximum of up to 6 grid cells of 5 m (ca. 30 m) man-made surfaces, making roads semi-permeable, but areas with a high density of man-made structures become impermeable.

Elevation data
Laboratory trials, using X. laevis individuals from Portugal, were used to quantify the effect of slope on dispersal. These trials showed an increasing difficulty to overcome slopes, with 60 degrees as the upper limit. An elevation layer with a spatial resolution of 30 m derived from the 'Advanced Spaceborne Thermal Emission and Reflection Radiometer' 'Global Digital Elevation Model' (ASTER GDEM) was obtained from the online database of the NASA Land Processes Distributed Active Archive Centre (LP DAAC) of the USGS/Earth resources observation and science (EROS) centre (https:// lpdaac.usgs.gov). We re-sampled it to a resolution of 5 m using bidirectional interpolation, available in the RASTER package (Hijmans 2015) for R.

Resistance kernels
The remote sensing derived resistance surfaces, in combination with the elevation data, were used to calculate resistance kernels that quantify permeability of the landscape after Compton et al. (2007), using the UNICOR package for python (Landguth et al. 2012; Fig. 2). The resistance kernel approach combines a kernel density estimator with a directional least-cost matrix to produce a multidirectional probability distribution representing variability in habitat quality (Compton et al. 2007). This measure considers land use and elevation derived from remote sensing data and equals higher permeability with higher connectivity as it suggests a higher probability of a dispersing frog to arrive at the water bodies (Compton et al. 2007).
Based on the laboratory trials, the slope function was defined so that an upward slope of 60 degrees is the maximum, while downward slopes were considered as generally permeable (upward slope function as determined, based on trials: y = 3.1051 e 0.038x , scaled to 0-1; downward resistance = 0; settings UNICOR: Type_Direction = Hiking; 6;-3).
Based on capture-mark-recapture data from South Africa (De Villiers 2015), we determined the maximum cumulative resistance, which is observed in the field using least cost paths, calculated with the same set of remote sensing derived resistance plus elevation data. For this, we used the satellite images with the title-IDs 3423406, 3423407, 3423307 and 3423207 and Card-IDs 26195795, 26195801, 26195803 and 26196899. The highest cumulative resistance detected was used to parameterise the dispersal models for the study area (maximum cumulative resistance = 308). Model accuracy was evaluated by extracting the values of the resistance kernel for the used presence records. Further, true absence records were used in the same way to check if they were in-or outside of the predicted kernel area.

results
UNICOR outputs show the cumulative density of optimal paths buffered by the kernel density estimation (Fig. 1C). As for model accuracy, the 201 presence records had a mean cumulative density of 51.21 ± 20.02 (range: 2.19-80.32), while, for the 19 absence records, this was 19.21 ± 28.13 (range: 0-70.52). Eleven of the 19 absences (57.9%) were outside of the kernels' range. The resistance kernels also located functionally well-connected water bodies and complexes in close vicinity to the existing populations ( Fig. 1B -invaded ponds). The known X. laevis populations, occupying the two rivers, are evidently constrained by landscape resistance and high permeability was attributed only to the valley bottoms around the river beds of Laje and Barcarena. Importantly, our results explain the current distribution of the species, including its absence from nearby streams and locate the probable contact route between the two invaded basins, supporting the hypothesis of a natural colonisation of Barcarena by overland dispersal. In fact, areas of low (but still possible) permeability connect the two valleys at two locations, but the isolated animals found upstream of Barcarena seem to have no connectivity with the main downstream population (Fig. 1B).

Discussion
With this work, it was possible to reconstruct the most probable past dispersal routes, terrestrial corridors for overland dispersal and water bodies that function as stepping stones, fostering the X. laevis invasion. Additionally, we found potential stepping stones into novel areas, now considered of high invasion risk.

Pace of invasion
Despite its dispersal abilities, which include terrestrial movements up to 2.36 km (De Villiers 2015), an ascertained maximum velocity of 0.2 km/h (Measey and Tinsley 1998) and the apparently ideal climatic conditions (Measey et al. 2012;Ihlow et al. 2016), the X. laevis invasion in Portugal was quite slow. This is probably the result of a mostly aquatic (lotic) dispersal route, rather than the terrestrial overland dispersal documented on the species' original range . Dispersal along flowing waters -and, in this case, opposed to the water flow in the first river to be colonised (Fig. 1A) -is affected by a combination of the species' dispersal ability, the location of the introduction site, the hydrological regime and landscape resistance.
The landscape of the Laje and Barcarena basins is hostile to a semi-terrestrial frog (see below). In fact, only a few isolated ponds were colonised (Fig. 1C); several ponds at 50 to 80 m from the streams were not reached by X. laevis. A few triggers of X. laevis terrestrial movement have been identified (De Villiers and Measey 2017). One of the most important triggers seems to be population density, which is here relatively low (Moreira et al. 2017), probably due to the low habitat quality -heavily polluted urban rivers. This low abundance also explains why the species went undetected for more than 25 years (Sousa et al. 2018). Although X. laevis can live, disperse through and successfully reproduce in flowing waters (Lobos and Jaksic 2005;Courant et al. 2017;Moreira et al. 2017), these seem not to be ideal for the species. Lotic habitats have been mostly identified as pathways for dispersal, while breeding is commonly referred to take place in lentic water bodies, like pools or ponds (Fouquet and Measey 2006;Faraone et al. 2008;Measey 2016). In Portugal and probably due to the poor habitat quality of the two streams and/or presence of predatory fish, the number of metamorphs is much lower in lotic than in lentic environments; metamorph size is also smaller, whereby reaching the reproductive size takes longer (Moreira et al. 2017).
Some features of the Mediterranean climate may have also contributed to the slow dispersal. The annual period where terrestrial dispersal could take place is not certain, as the mostly dry and hot summers seem too risky for terrestrial movement. The mild winters could be very suitable for dispersal overland, because these Mediterranean streams are typically subject to high water level variability; the rainy winters regularly cause river floods (Boix et al. 2010), spreading the species along the river valley. However, the site of initial introduction was located downstream, close to the mouth of the River Laje (Fig. 1A), meaning that, until the colonisation of Barcarena, only upstream dispersal was possible. Events like the 1979/80 flood, when the species escaped in Laje River, could boost its dispersal, but hardly upstream.

Distribution
We found that the modelled landscape connectivity correlates well with the distribution of this frog. In the areas of high connectivity along the river-beds of Laje and Barcarena rivers, the species' dispersal is hampered by > 22° slopes and > 60° slopes seem to be nearly unconquerable. Landscape connectivity along large parts of the river sections is further constrained by cement walls instead of natural riverbank. Further, the rivers have several physical barriers like waterfalls, hampering the connectivity amongst populations and reducing landscape permeability (Sousa et al. 2018). In fact, the lack of continuity upstream of the Barcarena basin (Fig. 1C) results from a ~250 m long tunnel through which the stream flows beneath two highways (detected as roads by remote sensing). Frogs were recently (June 2020) found just upstream and downstream of the tunnel and are very probably able to pass through the tunnel.
Away from the riverbeds, connectivity decreases very quickly along the steep, nonurbanised slopes to the very low connectivity of the highly-urbanised plateaus. If the frogs manage to leave the stream, they become hampered or blocked by traffic and buildings. This complex topology constrains connectivity amongst the invaded locations and the few other water bodies within the study area. According to our model, topography and urban areas are therefore sufficient to explain the non-colonisation of the three nearest streams -Sassoeiros to the west, Jamor to the east and Porto Salvo in between the two colonised rivers (Fig. 1C).

Past dispersal routes
Due to road and building constructions after the year 2000, the maps that we used for this work may not depict all the dispersal corridors that were available in the 1980s and 1990s. However, our models show that Barcarena was very probably invaded from Laje by frogs that dispersed overland and used the golf-course ponds as stepping stones. According to the model, dispersing frogs may have used two small tributaries to reach the golf-course ponds. The northernmost tributary is the strongest candidate as a past dispersal route, given the large population that was found there. The golf ponds were dug and filled in 2002 and are located exactly in the single suitable corridor identified. As noted by other authors on less hilly landscapes (cf. Faraone et al. 2008), small water bodies can be used as stepping stones during the rainy season to reach suitable habitat. However, we cannot fully exclude other non-accounted factors, such as the intentional release by amphibian keepers.

Areas of high risk
All factors, low habitat quality, restricted availability of water bodies and hampered dispersal ability, probably explain the comparatively slow invasion of X. laevis in Portugal. Still, the species has managed to colonise two rivers and this work suggests that it used artificial water bodies as stepping stones on a terrestrial pathway in a densely-urbanised area, highlighting the risks of further invasions.
Possible stepping stones for dispersal into other streams should be identified, monitored and, if possible, altered (e.g. by encircling them with walls, since there are no natural ponds in the area) to hamper further overland spread. The Jamor River basin and Carregueira Mountain are two semi-natural regions not yet invaded, located east and northeast of the currently-invaded area ( Figure 1A). Jamor River is apparently protected by a high landscape resistance (Fig. 1C) and is, therefore, not naturally reachable by terrestrially dispersing X. laevis. However, landscape resistance is low at Carregueira, a small forested mountain range with farms and golf-courses (including lakes), which also contains the headwaters of Jamor River. Currently, this area is not in the range of the connectivity kernel, but a X. laevis invasion is possible starting at its western edge. In the northern limit of the Barcarena population, there are also several small water bodies that could be used as possible stepping stones into novel areas.
The current eradication plan for this species in Portugal (Sousa et al. 2018) can be informed by this study. Regular monitoring of the water bodies that are within reach of the species, according to our model, is strongly advised. If possible, toad barriers could be built on the main pathways for overland dispersal, particularly at the edges of Carregueira and around the northern isolated population, which should effectively block further spread. As the number of colonised sites is reduced by the eradication programme, we recommend to verify (and if possible, to block) all paths within a radius of 5 km around the colonised water bodies (double the species currently known maximum terrestrial dispersal ability during a dispersal event) to minimise expansion risk. Our results are also relevant for other countries where X. laevis occurs, highlighting the importance of blocking strategic overland routes of dispersal, either by using toad fences or by draining ponds that may be used as stepping stones.

Advantages and limitations
The fine-scale remote sensing derived resistance surfaces, based on NDVI and NDWI, in combination with elevation layers, allowed us to reconstruct potential past dispersal routes between the two invaded rivers and highlighted areas at high risk of invasion. This provides a detailed map highlighting areas which are threatened by invasion and knowledge of potential corridors for the invasive species. However, the computational power and time needed for this method increases with the number of starting points and with the resolution of raster layers. Furthermore, this approach is based on speciesspecific knowledge about biology and physiology and model accuracy strongly depends on evaluation by experts. Some very fine scale dispersal barriers may remain undetected by remote sensing, such as waterfalls with seasonally varying intensities or smooth walls. These landscape features may further restrict the dispersal potential on a local scale.