Combining geostatistical and biotic interaction model to predict amphibian refuges under crayfish invasion across dendritic stream networks

Biological invasions are pervasive in freshwater ecosystems, often causing native species to contract into areas that remain largely free from invasive species impacts. Predicting the location of such ecological refuges is challenging, because they are shaped by the habitat requirements of native and invasive species, their biotic interactions, and the spatial and temporal invasion patterns. Here, we investigated the spatial distribution and environmental drivers of refuges from invasion in river systems, by considering biotic interactions in geostatistical models accounting for stream network topology. We focused on Mediterranean amphibians negatively impacted by the invasive crayfishes Procambarus clarkii and Pacifastacus leniusculus.


| INTRODUC TI ON
Biological invasions are pervasive in freshwater ecosystems, where they are major drivers of native species declines (Strayer, 2010;Walsh, Carpenter, & Vander Zanden, 2016). Addressing this threat is challenging, because once fully established the control of invasive species is often nearly impossible, which limits the management options to protect native species. In some circumstances, the impacts of biological invasions may be partly offset by the presence of ecological refuges, which are habitats where a species can retreat, persist in for up to a few decades, and eventually expand from under changing environmental conditions (Davis, Pavlova, Thompson, & Sunnucks, 2013).
Such refuges correspond to freshwater habitats unsuitable for invasive species, or areas where their spread is prevented by physical barriers such as waterfalls or culverts (Kerby, Riley, Kats, & Wilson, 2005;Rahel, 2013). Refuges may thus allow the persistence of at least some remnant populations of native species (e.g. Chapman et al., 1996;Grabowski, Bacela, Konopacka, & Jazdzewski, 2009;Habit et al., 2010;Radinger, Alcaraz-Hernández, & García-Berthou, 2019), making it a priority to understand where, why and how refuges can contribute to species conservation under biological invasion.
Species distribution models (SDM) incorporating biotic interactions provide a simple framework to quantify how one or more species influence the distribution of others (e.g. Wisz et al., 2013), making them useful to predict the location and drivers of refuges from invasions. A straightforward approach is to take the distribution patterns of invasive species together with abiotic variables to model the occurrence of native species, and then use the ensuing models to predict the distribution of refuges under current and future invasion scenarios (Araújo & Luoto, 2007;Wisz et al., 2013).
One problem is that SDMs for geographical range prediction assume equilibrium between species distribution and the environment, which is unwarranted when modelling range contractions by native species in face of biological invasion (Elith, Kearney, & Phillips, 2010;De Marco, Diniz, & Bini, 2008;Václavík et al., 2012). A native species may occur in areas that will latter become unsuitable due to the expansion of an invasive species, or it may eventually be able to coexist only temporarily with the invasive species due to time lags in negative impacts (Crooks, 2005). In either case, SDMs built using snapshots of species distributions may overestimate the extent of refuges, eventually misdirecting conservation efforts towards areas where native species persistence is unlikely.
Incorporating predictors describing spatial autocorrelations in invasive species occurrences to account for unmeasured dispersal and colonization processes may help mitigating, albeit not solving, problems associated with non-equilibrium conditions in SDMs (De Marco et al., 2008;Václavík et al., 2012;Filipe, Quaglietta, Ferreira, Magalhães, & Beja, 2017). For alien species invading river ecosystems, SDMs can be improved using geostatistical models accounting for spatial dependencies in physical and ecological processes across stream networks (Filipe et al., 2017;Lois et al., 2015;Lois & Cowley 2017). These models are similar to conventional mixed models, with species occurrence modelled in relation to environmental variables using a logistic function, and spatial autocorrelation considered in the random errors Peterson et al., 2013;Ver Hoef, Peterson, & Theobald, 2006).  . This approach can easily incorporate biotic interactions by including predictors describing the occurrence or abundance of potentially interacting species in the fixed component (Lois et al., 2015, Lois & Cowley 2017. Another possibility is to develop a geostatistical model for the invasive species itself, and then use the fitted response (i.e. the probability of occurrence) in the native species model. This should be useful for predicting the location of refuges, as it would consider not only the current distribution of the invasive species, but also suitable areas that will eventually be colonized during the expansion process.
This study investigates the location and environmental drivers of refuges in dendritic stream networks, combining biotic interactions and geostatistical modelling to predict their spatial distribution under current and future scenarios of invasive species expansion. We focused on interactions between amphibians and the exotic crayfish Procambarus clarkii and Pacifastacus leniusculus in the Iberian Peninsula, where there are no native crayfish (Clavero, Nores, Kubersky-Piredda, & Centeno-Cuadros, 2016). These crayfishes are among the most widespread and damaging aquatic invaders (Lodge et al., 2012;Twardochleb, Olden, & Larson, 2013), which have expanded widely in Iberia since the 1970s due to multiple introductions for commercial purposes and subsequent natural dispersal (Bernardo, Costa, Bruxelas, & Teixeira, 2011;Clavero, 2016;Gutiérrez-Yurrita et al., 1999). Invasive crayfish predate on amphibian eggs and larvae (Axelsson, Brönmark, Sidenmark, & Nyström, 1997;Cruz & Rebelo, 2005;Cruz, Rebelo, & Crespo, 2006b;Gamradt & Kats, 1996), and seem to have strong negative impacts on native amphibian populations in Iberian waters (Cruz et al., 2006a;Cruz, Rebelo, et al., 2006b;Cruz, Segurado, Sousa, & Rebelo, 2008) and elsewhere (Ficetola et al., 2011). The main amphibian refuges are probably temporary ponds far from permanent waters (Beja & Alcazar, 2003;Cruz, Rebelo, et al., 2006b;Ferreira & Beja, 2013), Alien invasive species, biological invasions, ecological refuges, frog, geostatistics, newt, species distribution models, stream ecology, toad where crayfish cannot persist (Cruz & Rebelo, 2007). Refuges may also exist in small and intermittent Mediterranean streams, which often hold rich amphibian communities (de Vries & Marco, 2017) and where crayfishes are usually absent (e.g. Cruz & Rebelo, 2007;Filipe et al., 2017;Gil-Sánchez & Alba-Tercedor, 2002). To investigate amphibian refuges from crayfish invasion, we (a) made a detailed survey of amphibian occurrence in a Mediterranean watershed; (b) developed geostatistical models relating the occurrence of each amphibian species to environmental variables, the probabilities of crayfish occurrence (Filipe et al., 2017) and spatial dependencies; and predicted the spatial distribution of refuges under (c) current and (d) future scenarios of crayfish expansion. Our study can help improve conservation strategies for amphibians negatively affected by crayfish invasions, and more generally, it provides a framework for identifying refuges from biological invasions in dendritic stream networks.

| Study area
The study was conducted in the river Sabor watershed (NE Portugal; N41º090-42º000, W7º150-6º150; Figure S1). Human population is low (8.5-28.7 inhabitants/km 2 ; https://www.porda ta.pt/Munic ipios) following a process of land abandonment since the 1970s (Azevedo, Moreira, Castro, & Loureiro, 2011). Land cover is dominated by extensive agriculture and pastureland, forest plantations, and natural vegetation (Caetano, Marcelino, & C. Igreja e I. Girão, 2018;Hoelzer, 2003). A large proportion of the watershed is included in the Natura 2000 network (Costa, Monteiro-Henriques, Neto, Arsénio, & Aguiar, 2007). The watershed covers a wide range of elevations (100-1,500 m above sea level), total annual precipitation (443-1,163 mm) and mean annual temperature (6.9-15.6°C). The climate is Mediterranean, with precipitation concentrated in October-March and virtually none in June-August. Most small streams dry out or become reduced to a series of disconnected pools during the dry months, though the main watercourse and the largest tributaries are permanent (Ferreira, Filipe, Bardos, Magalhães, & Beja, 2016). Two large hydroelectric reservoirs were built and filled just before this study, but otherwise the river is largely free-flowing. Crayfish were first reported in the Sabor watershed in the 1990s, with P. clarkii probably introduced by local people, while P. leniusculus was introduced in 1994 by Spanish authorities (Bernardo et al., 2011). P. clarkii is far more widespread than P. leniusculus, but both species seem to still be spreading, possibly through natural dispersal along the stream network (Anastácio et al., 2015;Bernardo et al., 2011).

| Study design
The study was designed to obtain a comprehensive snapshot of stream-dwelling amphibian distributions, considering both species such as Iberian frog (Rana iberica) and Iberian green frog (Pelophylax perezi) that occupy streams during their entire life cycle, and species such as fire salamander (Salamandra salamandra) and midwife toads (Alytes obstetricans and A. cisternasii) that have both terrestrial and aquatic phases, occupying streams mostly during the breeding and larval development periods. To adequately cover all species, surveys encompassed the main environmental gradients represented in the watershed, and they were carried out monthly during one year to account for differences in activity peaks and breeding phenology across species (e.g. Díaz-Paniagua, 1992;Ferreira & Beja, 2013). Yet, because it was logistically unfeasible to sample a sufficiently large number of sites each month, the sampling effort was distributed over the year, with different sites sampled in different months. Considering these constraints, we initially selected >200 potential sampling sites, based on previous studies (Ferreira et al., 2016;Filipe et al., 2017;Quaglietta, Paupério, Martins, Alves, & Beja, 2018) and new field surveys. Sites were constrained to cover all Strahler stream orders and to be at >1 km from each other. All sites were in wadable stream reaches (water depth < 1.20) to facilitate amphibian surveys. From the overall set of potential sampling sites, we selected each month a subset of 30 sites, following a stratified random procedure to guarantee that comparable environmental conditions were covered each month, and thus avoiding space × time interactions. To do this, we divided the Sabor watershed into three sub-basins ( Figure S1), and randomly selected each month two sites of each Strahler stream order represented in each sub-basin. A total of 168 sites were surveyed ( Figure S1), with each site visited at most twice, except sites in higher orders that were visited more often because they were relatively scarce in the watershed. The river network, sub-basins and stream orders were obtained from CCM 2.1 (Catchment Characterization and Modelling database; Vogt et al., 2007). A detailed workflow of the procedures used to analyse the data is provided in the Supplementary Methods in Supporting Information.

| Field sampling
Sampling was carried out monthly in 2015, except in May due to logistic constraints. At each sampling site and date, a 200-m stream reach was thoroughly surveyed for amphibians, including both adults and larvae. The survey was conducted by two observers walking slowly along the banks or wading in shallow water along the stream. Observers used dip nets to collect aquatic larvae and adults, and they systematically searched the stream banks for terrestrial adults, using torches where necessary to survey cavities and other shaded areas. All amphibians found were identified to species level in situ and released thereafter. In a few cases, small unidentified larvae were preserved and identified in the laboratory. During surveys, crayfish were also recorded and identified to species, and all individuals collected were eliminated following guidelines established by the Portuguese biodiversity conservation agency.

| Environmental and spatial variables
Sampling sites were characterized using variables potentially affecting the distribution of stream-dwelling amphibians that could be extracted from digital maps (e.g. Cruz & Rebelo, 2007;de Vries & Marco, 2017), making it possible to predict species distributions across the entire watershed in relation to potential expansions of crayfish ranges. Each site was characterized using four environmental variables (elevation [Alt], total annual precipitation [Prec], Strahler's stream order [SO] and the probability of water presence during the dry season [Water] and two variables describing the potential for biotic interactions between amphibians and either P. clarkii [Pclar] or P. leniusculus [Plen] (Cruz, Rebelo, et al., 2006b;de Vries & Marco, 2017). We also included a multiplicative interaction between elevation and Strahler's stream order (SOxAlt), which was used to distinguish between small streams in the lowlands and small streams in mountain areas. Initially, we also considered other climatic variables, but they were discarded because of strong correlations with precipitation and/or elevation. Although features of the surrounding landscape are known to affect stream-dwelling amphibians (Ficetola et al., 2011;Riley et al., 2005), these were not considered because preliminary analysis showed very minor effects of land cover variables in our study area, possibly due to the dominance of natural vegetation and low-intensity land uses that are generally suitable for amphibians. Furthermore, adding land cover variables often caused model instability and convergence problems, possibly due to redundancies with other environmental variables already included in the models.
All variables were computed in a Geographic Information System (GIS) using ArcGis (ESRI, 2016). Elevation was taken from a DEM built from 1:25,000 topographic maps as in Ferreira et al. (2016). Total annual precipitation was extracted from WordClim 2 with a 30' (≈1 km) resolution (Fick & Hijmans, 2017). Strahler's stream order was used as a proxy for habitat size and heterogeneity (Ferreira et al., 2016;Hughes, Kaufmann, & Weber, 2011), and it was extracted from CCM 2.1, which is based on a 100-m resolution digital elevation model (DEM; Vogt et al., 2007). The probability of water presence in the dry period was used because many Mediterranean amphibian species are associated with temporary water bodies (Beja & Alcazar, 2003;Ferreira & Beja, 2013;de Vries & Marco, 2017), and it was extracted from a model developed in a previous study (Ferreira et al., 2016). Variables describing biotic interactions were specified considering the probability of occurrence of either P. clarkii or P. leniusculus, extracted from a previous geostatistical modelling of crayfish distribution in the Sabor watershed (Filipe et al., 2017). These models were built using electrofishing data collected on 167 sites in summer 2012, and they showed that crayfish distributions were mainly associated with stream order, elevation and spatial dependencies across the stream network (Filipe et al., 2017). The models had reasonable predictive accuracy, for both P. clarkii (AUC = 0.963) and P. leniusculus (AUC = 0.823). Probabilities derived from distribution models were used instead of the actual crayfish presences/absences recorded at sampling sites, in order to project the distribution of each species across the entire watershed, as well as to build scenarios of future crayfish expansion. All variables were standardized to have a mean of zero and a standard deviation of one to improve the interpretability of model coefficients (Schielzeth, 2010), and we screened for outliers and influential points that might bias coefficient estimates.
Spatial data necessary to account for spatial autocorrelation (see below) were obtained in a GIS using the Sabor watershed stream network extracted from CCM2.1 (Vogt et al., 2007), and the layer of survey sites. Estimates included the Euclidean and hydrologic distances (total and downstream hydrologic distances) between every pair of sampling sites . To deal with confluences in tail-up models, we also estimated watershed areas to weight the relative influence of the branching upstream segments

| Distribution modelling
Distribution models (SDM) were developed considering for each species only the occurrence data from the sampling months encompassing its aquatic phase, and thus the periods when the species is detectable during stream surveys (Table S1). For instance, we considered data from all sampling months for R. iberica and P. perezi, because they are largely aquatic and strongly attached to stream habitats all year round, while we only considered data from April-November and October-April for A. obstetricans and S. Salamandra, respectively, corresponding to their aquatic phases. This approach aimed at avoiding false negatives caused by species not using potentially suitable stream habitats during the terrestrial phase.
For each species, we developed three logistic models relating its presence/absence to either (a) only environmental predictors, (b) only biotic predictors or (c) environmental + biotic predictors, and a (d) geostatistical distribution model (Peterson et al., 2013). The logistic models were used to evaluate how considering biotic interactions affected species perceived responses to environmental variables, and as preliminary steps in geostatistical model building. The geostatistical model for each species included a fixed component corresponding to a logistic function linking its probability of occurrence to environmental and biotic predictors, and random components accounting for spatial dependencies in the stream network Peterson et al., 2013;Ver Hoef et al., 2006). As random components, we considered the Euclidean model, assuming that spatial dependencies among sites can occur overland, and the tail-up and tail-down covariance models, assuming that spatial dependencies can also occur along the hydrological network independently of flow and/or only between flow-connected sites, respectively .
To build the logistic models for each set of predictors and species, we considered all combinations of predictors of each set and retained for inference the best subset model minimizing AIC (Murtaugh, 2009). Autocorrelation in model residuals was visualized using Torgegrams, depicting how semivariance in the residuals of the best logistic models between pairs of sampling sites changed in relation to their hydrologic distances (Peterson et al., 2013). In Torgegrams, increasing semivariance reflects declining spatial dependency between points. The fixed component of the geostatistical model for each species was then built considering the predictors included in the best environmental + biotic logistic model. The structure of the random components was assessed using the residuals of the environmental + biotic model following Quaglietta et al. (2018), by testing all combinations of alternative functions available in the R package "SSN" (Ver Hoef, Peterson, Clifford, & Shah, 2014) for the Euclidean and the hydrologic autocovariance functions, and retaining the model structure minimizing AIC. We then combined the variables selected for the fixed component with the best spatial structure to build the final model for each species. In all logistic models and in the fixed component of the geostatistical model, we considered significance level for individual predictors at p < .10, to reduce the likelihood of type II errors and thus the probability of missing true-negative effects of crayfish invasion.
The crayfish models previously developed by Filipe et al. (2017) were validated with the presence/absence data from the 2015 survey, using the area under the receiver operating curve (AUC) (Allouche, Tsoar, & Kadmon, 2006). Discrimination ability of amphibian models was estimated using predictions obtained by the "leaveone-out" cross-validation method, considering the overall prediction success, the AUC, Cohen's kappa and the true skill statistics (TSS) (Allouche et al., 2006;Václavíık, Kupfer & Meentemeyer, 2012).

| Distribution mapping under current and future scenarios
To map the predicted distribution of each amphibian species under current conditions, we projected the distribution models on the stream network of the entire Sabor watershed. First, we divided the stream network into segments of a maximum length of 1,000 m using ArcGIS desktop (ESRI, 2016), and we extracted the value of environmental variables from the centroid of each segment. We then predicted the probability of each species occurring in the segment using universal kriging within the "SSN" package . We considered segments occupied if the predicted occupancy probability was above the species prevalence threshold (Liu et al., 2005).
To simulate how crayfish expansion might affect amphibians, we changed the value of variables describing biotic interactions assuming a twofold, threefold and fivefold increase in the relative risk (Rr) of each crayfish species occurring at each site, using as baseline the predictions from the geostatistical models of Filipe et al. (2017). Future distributions of each amphibian species were predicted using either the non-spatial environmental + biotic logistic model or the spatial geostatistical model, which reflect different assumptions on range change processes (Record, Fitzpatrick, Finley, Veloz, & Ellison, 2013). The non-spatial model assumes that amphibian distributions will change along with changes in crayfish occurrence, irrespective of amphibian current distributions. However, spatial structure is still implicit in predictions, because probabilities of crayfish occurrence across the dendritic stream network were themselves predicted using geostatistical models (Filipe et al., 2017).
In the geostatistical model, spatial random effects act to draw the projected distributions back towards the observed distribution used to calibrate the model (Record et al., 2013). Therefore, the current and future distributions will be similar, unless there are strong negative effects of biotic interactions. Prediction of future distributions was only made for amphibian species showing significant negative effects of crayfish occurrence. We did not consider climate change effects due to uncertainties regarding how climate will change in our relatively small area and how crayfish and amphibians will respond to such changes, though this should be the subject of further research due to potential interactions between climate change and biological invasion (Hulme, 2017).

| RE SULTS
We detected a total of 10 amphibian species, the most widespread (>20% of sites) of which were P. perezi (69%), S. salamandra (28%) and R. iberica (26%) ( Table S1). The frogs Discoglossus galganoi and Hyla     The distributions of these species correspond to streams with low probability of occurrence of both invasive crayfish, but that were predicted to be progressively colonized under the invasion scenarios ( Figures S4 and S5). Expansion of crayfish through the stream network was predicted to decrease the length of potential stream habitat for amphibians, with each species becoming progressively more confined to first-and second-order streams ( Figure 2, Table   S2). Predictions using spatial models suggested that the length of habitat of A. cisternasii will decline up to about 30% due to the expansion of P. leniusculus, while a reduction of about 7% was predicted for A. obstetricans due to the expansion of P. clarkii. The potential habitat for R. iberica and S. salamandra was expected to decline by about 20%, due to the joint expansion of P. leniusculus and P. clarkii. Non-spatial models predicted even larger declines for A. cisternasii (up to 69.9%), A. obstetricans (51.3%) and R. iberica (53.9%), but not as much for S. salamandra (24.1%) (Figure 2, Table   S2). For the latter species, there were upstream areas in the far north that were predicted to be occupied by the non-spatial model but not by the spatial model (Figure 2), suggesting that environmentally suitable habitats may remain unoccupied due to spatial processes.

| D ISCUSS I ON
Our study showed the value of combining geostatistical and biotic interaction modelling to quantify the spatial consequences of biological invasions on native species in dendritic stream networks, and to predict the spatial distribution of ecological refuges under current and future invasion scenarios. Using this approach, we confirmed the strongly negative interactions between invasive crayfish and amphibians (Cruz, Pascoal, et al., 2006a;Cruz, Rebelo, et al., 2006b;Ficetola et al., 2011;Riley et al., 2005)   represent key refuges from crayfish invasion for many amphibian species, as these streams dry out for more or less extended periods during the dry season (Ferreira et al., 2016) and are thus expected to remain largely free from crayfish impacts (Cruz & Rebelo, 2007;Filipe et al., 2017). Overall, our study reinforces the conservation importance of stream headwaters in the Mediterranean region, which are increasingly perceived to play key roles as refuges from biological invasions and other human-mediated disturbances, both for amphibians (de Vries & Marco, 2017) and other vulnerable species (Quaglietta et al., 2018;Sousa et al., 2019). Our approach may be applied to other aquatic species, with major implications for conservation and management by permitting a better identification of areas acting as ecological refuges under biological invasion.

| Study limitations
Although our study had some limitations, it is unlikely that they affected our main conclusions in any significant way. One potential problem is that we sampled streams from first to sixth orders, and it may be argued that the prevalence of amphibian species may be underestimated in larger streams due to lower detectability. Although this can be addressed by modelling occupancy while controlling for detectability (e.g. MacKenzie et al., 2017), this was not possible in our case because occupancy-detection models accounting for hydrological tail-up and tail-down spatial dependencies have yet to be developed. To deal with this problem, we have surveyed only wadable streams and increased the sampling effort in larger streams, which should have contributed to achieve comparable detectability across stream orders. This is supported by the higher prevalence of P. perezi and the crayfish P. clarkii in higher than lower orders, which suggest that we did not miss species known to occur in larger streams (Cruz, Rebelo, et al., 2006b;Filipe et al., 2017). Moreover, species distribution patterns observed in our study were consistent with those of others focusing on stream-dwelling amphibians in Iberia (Cruz, Rebelo, et al., 2006b;de Vries & Marco, 2017), thereby suggesting that they were not artefacts shaped by sampling biases.
Another potential problem is that amphibian models were based on large-scale variables, while ignoring local drivers such as the structure and composition of riparian vegetation, the land uses surrounding streams and water quality (Crawford & Semlitsch, 2007;Guzy, Halloran, Homyack, & Willson, 2019;Riley et al., 2005). Missing these variables might have reduced the predictive ability of our models, but we believe that the variables considered are relevant to investigate large-scale distribution patterns, as observed for instance F I G U R E 2 Maps of potential distributions of four amphibian species under the worst case scenario of future crayfish expansion in the river Sabor watershed (NE Portugal). The scenario was built considering a fivefold increase in the relative risk of crayfish occurrence (both P. clarkii and P. leniusculus) at each stream segment in relation to the baseline scenario corresponding to the predicted distribution of each species in 2012 (Filipe et al., 2017). In each map, we indicate the waterlines where potential habitat will remain available (suitable habitat), and those where potential habitat will be lost (lost habitat) in relation to the baseline scenario. Maps were produced using predictions from either spatial or nonspatial (i.e. including only environmental effects and biotic interactions) models [Correction statement added on 05 May 2020 after first online publication: Figure  2 was previously merged and has been divided into two parts in this version.]

Alytes cisternasii
in other species modelled in the Sabor watershed (Ferreira et al., 2016;Filipe et al., 2017;Quaglietta et al., 2018). Nevertheless, while our models should be informative to understand broad changes in amphibian distributions in relation to crayfish invasion, they may be less useful to predict whether a given species will be present at any given site, which may be strongly affected by more local environmental conditions.

| Effects of biotic interaction and environmental effects on amphibian distributions
The negative responses of amphibians to invasive crayfish observed in the Sabor watershed were comparable to those reported elsewhere (Cruz, Pascoal, et al., 2006a;Gamradt & Kats, 1996;Gamradt, Kats, & Anzalone, 1997;Girdner et al., 2018;Nyström, Birkedal, Dahlberg, & Brönmark, 2002;Nyström, Svensson, Lardner, Brönmark, & Granéli, 2001). However, impacts of P. leniusculus are reported here for the first time in Iberia, though this crayfish was already associated with local amphibian declines in Sweden (Nyström et al., 2002(Nyström et al., , 2001. We also confirmed that crayfish impacts seem to be particularly strong on urodela (salamander and newts), possibly because their eggs and larvae are highly vulnerable to introduced predators (Cruz, Rebelo, et al., 2006b;Gamradt & Kats, 1996;Gamradt e al. 1997;Girdner et al., 2018). Mediterranean amphibians may be especially vulnerable to crayfish because many species are adapted to live in water bodies that dry out in summer and are naturally free from fish and other large predators (e.g. Beja & Alcazar, 2003;Ferreira & Beja, 2013), and thus may be less adapted to cope with crayfish  predation than species living in permanent waters (Cruz, Rebelo, et al., 2006b;Nunes et al., 2011). This is supported by the lack of negative effects on P. perezi, which is known to thrive in permanent waters where predators are abundant (Beja & Alcazar, 2003;Cruz, Rebelo, et al., 2006b;Ferreira & Beja, 2013). We also found no negative effects on B. spinosus, which is widespread in permanent waters and seems to be less vulnerable to predators due to its toxic eggs and larvae (Cruz & Rebelo, 2005), though a previous study reported negative impacts of P. clarkii (Cruz, Rebelo, et al., 2006b). In contrast, we found negative impacts on A. cisternasii, which is associated with more permanent water bodies and was found previously to be unaffected by crayfish (Cruz, Rebelo, et al., 2006b). However, that study was carried out in an area where only P. clarkii occurred, while in our study, we only found significant negative effects for P. leniusculus, suggesting that impacts may differ among crayfish species. P. leniusculus also showed strong negative effects on R. iberica, possibly because this species occurs primarily in mountainous streams largely unsuitable for P. clarkii (Filipe et al., 2017). These results suggest that invasion by multiple crayfish species may be more serious than invasion by a single species, by increasing the types of habitats invaded (Filipe et al., 2017) and the number of species vulnerable to predation.
Although our study showed that many amphibian species were associated with stream headwaters, different species occurred in different areas, probably due to differences in ecological requirements. For instance, environmental models suggested that while both A. cisternasii and A. obstetricans were mainly found in lower and middle stream orders, the former favoured areas with lower precipitation at low elevation, while the reverse was found for the latter. This probably explains their largely parapatric distributions in the study area, as observed at broader spatial scales (Reino et al., 2017). Models for R. iberica also confirmed their preference for small permanently flowing streams in mountainous areas (Bosch, Rincón, Boyero, & Martínez-Solano, 2006;Rodríguez-Prieto & Fernández-Juricic, 2005). Many of these environmental effects, however, were lost from the best models or became non-significant, once biotic interactions were included. In some cases, there were also changes in the significant environmental effects, with for instance the selection of permanently flowing waters by A. obstetricans only becoming apparent after controlling for the effects of P. clarkii. These results support the idea that the two exotic crayfish are key drivers constraining amphibian distributions in our study area, limiting the range of environmental conditions where they can be found.

| The role of spatial dependencies across the stream network
Incorporating spatial covariance structure greatly enhanced the distribution models, with both Euclidean and hydrologic distances often included in the best models. The Euclidean component was important for A. cisternasii, P. perezi, S. salamandra and T. marmoratus, suggesting that adjacent streams have more similar occupancy status than streams farther apart, which may be due to the dispersal of individuals overland (e.g. Semlitsch, 2008)  Crawford & Semlitsch, 2007;Guzy et al., 2019;Riley et al., 2005).
Future studies combining drivers operating at landscape and local scales should thus be developed, which would likely improve predictions on potentially favourable areas across the watershed.

| Predicting amphibian refuges under crayfish invasion
Despite the limited explanatory power of the spatial distribution models, the mapping of predicted distributions clearly showed that further crayfish expansions will likely result in amphibian range contractions, with populations becoming progressively more encroached in lower order streams. This could be inferred quantitatively for four species with geostatistical models with sufficient predictive ability (A. cisternasii, A. obstetricans, R. perezi and S. salamandra), but will probably occur also for the other two species showing negative associations with crayfish occurrence (L. boscai and T. marmoratus) in ours and other studies (Cruz, Rebelo, et al., 2006b). The non-spatial models (i.e. including only environment + biotic interactions) predicted even stronger declines in the availability of potential habitats, particularly for A. cisternasii, A. obstetricans and R. iberica. This is because spatial models anchor model predictions to the current distribution of each species, and thus may be regarded as conservative because they add inertia against abrupt changes in distribution driven by environmental factors (Record et al., 2013). In contrast, the non-spatial models are only driven by changes in crayfish occurrence irrespective of the spatial structure in current amphibian distribution, thereby disregarding possible spatially structured population processes shaping amphibian species distributions (Record et al., 2013). Therefore, we expect that the extent of range contractions will be somewhere in-between the predictions of spatial and non-spatial models, which may thus be substantial for A. cisternasii (up to 30.6%-69.6%), A. obstetricans (6.8%-51.3%) and R. iberica (20.6%-53.9%), though only moderate to S. salamandra (22.3%-24.1%). The consequences of such changes for population persistence should be evaluated in future studies, as it is likely that the risk of local extinctions will be high for small and isolated populations confined to headwater streams.

| Conservation and management implications
Invasion by alien crayfish is a major cause of concern for amphibian conservation (e.g. Cruz, Pascoal, et al., 2006a;Gamradt & Kats, 1996;Gamradt et al., 1997;Girdner et al., 2018;Nyström et al., 2001Nyström et al., , 2002. Addressing this problem is challenging, because once established invasive crayfish populations are virtually impossible to eradicate, and thus, remediation of ecosystems invaded by crayfish has met very limited success (Gherardi, Aquiloni, Diéguez-Uribeondo, & Tricarico, 2011;Stebbing, Longshaw, & Scott, 2014). Furthermore, many invasive crayfish species are still expanding within and across watersheds (e.g. Bernardo et al., 2011;Kouba, Petrusek, & Kozák, 2014), and so the problem is likely to get worse in the future. In this context, our study suggests that stream headwaters may be critical for the persistence of many stream-dwelling amphibian species, at least in the Mediterranean region, as they often hold diverse amphibian communities (de Vries & Marco, 2017) and are likely to provide refuges with minimal or no crayfish impacts (Filipe et al., 2017). These headwaters correspond not only to small order streams in mountainous areas, as those inhabited for instance by A. obstetricans and R. perezi, but also to small temporary streams at lower elevation, which seem to be preferred by species such as A. cisternasii. Overall, therefore, headwater streams should be regarded as priority targets for conservation, requiring the preservation of habitat conditions compatible with amphibian persistence. Although the ecological requirements of stream-dwelling amphibians are poorly known in the Iberian Peninsula (de Vries & Marco, 2017), it is likely that conservation efforts should target preserving water quality, natural flow regimes, well-developed riparian vegetation and suitable terrestrial habitats (e.g. Crawford & Semlitsch, 2007;Guzy et al., 2019;Riley et al., 2005). Furthermore, efforts should be made to avoid the colonization of headwater refuges by invasive crayfish such as P. leniusculus, which is fast expanding into new areas (Anastácio et al., 2019;Bernardo et al., 2011), and may be able to colonize mountainous headwater streams inhabited by endemic amphibians such as R. iberica (Filipe et al., 2017). This would require monitoring crayfish populations in key amphibian refuges, which should be used to trigger careful management programmes if the risk of negative impacts become unacceptably high, involving for instance the implementation of control or eradication programs, and/or the introduction of physical barriers to crayfish dispersal (Gherardi et al., 2011;Sousa et al., 2019;Stebbing et al., 2014). Although such conservation measures may require considerable efforts and may only be applicable in some areas, maintaining stream headwaters free of invasive crayfish should have major conservation benefits for a range of endangered species (Quaglietta et al., 2018;Sousa et al., 2019;de Vries & Marco, 2017) and for aquatic biodiversity in general (Finn et al., 2011;Meyer et al., 2007).

| CON CLUS IONS
Under biological invasion, many native species are becoming confined to refuges where invasive species are still absent or scarce, and thus may hold remnant populations of high conservation value (e.g. Chapman et al., 1996;Grabowski et al., 2009;Habit et al., 2010;Radinger et al., 2019). This study provides a framework to predict the location and environmental drivers of such refuges, using geostatistical tools to model native species responses to exotic species while controlling for environmental effects and spatial dependencies across dendritic stream networks (Filipe et al., 2017;Peterson et al., 2013).
This approach is relatively simple and can be used where only snapshot surveys on the occurrence patterns of native and invasive species are available, though it can be easily extended to deal with data on distributional dynamics (Quaglietta et al., 2018) and additional complexities such as climate change (Peterson et al., 2013). This framework may be generally useful to understand the distributional consequences of interactions between native and invasive species, providing information on the location of potential refuges where conservation efforts should concentrate, and on management actions required to enhance the persistence of remnant populations within refuges.

ACK N OWLED G EM ENTS
This work was carried out in the scope of the Baixo Sabor Long Term

DATA AVA I L A B I L I T Y S TAT E M E N T
The SSN files including information on species presence/absence and on the spatial topology of the stream network and sampling sites are provided in Dryad (https://doi.org/10.5061/dryad.mw6m9 05tb), together with the R scripts and additional files needed to repeat all analysis using the SSN package. The raw data on amphibian occurrences are available from the authors upon request.