Reconciling cities with nature: Identifying local Blue-Green Infrastructure interventions for regional biodiversity enhancement

Increasing urbanization degrades quantity, quality, and the functionality of spatial cohesion of natural areas essential to biodiversity and ecosystem functioning worldwide. The uncontrolled pace of building activity and the erosion of blue (i


Introduction
Urban sprawl and agricultural intensification are major drivers of habitat degradation, loss and fragmentation of ecosystems resulting in declines of biodiversity and ecosystem services worldwide (Butchart et al., 2010;Horvath et al., 2019).Habitat fragmentation increasingly challenges wildlife persistence and movement in human-dominated landscapes, due to urbanization and the erosion of blue (i.e., aquatic)green (i.e., terrestrial) spaces to accommodate a growing human population.Consequently, key biological processes such as breeding, dispersal, migration and resource utilization are disrupted (Fletcher et al., 2018), ultimately reducing gene flow between populations with profound implications for genetic diversity and adaptive processes (Cayuela et al., 2018).
Protected areas (designated as 'ecological infrastructures' in Switzerland) are recognized conservation measures (Jongman et al., 2011) suitable to counteract such negative processes.More broadly, many of these are commonly referred to as Green infrastructure (Benedict and McMahon, 2012) or, increasingly, Blue-Green Infrastructure (BGI) -strategically planned, interconnected, nature-based solutions (encompassing 'green' vegetation and 'blue' water spaces and technologies) that not only provide biodiversity enhancement, but a range of other ecosystem services to cities and their hinterlands (Bozovic et al., 2017;Infrastructure, 2013;Oral et al., 2020).Yet, the protective effectiveness of the often spatially distributed areas may be limited in the long-term, as these areas are challenged by the steadily growing cities which act as regional barriers to the connectivity of natural habitats (Bolliger and Silbernagel, 2020;Nori et al., 2015), indispensable to maintaining biodiversity goals (Rodrigues et al., 2004).BGI help make cities less vulnerable to impacts of climate change such as preventing flooding, stormwater quality, drought, heat, air quality and quality of life for people (Kuller et al., 2017;Pauleit et al., 2017).Such instruments also offer interesting synergies for biodiversity and ecosystem functions (EEA, 2011;Bozovic et al., 2017;Oral et al., 2020).Yet, biodiversity enhancement aspects of BGI are frequently relegated as a secondary benefit among various ecosystem services in place of addressing more anthropocentric concerns such as water quality (Bach et al., 2020;Walsh et al., 2005;Zhang et al., 2020), urban heat (Gill et al., 2007;Wong et al., 2021) and stormwater protection (Zölch et al., 2017).This is in spite of their potential to create a wide variety of aquatic and terrestrial habitats.In particular, parks, cemeteries, ruderal areas, roadside rain gardens, canals, rivers and their riparian zones and even industrial infrastructures, such as wastewater treatment and recovery plants (that use natural treatment processes and storage ponds) can also make an important contribution, with the potential to form a spatially explicit, functional network of (semi-)natural blue and green habitats that maintain biodiversity and optimize ecosystems to support protected areas (EEA, 2011;Pauleit et al., 2017).
Urban-rural landscape planning research is increasingly focusing on finding strategies and tools that can better support the design of local measures while accounting for regional landscape and balancing conflicts between nature and humans (De Montis et al., 2016;Modica et al., 2021).Few examples, which do address biodiversity are fraught with spatial and taxonomic biases (McIntyre, 2000;Meyer et al., 2016) and a predominant emphasis on the 'green' component of BGI (Fontana et al., 2011).A prominent strategy and active consideration of biodiversity needs in BGI planning, to satisfactorily meet these needs, should: (1) identify actually occupied habitat ranges and secure their quality through the establishment of 'green' and 'blue' measures, and, (2) identify, protect or sustainably manage ecological networks i.e., systems of existing natural habitat patches interconnected through continuous corridors or "stepping stones" between existing landscape patches (De Montis et al., 2016;Modica et al., 2021;Moilanen et al., 2005).Enabling this is, however, difficult as we lack significant knowledge on the permeability of rural-urban landscapes, for vertebrates living at the interface of aquatic and terrestrial ecosystems and thus, are unable to select and spatially prioritise BGI measures across such large human-dominated landscapes.Several studies have used expert opinions and field data on few individual species to provide knowledge about within-city scale factors (Braaker et al., 2014;Finch et al., 2020), while a regional-scale generalization remained difficult due to the large amount of data required.Furthermore, considering the cost implications and the lack of equipment to be able to exhaustively tag small and threatened species safely and ethically with global positioning system (GPS) technology, we highlight the need to develop non-invasive methods for examining conservation issues surrounding landscape fragmentation at multiple spatial scales as well as approaches to spatially plan BGI coherently and, more strategically, in urban environments (Kuller et al., 2018;Schmidt et al., 2019).
Monitoring and species distribution models (SDMs) as well as connectivity models are promising tools to support conservation planning regionally (Clauzel and Godet, 2020;Dondina et al., 2020;Khosravi et al., 2018), yet in practice, barriers to connectivity, as well as conservation actions, frequently operate at much smaller and artificial (e.g.cities) spatial scales, which are seldom considered.For example, decisions must be made about the probable effect of a blue-green space erosion because of a redevelopment plan on the ability of a local wildlife population to access parts of its habitat and, hence, what (if any) mitigation is required.
Globally driven to the brink of extinction by urbanization, Amphibians, may particularly benefit from and serve as representative model system in supporting the implementation of BGI at strategic locations (e.g., Churko et al. 2020aChurko et al. , 2020b) ) enhancing their movability at broader spatial scales (Arntzen et al., 2017;Hamer and McDonnell, 2008;Holzer, 2014;Stuart et al., 2004).Amphibians spend their early life-stages in streams, ponds and wetlands (Cayuela et al., 2020), while juveniles often disperse from their natal pond to other spawning sites and adults alternate seasonally between terrestrial and wetland areas (Cayuela et al., 2020).Hence, the preservation and restoration of the "blue" and "green" attributes in core habitats (e.g., spawning sites), as well as the permeability of the landscape between different spawning sites, and between spawning and non-spawning habitats (e.g., feeding and hibernation sites), are essential determinants of their persistence in highly fragmented landscapes (Parris, 2006), making them a suitable model organism for BGI planning to enhance biodiversity in an urbanizing world.
The aim of this paper is to develop and demonstrate a rigorous and integrated approach to model opportunities for BGI to enhance biodiversity in human-dominated landscapes.Herein, we provide a multiscale, multi-species integrated assessment framework for BGI planning as baseline information to enhance blue-green biodiversity in humandominated landscapes.To this end, we combined biological and spatially, highly resolved urban-rural land-cover data to (1) determine species-specific suitable habitats and hotspots of amphibian biodiversity along with an evaluation of their current state of overlap with existing protected areas (i.e., ecological infrastructure), (2) identify major landscape elements acting as continuous ecological corridors and determine the role of the urban form in facilitating or impeding movability.Lastly, we combine (1) and (2) to identify interventions for BGI planning within urban settings to support connectivity and enhance broader regional landscape connectivity, essential to the maintenance of essential biodiversity processes and ecosystem functioning in humandominated landscapes.

Study area
The study area (3133 km 2 ) comprises the administrative Cantons of Aargau (AG) and Zürich (ZH) in the Swiss lowlands (264 m-1212 m a.s.l.), characterized by distinct, human-dominated settlement patterns.AG exhibits rather distributed small settlements whereas ZH comprises several major urban centers, including the cities of Zürich (approx.415,000 inhabitants) and Winterthur (approx.110,000 inhabitants; Fig. 1).The study area includes 42.9% of agricultural land (predominantly meadows and pastures), 31.1% of forests, 22.2% of urban settlements (of which: 12.2% are buildings and other impervious landcovers, and approx.10% are urban green spaces), 3.8% of water bodies, and <1% of landcovers, such as raised bogs and fens, rocks or rubble-sandy grounds.
Various protected blue-and green-areas contributing to ecological infrastructures (Swiss term) include the regional parks (e.g., the Jura Park of Canton Aargau, the Sihlwald of Canton Zürich), Ramsar and Emerald sites, alluvial zones-dry meadows and pastures-, raised bogs and fens-, and amphibian breeding regions-of national importance (Ryser et al., 2002), as well as other natural reserves (Pro Natura sites; see Fig. 1).Green spaces within settlements encompass approximately 10% of the study area and included vegetated areas, such as gardens, house turnarounds, public parks, cemeteries, and sports facilities.

Overview of integrated assessment framework
The concept of our multi-scale, multi-species integrated framework for the identification of BGI opportunities in support of amphibian movability in human-dominated landscapes is shown in Fig. 2. The framework integrates ensemble models (EM) of species distributions (Di Cola et al., 2017;Thuiller et al., 2009) and circuit theory (Anantharaman et al., 2019), results of which can then be spatially analyzed to, firstly, understand actual/potential biodiversity hotspots and coverage by existing protected areas referred to as ecological infrastructure (Swiss term), secondly, identify landscape elements acting as essential ecological corridors through the study area, that may be protected/restored or further enhanced through BGI measures.The framework encompassed a data assembly for ten species and amphibian whole-life cycle environmental predictors (step1), habitat suitability modelling and identification of actual/potential species occupancy ranges (step 2), modelling their regional ecological network (step 3), and identifying opportunities for strategic BGI consideration in urban-rural settings that could further enhance their regional connectivity (step 4).

Amphibian selection and summary of datasets
Amphibian species selection.We investigated 19 candidate amphibian species (Fig. S1) and selected 10 species using a multidimensional trait analysis (Villeger et al., 2008).We compared candidates against six traits according to Trochet et al. (2014) to cover the overall ecological and life-history trait diversity related to dispersal for amphibians (see Supplementary Material -SM Table S1 and Table S2 for more information on relationship with dispersal): (1) Adult body size: snout-to-vent length in adults [mm]; (2) Metamorphosis size: snout-to-vent length of juveniles before metamorphosis [mm]; (3) Parental care: simplified binary variable; (4) Egg/offspring number: number of eggs/offspring laid; (5) Juvenile diet composition: three categories, Her (specialised herbivores), Herb-Detr (species feeding on vegetation and detritus), and higher-TL (species having carnivorous, insectivorous, molluscivorous and cannibalistic diet composition); and (6) Displacement mode: three categories, one-medium for species displacing either in terrestrial or aquatic habitats, multiple for aquatic and terrestrial displacement, multiple vertical for aquatic, terrestrial and vertical (i.e., climbers) displacement.
Species were also selected to attain a good representation of the phylogenetic diversity of amphibians (Pyron and Wiens, 2011) by selecting at least one representative species per family (Fig. S3).Precisely, the ten species retained in the analysis were Alytes obstetricans, Bombina variegata, Bufo bufo, Hyla arborea, Ichtyosaura alpestris, Lissotriton helveticus, Lissotriton vulgaris, Rana temporaria, Salamandra salamandra, and Triturus cristatus.
Amphibian occurrence data.Observations of amphibians were gathered from the Coordination Center for the Protection of Amphibians and Reptiles of Switzerland (http://www.karch.ch).To cover the most recent amphibian distribution in the study area, we included observations from 2017 to 2019 (Table S1), ranging from 108 observations for one of the species and up to 813 observations for another (averaging around 350 observations per species).All the presence records were georeferenced to a 10 m × 10 m grid system.Observations were disaggregated to 10m using the ecospat.disaggregation()function from the R package ecospat (Di Cola et al., 2017).Amphibian whole-life cycle environmental predictors.To quantify the ecological niche encompassing the whole life cycle of each amphibian species, 13 predictors across six environmental variable types (prepared at 10m resolution) describing all amphibian life stages were considered (Table S3): (1) Topographicslope [ • ]: derived from the Swiss Topographical Landscape model (swissTLM 3D 2015) at 2 m resolution (elevation uncertainty ± 0.5m) and, subsequently, aggregated to a median value; (2) Hydrologic: two predictors including nearest distance to water [m] (all water-related objects and rasterized lotic, lentic and reedbeds from the regional cadastral maps -AV: Daten der Amtliche Vermessung) (Kanton Aargau, 2020;Kanton Zürich, 2020) and runoff coefficient [C] (between 0 and 1 relating amount of runoff to amount of precipitation, mapped AV: Daten der Amtliche Vermessung data (Kanton Aargau, 2020; Kanton Zürich, 2020) against literature values for different land coverssee Table S4); (3) Edaphicsoil moisture variability: gathered from the spatial modelling of Ecological Indicator Values in Switzerland (Descombes et al., 2020) and resampled from 100m to 10m; (4) Vegetation: we calculated the median and standard deviation of the Normalized Difference Vegetation Index (NDVI) from all available cloud free Sentinel-2 pixels available between April and October for the years 2016-2019 (Drusch et al., 2012)  To avoid multicollinearity and to maximize the accuracy and predictive power of our models (Guisan and Zimmermann, 2000), we calculated the variance inflation factor (VIF, Quinn and Keough, 2002) for all pairs of environmental predictors (Table S5) through the vistep() function from the R package usdm (Naimi 2015).Since we found low (VIF <3) collinearity for all our predictors (Zuur et al., 2010), we retained them in the subsequent analyses.

Modelling species distributions using ensemble modelling
The assessment of the amphibian-species distribution relied on an ensemble modelling approach (Thuiller et al., 2016).By fitting several species distribution model (SDM) techniques using various types of regression models, EM decreases the uncertainty associated with using a single modeling technique and, therefore, increases the accuracy of model predictions (Bolliger et al., 2017;Kozak et al., 2017).To ensure a more robust output we considered a wide range of modeling techniques: generalized additive model (GAM); generalized boosting model (GBM); generalized linear model (GLM); random forest (RF), and Maximum Entropy model (Maxent).Multi-species potential distribution modelling in the complete geographical space was performed using biomod2 platform (Thuiller et al., 2016).Datasets included species data (presence and pseudo-absences) and the 13 environmental predictors.Furthermore, for all species, we used 1000 pseudo-absences and a random pseudo-absence strategy, except for the most representative species, R. temporaria and I. alpestris for which we generated 2000 and 2,500, respectively.This ensured a more balanced "presence to pseudo-absence" ratio.To avoid the effect of sampling bias, we generated five random pseudo-absences datasets for each species across the study area following biomod2 recommendations.
We ran five replications for each model where 70% of occurrence points were used as a training set, while the remaining 30% was used for the model evaluation (Thuiller et al., 2016).The relative contribution of each environmental predictor was estimated by five permutations.All models were evaluated for their accuracy using the areas under the curve (AUC) of the receiver operating characteristic (ROC) and the true skill statistic (TSS; Allouche et al. (2006)).To integrate model outputs and develop the EM predictions, we used a weighted-averaging approach through which each single distribution model was weighted according to its predictive accuracy.The performance of the ensemble models was evaluated using the AUC, TSS and the Boyce index (Hirzel et al., 2006), a presence-only evaluator.Projections were performed under current conditions (2017-2019) and selection of the optimal threshold for binarization using the 10th percentile of training presence points (or "mpa" threshold -"minimum predicted area") (Preau et al., 2020).This threshold, computed with the mpa.ecospat function in the "ecospat" R package (Di Cola et al., 2017), allows the removal of locations where habitat suitability values from the EMs outputs are lower than the suitability values of the bottom 10% of the species occurrence (Pearson et al., 2007), which could be affected by errors in the data or represent sink populations (Preau et al., 2020).
For each species, from the binarized EMs output, we extracted the actual occupied habitat patches in the geographic space (i.e., referred to as occupancy areas) based on their respective occurrence data (2017-2019) using the ecospat.occupiedpatch() function implemented in the "ecospat" package (Di Cola et al., 2017).The difference between the binarized EMs outputs and the occupancy areas resulted in the inoccupancy areas.We also calculated the percentage of occupied over the total potential suitable habitats (% occupied cells/total potential cells) to assess the status of occupancy of suitable habitats and inform on regional species conservation prioritization.Similarly, we calculated the % of occupied-and potential-habitat patches which overlapped or not with ecological infrastructures to assess the actual and potential degree of protection of habitat ranges.
Finally, by overlaying the binarized actual and potential distribution maps for all the species we also assessed the actual and potential amphibian biodiversity "hotspots" (calculated as Species Richness, SR) across the study area (Wiens, 2007).

Connectivity models and identification of key ecological indicators
For all the species, we performed inter-breeding site connectivity analyses to identify species-specific and multispecies ecological corridors.We used the Circuitscape model's Julia implementation (version 1.5) (Anantharaman et al., 2019).We set up an all-to-one mode for raster format, run on a High-Performance Cluster.The software requires a resistance map layer and a focal point file containing the geographic coordinates of each source location.The resistance map represents the degree with which a landscape feature facilitates or hinders species-specific movement across the landscape.For each species we aggregated the weighted average EM suitability map to 30m to avoid computational issues given the size of the region.We then converted it to a resistance map using an exponential decay function following (Duflot et al., 2018) using the species-specific mpa threshold previously computed (Preau et al., 2020).The function assigns a resistance value of 1000 when habitat suitability equals zero, and a value of one when habitat suitability is greater than or equal to the binarization threshold (see Table S6).The exponential shape strengthens the barrier effect of less favorable areas.As focal nodes (i.e., source and ground nodes representing start and end points for movement) we used all amphibian spawning sites of national importance, to predict the maximum inter-spawning sites movement potential of amphibians across the study area.All resulting connectivity maps were normalized (between 0 and 1) using the raster.transformation()function implemented in the "spatia-lEco" R package (Evans and Ram, 2021).Ecologically, the modelled spatially explicit surfaces represented potential amphibian connectivity and movement across the landscape.The amount of movement can also be referred to as current, i.e., high current indicates high movement potential, whereas low current refers to low movement potential.
To verify the validity of our connectivity maps and evaluate the congruency between simulated and observed movement patterns of amphibians in human-dominated landscapes we intersected current maps (i.e., degree of potential amphibian movement) with an independent dataset of sites of observed amphibian migration, available from the Coordination Center for the Protection of Amphibians and Reptiles of Switzerland (http://www.karch.ch),and collected between 2017 and 2019.Such migration sites consist of yearly monitored and assisted (temporarily or permanently) locations where amphibians are known to cross a road to move from their summer-winter habitats to the spawning sites.Using this data set, we computed the percentages of migrations sites, which overlapped within 1st to the 99th percentiles of the normalized cumulative current maps (i.e., all potential routes of species movements across all spawning sites in the study area) for each species and at the multi-species level.

Analysis of BGI opportunities in human-dominated landscapes
Using the 90th percentile of the normalized cumulative current maps (i.e., the highest current values), we binarized the connectivity maps (0 or 1) to quantify to which degree human settlements overlap with individual species amphibian dispersal pathways across the study area.Hence, we quantified the amount of the binarized connectivity pathways that intersect with urban "green" (e.g., private and public parks and gardens) and "grey" spaces (e.g., roads, other impervious surfaces).Based on this, we identified locations where such overlaps occur and thereby, where opportunities for BGI may be achieved either through a transformation of already existing green-(rather than a physical transformation of green space that would for example involve management transformations of vegetated landcover) or from impermeable greyareas to support amphibian connectivity.

Overall predicted amphibian habitat suitability
We used all 3537 amphibian observations collected between 2017 and 2019 to build models for our ten amphibian species.Amphibian occurrence points included both, sites insides and outside of existing ecological infrastructure, providing a comprehensive representation of species distribution patterns in the study area.Comparisons of the different modelling types (i.e., GLM, GBM, GAM, RF, MAXENT) included in the EMs revealed high correlation of the predicted suitable habitats (Table S7).All modelling approaches performed reasonably well at predicting occurrence, indicating that environmental variables were good predictors of potential habitats of individual amphibian species (Table 1).Among modelling techniques, RF's predictive power outperformed other algorithms (Table S7).However, based on the spatial agreement between predictions of the different algorithms and the robustness of the EM (Table 1), we consider the results from the EM through the rest of the integrated modelling framework.
Predicted area currently occupied indicate major amphibian biodiversity hotspots along the rivers Rhine, Aare, Reuss and Sihl, with scattered hotspots also found in peri urban areas (Fig. 3a).We also identified substantial multispecies (re-)colonization potential of suitable habitats beyond the current occupancy hotspots throughout the whole study area (Fig. 3b).The average importance of the environmental predictors among the models indicated that species select for specific habitat types within the study area (Table 1, Fig. S4).However, the amount of urbanization, proximity to water bodies and to forest were the most important predictors of amphibian's potential suitable habitats (Fig. S5).Major gaps between amphibian richness hotspots and existing ecological infrastructure occurred along the following watercourses: the Limmat, the Rhine, the Glatt, the Töss and the eastern part of the Thur, as well as in the southeastern region of the study area (Canton Zurich; Fig. 3).

Spatial variability of species-specific suitable habitats
The species-specific projection of the EMs allowed the prediction of the individual suitable habitats as well as the extraction of 2017-2019 occupied habitat patches important for regional species conservation prioritization.More precisely, the proportion of truly occupied and potentially suitable habitats differed widely between the investigated species.Three species were found to occupy less than 20% of their potential suitable habitat ranges, namely L. vulgaris (15.5%), A. obstetricans (16.2%) and H. arborea (17.1%;Table S8).For example, A. obstetricans was mainly observed in broad areas of the regional Jura Park (northwest of the study area) and in proximity to the major rivers Aare and Reuss (Fig. 4a), with more than 63% of its occupied habitats lacking protection by ecological infrastructures.L. vulgaris and H. arborea showed similar potential distributions and habitat occupation ranges.Both species occurred along the Reuss river and in the northern periurban area of the city of Zürich (H.arborea shown as an example in Fig. 4b).For all three species, on average, half of their occupied areas (i.e., mean ± SD: 9.4% ± 0.8%) and the majority of suitable, yet unoccupied, habitats (i.e., mean ± SD: 71% ± 2.4%) were not covered by existing ecological infrastructure (purple areas in Fig. 4a and b, Table S8).
Intermediate ranges of occupied, potentially suitable habitat (20-50% of the modelled habitat) were found for T. cristatus (23.56%),B. variegata (32.84%, shown in Fig. 4c), L. helveticus (42.13%) and I. alpestris (48.19%).B. variegata showed an extensive proportion of occupied habitat over the total potentially modelled habitat (24.86%) occurring outside from any ecological infrastructure (red areas in Fig. 4c).In addition, accounting for 57.94%, several clusters of suitable yet apparently unoccupied nor protected habitats appeared in the southern regions of the study area (purple areas in Fig. 4c).

Predicted connectivity hotspots
We identified four landscape elements contributing to the regional connectivity of amphibians between spawning areas: (1) forest edges, (2) soils with variable moisture, (3) wet forest habitats and (4) riparian areas (variable importance, Fig. 5).Additionally, areas with the highest (i.e., >90th quantile) values of cumulative current for all amphibians overlapped with 11.3% of urban "green" spaces (e.g., private gardens, parks, sport courses) and 3.4% of "grey" spaces (i.e., impervious surfaces).This suggests that movement may also occur within the urban

Table 1
Performance of the Ensemble Models (EM) built with five algorithms (GLM, GBM, GAM, RF, and MAXENT) for the ten amphibian species.The three most important variables in the development of ensemble model potential habitats are indicated.The AUC derived from signal detection theory is the area under the receiver operating characteristic curve (ROC).Evaluation criteria for the AUC statistic are as follows: excellent (0.90-1.00), very good (0.8-0.9), good (0.7-0.8), fair (0.6-0.7), and poor (0.5-0.6).
The TSS considers the true positive and negative prediction rates (Allouche et al. 2006).TSS ranges from − 1 to +1 with scores >0 denoting SDMs performing better than random.

The Boyce index measures the deviation between the model predictions and the random distribution of the observed presences across the range of prediction values and is analogous to a correlation value, varying between -1 (counter-prediction
) and 1 (perfect prediction), with 0 indicating random predictions.NDVI: Normalized Difference Vegetation Index, med: median, sd: standard deviation.setting.More precisely, B. bufo (16.8%) and S. salamandra (18.9%) were the species predicted to move the most through artificial landscapes (Table S8).We performed a visual comparison of the multispecies predicted connectivity map between spawning sites with an independent dataset of locations where amphibian movements during the migration period are assisted by human interventions for the same temporal extent (2017)(2018)(2019).This revealed high spatial congruency (Fig. S6) between predicted and observed amphibian movements, supporting the wholelife cycle predicted ecological corridor pathways obtained in our analyses.

Clusters of ecological corridors
To generalize results, we evaluated pairwise Spearman's correlations between the species-specific predicted connectivity maps.This resulted in four major movement patterns adopted by amphibians, allowing us to group the different species (Fig. 6, Fig. S7 and correlation results in Fig. S8).Two species presented distinct connectivity maps, namely A. obstetricans and S. salamandra.For the former, the most important explanatory variables in predicting suitable habitats and movement habitats included the standard deviation of NDVI, distance to forest and to water (Fig. 6a).For the latter, predominant variables predicting preferred habitats were distance to water and forest (Fig. 6b).Furthermore, we identified two groups of species moving similarly through the landscape: (1) H. arborea, L. vulgaris and T. cristatus, for which the variability of the soil moisture, the distance to the aquatic medium, permeability of the soil estimated through the runoff C and vegetation (NDVI median) were important predictors (Fig. 6c) and (2) a large number of species such as B. bufo, B. variegata, I. alpestris, L. helveticus and R. temporaria.For the latter group, distance to water and forest remained as predominant predictors, but were complemented by additional vegetation-related predictors (Fig. 6d).

Identifying Blue-Green Infrastructure opportunities
To enhance amphibian connectivity between the occupied areas and/or towards new suitable habitats, we localized where BGI could be best planned.This may be achieved either through a transformation of already existing green areas (rather than a physical transformation of green space it would for example involve management transformations of vegetated landcover) or from impermeable grey areas to support amphibian connectivity.Among the species with the narrowest occupancy and most distinct movement pathways, we found that suitable locations for BGI implementation for A. obstetricans could be planned within the urban setting of Döttingen (AG, Fig. 7a) in proximity of the Surb river ("green transition"), which could specifically support the connectivity between two major occupancy areas of the latter species (i.e., in northwest and northeast of the river Aare) and along the Aare river from a transformation of an existing gravel pit area ("grey transition").Additionally, a major pinch point of connectivity was found in proximity of the city Winterthur, specifically along the river Töss in the Eschenbergwald forested area, which emerged as another unique connectivity pathway.BGIs in this area may enhance the connectivity between the occupied habitats in the northern part of ZH to extensive suitable areas currently only poorly occupied in the southeast (Fig. S9).
For S. Salamandra, also presenting distinct movement pathways, but a comparatively widespread occupation of potential suitable habitats, two potential major pinch-point zones (narrow high current) were identified within the largest urban settlement of the study area, potentially enhancing connectivity of the east and west occupancy areas of Lake Zürich, but implying substantial grey area transformations (Fig. 7b).The first consisted of riparian areas around the Sihl river, while the second was a stepping-stone pathway from the eastern peri-urban forest bordering the city (i.e., Aldisberg region), through scattered green urban spaces (e.g., in Hottigen and Hirslanden municipalities) (re) connecting to the Limmat river.Further pinch-points occurring in comparatively rural landscapes were identified, connecting minor and scattered occupancy areas such as for example in the northwest, where two minor streams (Mohlinbach and Magdenerbach) emerged as ecological corridors through the settlements of Möhlin and Rheinfelden (Fig. S10).
For H. arborea, T. cristatus and L. vulgaris, all species occurred at four major locations: (1) in scattered habitats of the peri-urban area of Zürich, (2) in AG along the rivers Aare and Reuss, and (3) in ZH along the river Thur (Fig. 8a).Our results suggest that opportunities for BGIs exist, for example: along the Limmat river, particularly within wide existing green spaces encompassing sports fields, such as golf courses (see Fig. 8a), along the Glatt river (Fig. S11) and, in the riparian zones of the Aare within the urban settlements of Brugg (Fig. S11).
In the peri-urban areas of Zürich, we found opportunities for BGIs along the Glatt and Chriesbach rivers for more generalist species such as B. variegata (Fig. 8b).We also visualized substantial predicted movement along forest edges, which extensively border agricultural land.Such patterns were also found, for example, in the least urbanized regions of the study area, such as in the poorly populated rural settlement of Böbikon (AG) (Fig. S12).

Evaluation of the integrated assessment framework
Our integrated assessment framework aimed at mapping options for BGI using amphibian habitat suitability hotspots and movement corridors with the aim of rigorously assessing opportunities for BGI in human-dominated landscapes in the Swiss lowlands.Such opportunities for BGI are required because urban sprawl, intensive agriculture and a dense traffic infrastructure increasingly fragment the natural habitat (Jaeger et al., 2008).Our framework offers a promising strategy not only for regional conservation actions, on priority areas and species to efficiently direct conservation efforts, encourages the involvement of urban planners to strategically design BGI at the local scale, support broader landscape connectivity (Hodgson et al., 2011) essential for biodiversity protection and enhancement.
Our outcomes were built on predictive models based on species occurrence, habitat suitability and landscape resistance to predict interbreeding site connectivity.A challenge such a framework is to build realistic and robust models to guide relevant conservation actions.This point benefitted from the integration of species observation data (i.e., occurrence and migration) as well as highly resolved whole-life cycle environmental predictors as input into the modeling (Bani et al., 2015).Our EM provided reliable predictive maps of amphibian species habitats and the environmental predictors' effects aligned well with amphibian ecology (Clauzel and Godet, 2020;Preau et al., 2020).Furthermore, our connectivity modelling method relied on resistance surfaces Fig. 5. Potential multispecies amphibian movement between spawning sites (light green polygons).The colour gradient for corridors represents predicted connectivity (i.e., the normalized cumulative current) between core patches from weak (white) to strong (dark red).Forests are shown in dark green, ecological infrastructures in dark grey, meadows and pastureland (yellow) and human settlements (building polygons) are shown in black.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)transformed from the EM (Duflot et al., 2018) and considered multiple paths simultaneously (Anantharaman et al., 2019), which are recognized to be adapted to model connectivity in complex habitat matrices such as in urban settings (Braaker et al., 2014) and an advancement (McRae et al., 2008) over conventional least-cost-path methods (e.g., Kong et al. (2010)).Such an approach also offers a quantitative alternative to replace the subjective use of expert-informed resistance maps (Churko et al., 2020a), which were found to be suboptimal and underestimated compared to empirical approaches (Keeley et al., 2016;Sawyer et al., 2011;Stevenson-Holt et al., 2014).Also, Zeller et al. (2018) concluded that resistance maps derived from habitat suitability models may be an alternative to costly, time-consuming and more invasive methods such as radio-tracking and mark-recapture protocols.

Implications for landscape management and conservation policy
Overall, EM outcomes within this study area were able to pinpoint actual and potential amphibian biodiversity hotspots that, upon further inquiry, currently lack conservation measures.For the considered landscape and temporal extent (2017-2019), we also identified three species of highest regional conservation priority (2017-2019), namely the endangered L. vulgaris, A. obstetricans and H. arborea, all of which exhibited the narrowest and most scattered occupancy ranges (see Fig. 4).For these species, we found the broad lack of protected area coverage in core habitats.
For effective conservation of these species, we recommend an integrated landscape management approach at the regional scale, supported by the urban/local measures that especially allow for genetic interchange of individuals between natural patches along ecological corridors (Le Lay et al., 2015;Rannap et al., 2009;Sievers et al., 2019).In this regard, four landscape elements of fundamental importance for amphibian movements across the landscape could be identified through this study (Fig. 5).In agreement with Clauzel and Godet (2020), the first comprised forest edges, found to be essential pathways for most species (n = 6) and particularly vital for the A. obstetricans which declined by 50% in the past quarter century in Switzerland (Cruickshank et al., 2016).Furthermore, these habitats were generally adjacent to agricultural lands, implying possibilities for co-management plans between forestry and agriculture for an effective preservation and restoration of forest edges.
A second ecological corridor emerged from habitats characterized by highly permeable grounds with variable soil moisture, especially in rural areas supporting the promotion of wet arable land practices for the enhancement of wetland-associated biodiversity (Churko et al., 2020b;Knutson et al., 2004).These practices could increase the prevalence of seasonal temporary wetland habitats for species such H. arborea, L. vulgaris and T. cristatus and be especially significant given the additional continued pressure of ongoing climate change (Preau et al., 2020).
Thirdly, for more widespread species such as, S. Salamandra, we found that, although the majority of occupancy areas are not designated protected areas, their core habitats and corridors are (and guaranteed in the long-term), as long as forest (particular wet-forest habitats with natural hydrology and structural complexity) remain preserved (Manenti et al., 2009), well connected (Fig. 3d) and exempt of predators (e.g., fish) (Baumgartner et al., 1999).This is key for the conservation of the S. Salamandra, characterized by low dispersal capability (Schulte et al., 2007) and recently reported to be in continuous abundance decline in the present study area (Bänziger et al., 2017).
Finally, for more ubiquitous species (e.g., B. bufo, R. temporaria) our outcomes highlight the importance of preserving riparian zones, forest edges and promoting the establishment of BGI at strategic locations in urban and peri-urban regions to increase the permeability and availability of "stepping stone" habitats in highly artificial landscapes (Grant et al., 2019).Indeed, the major contribution of our approach is the identification of BGI planning opportunities within urban settings, that may enhance broader connectivity, for multiple species while accounting for their individual ecological requirements.

Study limitations and further work
Although comprehensive, our approach is not exempt from limitations.Our habitat suitability models could have benefitted from presence-absence instead of presence-only data, as they contain more ecological information (Elith et al., 2006).Similarly, the inclusion of abundance data may also have provided important information on species-environment relationships (Howard et al., 2014), especially for Fig. 7. Potential amphibian movement represented by connectivity maps for (a) A. obstetricans and (b) S. salamandra overlaid with occupied habitats (dark green) as well as with predicted unoccupied suitable habitats (transparent gold).Identification of BGI through overlaps of binarized connectivity maps (90th quantile of normalized current values) with "green" (blue color) and "grey" (red color) urban areas.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)species such as S. Salamandra, which was predicted to be rather widespread, yet shown to experience population declines (Bänziger et al., 2017).Although valuable, abundance data need to be collected using comparable sampling designs in terms of time and number of visits, and therefore is seldom available at broader spatial scales, even between rather small administrative boundaries as between AG and ZH.
As input for the connectivity modelling, we used the species-specific transformed map from the EM outputs.However, this transformation assumes that animal movement behavior reflects similar factors to habitat selection, which may not always be the case (Zeller et al., 2012).Nevertheless, the inclusion of predictors related to movement behavior such as, distance to nearest water body, distance to forest and road, enabled us to build whole life cycle suitability models, which account for different processes such as dispersal, encompassing from the aquatic to dispersal phase (Bani et al., 2015), and migration from the aquatic to terrestrial feeding or sheltering sites (Cayuela et al., 2020).A final major limitation in our analysis is that we considered maximal dispersal potential by using all spawning sites as focal points in the connectivity analysis.This might overestimate the connectivity of species in certain areas of the landscape, especially if they are restricted to narrow and scattered occupancy areas (see e.g., H. arborea).However, given that the surface of the study area is rather limited (3133 km 2 ) in comparison to other case studies (Bani et al., 2015;Clauzel and Godet, 2020) and to the amphibians species distributions across Europe (Araújo et al., 2006), and that the understanding of dispersal distance remains relatively poor with high intra-specific variability observed (Cayuela et al., 2020), this approach allowed us to better understand the maximal connectivity potential across the entire study area, which can become relevant should new habitats form along current corridors or long-term urban planning require identification of future suitable habitats for active protection measures.
The establishment of BGI in urban settings may further support regional connectivity through the urban form.To date, the authors do not know of any comprehensive study on BGI for biodiversity enhancement, which has used such an approach.Previous work has frequently relied on key indicators or coarse approaches (e.g., Meerow and Newell (2017); Nguyen et al. (2021)), but has not integrated biodiversity hotspots and connectivity modelling for this purpose.Enhancing our work in future could integrate new tools such as the Urban Biophysical Environments and Technologies Simulator (Urban-BEATS) (Bach et al., 2020) to support the planning and design of ecological and liveable cities and other multi-functional aspects of BGI Fig. 8. Potential amphibian movement represented by connectivity maps for (a) Hyla arborea and (b) Bombina variegata overlaid with occupied habitats (dark green) as well as with predicted unoccupied suitable habitats (transparent gold).Identification of BGI through overlaps of binarized connectivity maps (90th quantile of normalized current values) with "green" (blue color) and "grey" (red color) urban areas.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)(Lafortezza et al., 2013;Meerow and Newell, 2017) such as urban heat mitigation (Koc et al., 2018) and broader ecosystem services (Kuller et al., 2017).For example, structurally diverse green roofs were found to support a broad species richness, including rare species (e.g., arthropods, beetles and spiders among others) through the provisioning of a variety of niches (Knapp et al., 2019).Further, constructed wetlands, including animal-aided designs (e.g., barrier-free design of shores provides a diverse aquatic and terrestrial vegetation, temporal drying-out of parts of the structure among other aspects) found to benefit wetland-associated biodiversity (Knapp et al., 2019).

Conclusions
This study demonstrated an integrated assessment approach to improving biodiversity assessment in human-dominated landscapes, aimed particularly at the strategic planning of Blue-Green Infrastructure.We used biological and highly resolved urban-rural land-cover data, ensemble models of habitat suitability and connectivity models based on circuit theory to identify single-and multi-species habitats of importance for resource utilization and dispersal from regional, to highly artificial and complex local landscapes.Overall, this study demonstrated that the integration of habitat suitability models and connectivity analyses is a valid strategy to identify core areas and ecological corridors and to inform regional and local conservation management plans in order to reduce fragmentation in humandominated landscapes and to locate opportunities for BGI.Particularly, we identified that, within the present study area of the Swiss lowlands, forest edges, wet forest habitats, wet arable lands and riparian zones are key ecological corridors and support multiscale preservation and restoration measures of these habitats.We also identified candidate locations for the design of BGI that may significantly contribute to enhancing the regional functional network of (semi-) natural habitats to maintain biodiversity and optimize ecosystem functions.For effective, long-term conservation of wetland-associated biodiversity our work encourages increased collaboration between regional and local stakeholders, from federal environmental agencies to local urban planners to preserve and enhance the quality of the spatial cohesion of ecological networks.

FUNDING sources
We thank the ETH Board for funding through the Blue-Green Biodiversity (BGB) Initiative (BGB 2020).

Fig. 1 .
Fig. 1.Study area.Cadastral map of the AG and ZH along with maps of the major settlements (names in red), water bodies (names in blue) and ecological infrastructure (grey shade).The extremes of the human-dominated landscape are represented by the densest (Zürich city) and the least (Böbikon) populated settlements in the study area.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) and the median vegetation height [m] gathered from Ginzler and Hobi (2015); (5) Land-use derived: includes three predictors, grassland density [grassland/cell] -data based on Pazúr et al. (2022), total grasslands area evaluated using a spatial weight matrix and a 250m neighborhood radius around each cell; urbanization [buildings/cell] -average building density for a 250m neighborhood around each cell based on the AV: Daten der Amtliche Vermessung (Kanton Aargau, 2020; Kanton Zürich, 2020); and traffic intensity [number of daily vehicles/cell] on roads -Movement-ecology related: includes three predictors of proximitynearest forest distance [m]; nearest rock-gravel and sandy area (RGS) distance [m] and nearest road distance [m] (considering any road type) -all calculated based on the AV: Daten der Amtliche Vermessung maps (Kanton Aargau, 2020; Kanton Zürich, 2020).

Fig. 2 .
Fig. 2. Overview of multi-scale, multi-species integrated modelling framework to assess biodiversity enhancement and Blue Green Infrastructure opportunities in urban landscapes.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 3 .
Fig. 3. a) distribution of predicted area occupied by amphibians, b) predicted amphibian distributionsboth refer to species richness (SR) and illustrate overlap with existing ecological infrastructure in the study area.

Fig. 4 .
Fig. 4. Comparative overview of species-specific occupancy and inoccupancy areas according to the binarized EMs and the occurrence data, along with overlaps with existing ecological infrastructure.a) A. obstetricans and b) H. arborea represent narrowly-distributed species, c) B. variegata is a species with intermediate occupation range and d) S. salamandra is an example for a widespread species.Colors indicate whether the occupied (red) or suitable but currently unoccupied habitat (purple) overlaps or not with an existing ecological infrastructure.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6 .
Fig. 6.Boxplots of variable importance describing the modelled suitable habitats of the ecological corridors clusters (right).Panels with boxplots (a-d) indicate averaged scores assigned to the environmental variables used to develop the EM predictions and the subsequent resistance maps (from SDM outputs) for the four connectivity clusters.Icons represent the species found in the clusters (a: Alytes obstetricans; b: Salamandra salamandra; c: Hyla arborea, Lissotriton vulgaris, and Triturus cristatus; d: Bombina variegata, Lissotriton helveticus, Rana temporaria, Bufo bufo and Ichthyosaura alpestris).
Ensemble model weighted mean by TSS merged algo-run-data