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 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².

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).

We obtained high resolution multispectral satellite imagery containing the invaded Portuguese distribution range (625 km² × 4 title IDs = 2500 km²) 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.

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.

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.

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 (De Villiers and Measey 2017). 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.

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, non-urbanised 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).

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.

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.

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 species-specific 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.