Network-scale effects of invasive species on spatially-structured amphibian populations

year. Hence, SSP-level analyses can identify effects that would remain unnoticed when focussing on single patches. Invasive species can affect population dynamics even in not invaded patches, through the degradation of subpopulation networks. Patch-scale assessments of the impact of invasive species can thus be insufficient: predicting the long-term interplay between invasive and native populations requires landscape-level approaches accounting for the complexity of spatial interactions.


Introduction
The idea that population processes heavily depend on spatial connections has revolutionized our understanding of population dynamics. Levins (1969) first posited that the colonisation -extinction dynamics within networks of subpopulations are affected by their migration pattern (metapopulations dynamics). In turn, migration is influenced by the spatial configuration of patches, thus the metapopulation concept has boosted research on the effects of fragmentation on the demography of populations, leading to models that help to predict the distribution patterns and dynamics of species in real landscapes (Hanski 1998, De Roissart et al. 2015. In most ecosystems, the variation of patch-scale characteristics has significant effects on local demographic rates (Hanski 1998), still the dynamics of subpopulations living in discrete habitats are also determined by processes acting at the landscape level, such as connectivity and dispersal (Friedenberg 2003, Semlitsch 2008, De Roissart et al. 2015. This can lead to source-sink dynamics, where the viability of a subpopulation often depends by the immigration/emigration from nearby populations, and some patches cannot sustain local populations without immigration from connected patches (Hale et al. 2015, Seward et al. 2019. For instance, sites close to populations with high net reproductive rate (which act as sources) can have more positive population dynamics than expected solely on the basis of their in-site recruitment (Hanski 1998, Figuerola and Green 2002, Hanski and Gaggiotti 2004. These ideas have prompted research on the effects of spatial relationships within spatially-structured populations (SSP) (Fortin andAgrawal 2005, Revilla andWiegand 2008). Spatially-structured populations are groups of local populations connected by the exchange of dispersing individuals and comprise a continuum range of categories classically considered within the metapopulation framework, such as patchy populations, metapopulations and sourcesink networks (Revilla andWiegand 2008, Cayuela et al. 2019). Immigration and connectivity are key determinants of population processes at multiple scales. Therefore, the analysis of SSP and the importance of their spatial connections have been a central theme in conservation biology and landscape management (Saunders 1991, Sadeh 2012, Oliver and Morecroft 2014. SSP analyses have greatly helped to assess the consequences of habitat fragmentation and degradation, and to evaluate the role of the matrix between habitat patches (Blixt et al. 2015, Davies-Mostert et al. 2015, Miller et al. 2015, Hamer 2018. However, habitat availability and fragmentation are not the only processes influencing the dynamics of SSP. First, there is increasing recognition of a major role of density-dependent factors, which can influence multiple parameters such as survival, emigration and subpopulation dynamics (Cayuela et al. 2019). Furthermore, all the environmental factors affecting fitness in subpopulations can impact the dynamics of nearby patches; the SSP approach can thus be successfully applied to a broad range of environmental stressors. For instance, the SSP concepts can help to understand and mitigate the impacts of climate change on populations (Oliver and Morecroft 2014) and could also allow better assessing the impact of invasive alien species (IAS).
Invasive species are a frequent disturbance factor for patchily-distributed native species, as they can affect the fitness of native species both directly (e.g. through predation and competition) and indirectly (e.g. habitat modifications) (Nentwig 2007, Vilà et al. 2010. IAS influencing the net reproductive rate within subpopulations are expected to ultimately affect the flow of dispersers that can emigrate to nearby populations (Hejda et al. 2009). Assessing the impact of IAS is a major task for ecologists and conservationists, still, the effects of alien species on the overall dynamics of SSP are rarely considered. In fact, even though several studies recognize that invasive species often exist within SSP, the impacts of IAS are generally measured (and consequently managed) only at the subpopulation level. Pond-breeding amphibians are excellent models to infer the effects of alien species on native SSP: their populations are structured into discrete patches where reproduction occurs (here, a pond = a patch) and connected by the migration of post-metamorphic individuals (Marsh and Trenham 2001, Smith and Green 2005, Cayuela et al. 2019). This 'ponds-as-patches' view has some limitations because in some cases nearby ponds are connected by a very high rate of dispersal (Smith and Green 2005); nevertheless, the SSP concepts remain essential to understanding the dynamics of amphibian populations. The analysis of factors acting jointly at the patch and landscape scales might indeed provide key insights on the processes regulating population persistence and allows planning robust management strategies.
The aim of this study was to evaluate the impact of invasive species at the SSP level, assessing whether IAS affect the dynamics of subpopulations through their influence on potential source populations. To achieve this task, we evaluated the long-term effects of the invasive crayfish Procambarus clarkii on the threatened Italian agile frog Rana latastei, considering negative impacts at both the subpopulation and the network levels. Despite multiple evidence of heavy predation pressure of the crayfish on amphibian populations (Cruz et al. 2006, Ficetola et al. 2011a, only a part of the studies detected a clear negative relationship between crayfish presence and amphibian occurrence, and some even suggested long-term coexistence between some amphibian species and the crayfish (Bélouard et al. 2019). Most of these works are snapshot studies describing the distribution of the alien crayfish and native amphibians in a single year across patches. However, species distribution is the results of complex dynamics occurring in the long-term and at the landscape scale. Understanding the impact of IAS on the distribution of native amphibians thus requires long-term information on the abundance of both native and invasive species. Unfortunately, such data are rare (Strayer 2012, Ficetola et al. 2018, Denoël et al. 2019). Here we take advantage of long-term data on the threatened Italian agile frog to test whether the impact of IAS on subpopulations has effects on the whole SSP dynamics. If the impact of crayfish on population dynamics mostly occurs at the subpopulation level, we predict a negative relationship between subpopulation size (i.e. the number of breeding females in a wetland, hereby abundance or subpopulation size) and the crayfish presence at the patch level. Conversely, if crayfish impact occurs over the whole SSP, we predict that the subpopulation size is strongly determined by the crayfish presence in nearby wetlands. Our study provides evidence that the occurrence of IAS in a landscape may affect local population size even in not invaded patches. This highlights the importance of incorporating SSP dynamics in the assessment of the impact of invasive species.

Study system
We focused on the effects of the invasive American red swamp crayfish Procambarus clarkii, on the threatened Italian agile frog Rana latastei. The Italian agile frog is a brown frog, endemic to the Po river lowland (northern Italy) and adjacent areas. This frog is a forest specialist; it is an explosive breeder and each female produces a single egg mass (~1000 eggs per clutch; Ambrogio and Mezzadri 2018) during late winterearly spring. Clutches are globular and are deposited onto submerged woods and tree roots within 1 m of the water surface (Lanza et al. 2007), therefore, they are well correlated to the number of breeding females (Ficetola et al. 2010). Previous studies have demonstrated that Italian agile frogs live in complex networks of spatially-structured populations and that wetlands at distances < 1500 m have connected population dynamics De Bernardi 2004, Stankovic andPobolisaj 2013). After metamorphosis, froglets disperse into the landscape, sexual maturity occurs usually at one year and young adults can colonize new wetlands. Maximum lifespan is four years, and adults generally show high site fidelity (i.e. after the first breeding season they continue to breed in the same pond, Glandt 1986, Berven and Grudzen 1990, Bernini et al. 2000, Guarino et al. 2003, Smith and Green 2005. The average age of breeding females is 1.7 yr, and just 50% of females survive until two years of age (Guarino et al. 2003). Given the very strong annual mortality, it is unlikely that individuals skip one breeding season, as this would mean a high risk of dying without breeding. Therefore, the number of breeding females likely corresponds to adult female population size.
Since the early 2000s, wetlands of the Po lowland have been invaded by the American red swamp crawfish. This invasive crayfish exerts a heavy predation pressure on the larvae of R. latastei and other amphibians, strongly reducing their breeding success (Ficetola et al. 2011a). The crayfish invasion is favoured by the colonization of permanent, large ponds where it can form numerous and stable populations. Hence, only part of wetlands in the landscape is colonized by the crayfish, thus determining spatially heterogeneous impact across amphibian subpopulations (Siesa et al. 2011).

Study area, surveys and environmental parameters
From 2010 to 2018, we regularly surveyed a network of 58 freshwater sites (mostly ponds) in the northern Lambro River basin (Lombardy, NW Italy; Supplementary material Appendix 1 Fig. A1). The average distance between ponds was 1227 m (SD = 752). The first detection of the invasive crayfish in this network dates back to 2008. Forty sites were monitored every year, and all sites were monitored at least every second year (average ± SD 7.27 ± 0.9 monitoring seasons per site). Each monitoring season consisted of five surveys between the end of February and May. Three surveys were performed by night and two by day. During surveys, we counted the number of frog clutches; the maximum number of detected clutches was taken as an indicator of adult population size (Ficetola et al. 2010). Furthermore, we used visual encounter surveys to assess crayfish occurrence. Pervisit detection probability of the crayfish is very high during night surveys (> 90%), thus the number and timing of our surveys allowed to assess crayfish occurrence with high confidence (Manenti et al. 2019a). We did not perform an exhaustive assessment of the occurrence of native predators such as fish or dragonfly larvae, which requires specific techniques. Nevertheless, previous studies showed that fish occurrence does not have strong impact on the distribution of R. latastei breeding sites (Ficetola et al. 2011a).
We recorded two wetland features known to affect wetland use by the Italian agile frog for reproduction: surface area and sun exposure. Light regime strongly influences wetland suitability for tadpoles (Halverson et al. 2003), and R. latastei is a specialist of shaded wetlands (Ficetola et al. 2011a). During the central hours of the day, we measured the percentage of aquatic site surface covered by shade according to the following rank scale: 1) shade covering < 10% of the surface; 2) shade covering between 10 and 30% of the surface; 3) shade covering 30-60% of the surface; 4) shade covering 60-90% of the surface; 5) shade covering > 90% of the surface. Furthermore, adult frogs live in forest environments, and are strongly associated with wetlands with a high forest cover in a radius of ~400 m (Ficetola et al. , 2011a, thus we measured forest cover within 400 m from each pond, using the 2012 vector map of the Lombardy region (< http://www.cartografia.regione.lombardia.it >; ground resolution: 3 m). Finally, for each pond, we considered two variables representing the processes occurring in nearby populations, to take into account the factors acting at the SSP level: surrounding clutch incidence and surrounding crayfish incidence. First, to take into account factors affecting population dynamics via migration, we considered the spatial connections between each pond and all surrounding ponds (Moilanen and Nieminen 2002). The connectivity S of a patch i with clutches in nearby ponds (hereafter: surrounding clutch incidence) was calculated using an incidence function model (Moilanen and Nieminen 2002) where d is the distance between the focal patch i and each of the remaining patches j, N is the number of laid clutches observed during the previous year in the j-th pond and α = 1/1500 m, i.e. the distance at which subpopulations are known to be connected by migrants (Moilanen andNieminen 2002, Ficetola andDe Bernardi 2004). Second, we calculated the incidence of P. clarkii in surrounding wetlands, weighted by distance as in Eq. 1 (hereafter: surrounding crayfish incidence). For both clutch and crayfish incidence, we only considered wetlands where the frog bred at least once during the study period, to avoid unsuitable sites.

Statistical analyses
Statistical analyses were performed on all sites for which in at least one season we found frog egg-clutches (n = 25). We assessed the factors determining the dynamics of subpopulations along the study period across wetlands, taking into account processes occurring both at the subpopulation level and within the SSP network. In so doing, we fitted hierarchical N-mixture models to the abundances observed across subpopulations ( We modelled both patch-scale and SSP dynamics following a Ricker-logistic model, which is a model of population dynamics combining an exponential growth model with density-dependent effects (Hostetler and Chandler 2015). The abundance at the first year of sampling was considered to be only affected by wetland features and forest cover, assuming a negative binomial distribution. When modelling population dynamics, we considered that the expected abundance E(N i,t ) at locality i and year t can be represented by: where ρ is a fixed effect representing the intrinsic growth rate of the metapopulation and η represents the strength of density dependence (Hostetler andChandler 2015, Mathieu-Begne et al. 2019). var 1 , var 2… var n are independent variables potentially affecting population dynamics. In a first step (patch-level analysis) we considered as independent variables crayfish occurrence in the patch at year t − 1, and the three variables describing the environmental features of breeding sites and the surrounding landscape: wetland surface area, sun exposure, and forest cover within 400 m from the wetland. In a second step (SSP-level analysis), we considered the same variables of the patch-level analysis, plus two variables representing the role of processes acting at the SSP level: the number of clutches laid in nearby wetlands (surrounding clutch incidence), and the occurrence of P. clarkii in surrounding wetlands (surrounding crayfish incidence). For both processes acting at the subpopulation and SSP levels, we assumed that subpopulation size in a given year can be affected by crayfish occurrence and clutch abundance during the previous year, because R. latastei generally reaches sexual maturity one year after metamorphosis and then shows site fidelity (Guarino et al. 2003). The density dependent component was only considered if population size was > 0. Given that detection probability is generally < 1, we also included the observation process in the model, using the Royle (2004) binomial model: where Y i,t is the number of frog clutches observed at pond i in year t, and p is the detection probability of frogs (Hostetler and Chandler 2015). We also included two random effects representing stochastic differences between ponds and years, drawn from normal distributions. For all the parameters, we used noninformative priors. Specifically, for regression coefficients, growth rate and density dependent effects we used a vague normal prior with mean = 0 and variance = 100; for the standard deviations of the random effects we used a half-Cauchy prior (Gelman 2006). A uniform prior bounded between 0 and 1 was used for detection probability, and a uniform distribution, bounded between 0 and 50, was used as prior for the overdispersion parameter of the negative binomial distribution. See Mathieu-Begne et al. (2019) for additional details on the metapopulation model. Before running analyses, pond area, surrounding clutch and surrounding crayfish incidence were log-transformed to reduce skewness and improve normality. We used a Bayesian approach to obtain the posterior distribution of parameters, implemented in JAGS 4.3 (Plummer 2003) using the R package R2jags (Yu-Sung and Masanao 2015). We run three MCMC chains, each with 300 000 iterations and a thinning rate of 100; we discarded the first 30 000 iterations as a burnin. Gelman and Rubin diagnostics were ≤ 1.1 for all parameters, indicating convergence (Kéry 2010). We used deviance information criterion (DIC) to compare the patch-and SSP-level models. DIC is a measure of model complexity and fit that is particularly useful for the comparison of Bayesian models; models with lowest DIC values are generally assumed to have a better performance (Spiegelhalter et al. 2002). All data and scripts are available in Supplementary material Appendix 2.

Results
During the nine years of study, we detected at least one breeding R. latastei in 25 sites out of the 58 monitored ones (average ± SD number of clutches for each site/year = 6.2 ± 25.8).
Despite some fluctuations, the total abundance across all the wetlands did not significantly change through the study period (Spearman's correlation between total number of clutches and year: r s = 0.367, p = 0.336; Supplementary material Appendix 1 Fig. A2).
Considering the whole study area, the invasive crayfish increased its occurrence from 13 sites in 2010, to 21 sites in 2018. Within the 25 sites where we detected breeding R. latastei, the crayfish increased its occupancy form 10 sites in 2010 to 16 sites in 2018 (Supplementary material Appendix 1 Fig. A1). When we analysed subpopulation size only considering the patch variables, frog abundance was negatively related to sun exposure and pond area, while we did not detect any significant impact of the crayfish (Table 1).
When we modelled subpopulation size as a function of both patch and SSP variables, hierarchical models showed that frog abundance was jointly determined by habitat features and by processes acting at the SSP scale. First, the analysis confirmed the importance of habitat features (Table 2; Fig. 1), with the highest abundances in the most shaded wetlands; for pond area, 95% intervals overlapped zero. At the subpopulation scale, the number of clutches laid in a year was unrelated to the occurrence of crayfish at that site. Conversely, at the SSP level, we found a strong and negative relationship between frog abundance in a given patch and crayfish incidence in the surrounding landscape, and a positive relationship between frog abundance and the abundance of clutches in surrounding patches (Table 2; Fig. 1). The SSP-level model showed a much lower DIC value than the patch-level one (SSP model: DIC = 1566.6; patch model: DIC = 1786.5).

Discussion
Long-term analyses of population networks are a powerful approach to understand the consequences of IAS, as they allow to tease apart processes occurring at the patch and at the whole SSP scales. When subpopulations are connected, the impact of IAS can extend beyond the simple patch-level effects, affecting subpopulation size even in not invaded sites. Our study focused on the Italian agile frog, a threatened amphibian that has been extensively studied before the onset of the crayfish invasion. This frog typically breeds in permanent ponds within lowland forests (Ambrogio and Mezzadri 2018), and such wetlands are also suitable for the invasive crayfish (Siesa et al. 2011, Manenti et al. 2014, Manenti et al. 2019a. It is widely recognized that the invasive crayfish heavily preys on amphibian larvae, alters wetland food webs and can be responsible for local populations extinctions (Cruz et al. 2006, Ficetola et al. 2011a, Grillas et al. 2018), still our patch-level analyses did not detect a significant effect of crayfish occurrence on adult frog abundance within focal ponds (Table 1). This result mirrors previous patch-level analyses, which failed at detecting negative relationships between crayfish occurrence and local occupancy of several native amphibians (Ficetola et al. 2011a, Rodriguez-Perez et al. 2014, Bélouard et al. 2019. Such a lack of negative relationships suggested possible coexistence between the invasive species and some amphibians (Rodriguez-Perez et al. 2014, Bélouard et al. 2019) and the same conclusion could have been drawn from our study, had we limited analyses at the patch level.
However, the demography of a large number of species (including many amphibians) is best described considering the spatial connections between (sub)populations, as immigration and emigration are major determinants of local abundance. Our SSP-level analysis revealed that the subpopulation size was jointly affected by processes acting at both the local and the network level (Fig. 1). Local effects were important also in our study system in term of habitat features. In fact, abundance was higher in shaded ponds within forest, confirming the key role of habitat variables for Italian agile frogs (Ficetola et al. , 2011a. Furthermore, local abundance was strongly related to the abundance of clutches in surrounding patches (Fig. 1), highlighting the importance of immigration from nearby ponds. Immigration is possible  for both froglets and for adults breeding in nearby ponds. Nevertheless, amphibians often show fidelity to breeding sites, as adults generally continue to breed in the same wetland year after year (Wells 2007). Moreover, the short lifespan of the Italian agile frog suggests a primary role of the immigration by froglets, and highlights that breeding sites can be attractive for immigrant frogs metamorphosed in surrounding patches, regardless to red swamp crayfish occurrence. Adult Italian agile frogs are probably unable to detect the crayfish cues when selecting breeding sites, therefore they continue to breed in invaded wetlands (Ficetola et al. 2011a). Explosive breeders in particular, such as agile frogs, can be unable to detect predatory crayfish because females remain in breeding ponds only for a few hours. Furthermore, even if the crayfish preys on frog tadpoles and can reduce breeding success (Cruz et al. 2008, Ficetola et al. 2011a, Milligan et al. 2017, some pond features may allow the coexistence between amphibians and the crayfish (Bélouard et al. 2019). Aquatic vegetation and other shelters often favour amphibian larvae escape from predators (Teplitsky et al. 2005, Manenti et al. 2017, and the survival of tadpoles is generally density-dependent (Vonesh and De la Cruz 2002), thus some tadpoles can be able to reach metamorphosis even in invaded ponds.
When we shift to the SSP perspective, the clear relationship between frog abundance and crayfish occurrence in nearby populations suggests that IAS can impact the whole SSP, with effects that can remain unnoticed when not considering the complex network within which each subpopulation persists. After metamorphosis, froglets move to their terrestrial habitat (forests) where they remain until sexual maturity and then move toward potential breeding sites, which are sometimes different from birth sites. In our study, the negative effect of the crayfish was only detectable when considering its occurrence over the whole SSP. This confirms the importance of recruitment from nearby ponds for the dynamics of subpopulations and shows that factors influencing reproduction at potential source sites have strong effects on all the connected patches (Seward et al. 2019). When dealing with a predator IAS invading new habitats, one may hypothesise that the new predator just replaces the other sources of mortality that already occurred for native species. The red swamp crayfish effectively reduces the density and diversity of aquatic insects preying on tadpoles (Ficetola et al. 2012). However, the direct predation pressure exerted by the crayfish on amphibian tadpoles is stronger than the one exerted by the whole community of native predators in uninvaded ponds, suggesting that the reduction of predatory insects is not enough to compensate for the negative impacts of the crayfish (Ficetola et al. 2012).
The state-space model also allowed to estimate important parameters, such as density dependence values and detection probability of frogs (Table 1). The density dependence parameter was very close to zero, suggesting very weak density dependence. However, it should be remarked that the total population size of this threatened frog was very small considering the size of study area (roughly 500 ha), therefore frog density likely is much lower than the carrying capacity commonly observed in amphibians (Santini et al. 2018) generally found densities of 100-> 1000 individuals ha -1 ). The estimated detection probability of frogs was low, suggesting that a substantial number of clutches remained undetected, still the observed detection probability matched the finding of other studies using N-mixture models and capture-markrecapture to estimate the detection probability of amphibians (Lampo et al. 2012, Brown et al. 2013, Courtois et al. 2013). Contrary to our expectations, we found a weak negative relationship between frog abundance and forest cover. However, in our dataset the largest pond were outside the forest, and the negative correlation between pond area and forest cover (r = −0.63) probably complicates the identifiability of the effect of these environmental features.
The correct identification of the scale at which IAS impact native species has strong management consequences. Had we only considered the patch level, we might conclude that the management of alien crayfish is less important than the conservation of vegetation cover around breeding sites, since the negative role of the invasive species was only evident when considering the processes affecting the whole network of subpopulations. If we focus on subpopulations processes, our results indeed stress the importance of maintaining shaded ponds within forests. However, the emigration/immigration processes are pivotal for the long-term persistence of population networks (Cayuela et al. 2019). Current recruitment at non-invaded ponds seems to be enough to maintain the SSP (Supplementary material Appendix 1 Fig. A2), but extinction thresholds often occur in the relationship between habitat loss and species extinction at the landscape scale. Assessing in a further study migration or colonisation/extinction rates in the SSP framework of our study area could help to determine if such thresholds really occur. Emigration from source patches can allow species persistence until a given level of habitat loss, but the risk of extinction abruptly rises when the frequency of unsuitable landscape reaches a threshold value (Fahrig 2003, Ficetola andDenoel 2009). IAS could cause a similar process: in our study system, the frequency of invaded ponds is steadily growing (Supplementary material Appendix 1 Fig. A2), and the whole SSP could collapse if the number of source ponds will decrease too much. Multiscale management, with actions at both the pond-and the landscape-level, is essential in these circumstances. Even if the eradication of the crayfish from the whole landscape is not feasible, actions targeted on a subset of ponds (e.g. manipulation of the hydroperiod; Ficetola et al. 2011b) could limit crayfish frequency, thus allowing the persistence of a sufficient number of source ponds and the maintenance of the whole SPP. Furthermore, the structure of landscapes in which patches are embedded is pivotal, and actions of IAS eradication and control may only have limited effectiveness if directed toward single patches, without considering their spatial connections (Morel-Journel et al. 2019). More studies are needed to understand if also alien crayfish population densities are important, beyond just presence/absence information; assessing crayfish density can be challenging and to compare multiple wetlands at the SSP amphibian level, specific protocols need to be tested.
There is growing evidence that predators can have an impact even beyond the boundaries of their distribution, affecting the immigration/emigration dynamics of their prey even in predator-free patches (Resetarits 2018, Trekels andVanschoenwinkel 2019). For instance, some amphibians avoid ponds invaded by alien predators, and this has cascading consequences at the landscape scale, with modification of emigrations dynamics, strong changes in local densities and fitness decrease due to density-dependence (Winandy et al. 2017). Furthermore, experiments on mosquitos showed that oviposition decreases in safe patches near risky patches, while it increases in risky patches if they are close to predator-free patches (Trekels and Vanschoenwinkel 2019). The relationships between factors acting at multiple (local, regional) spatial scales can be particularly complex, underlying the importance of a multi-scale approach to understand spatial variation of species distribution and fitness (Denoël andLehmann 2006, Trekels andVanschoenwinkel 2019). Our results clearly show that the effects of IAS may affect population dynamics even in not invaded patches through their degradation of the whole SPP. Patch-scale assessments of the impact of invasive species are thus insufficient to understand their actual effects, and predicting the long-term interplay between invasive and native populations would require a landscape-level approach which also considers the complexity of spatial interactions.