Landscape genetics identifies barriers to Natterjack toad metapopulation dispersal

Habitat fragmentation and loss reduce population size and connectivity, which imperils populations. Functional connectivity is key for species persistence in human-modified landscapes. To inform species conservation management, we investigated spatial genetic structure, gene flow and inferred dispersal between twelve breeding sites of the Natterjack toad (Bufo calamita); regionally Red-Listed as Endangered in Ireland. Spatial genetic structure was determined using both Bayesian and non-Bayesian clustering analysis of 13 polymorphic microsatellite loci genotyping 247 individuals. We tested the influence of geographic distance, climate, habitat, geographical features, and anthropogenic pressure on pairwise genetic distances between breeding sites using Isolation-by-distance and Isolation-by-resistance based on least-cost path and circuit theory models of functional connectivity. There was clear spatial structuring with genetic distances increasing with geographic distance. Gene flow was best explained by Isolation-by-resistance models with coniferous forestry plantations, bog, marsh, moor and heath, scrub, anthropogenic presence (Human Influence Index) and rivers (riparian density) identified as habitats with high resistance to gene flow while metapopulation connectivity was enhanced by coastal habitats (beaches, sand dunes and salt marshes) and coastal grassland. Despite substantial declines in census numbers over the past 15 years and its regional status as Endangered, the Natterjack toad population in Ireland retains high genetic diversity. If declines continue, maintaining habitat connectivity to prevent genetic erosion by management of coastal grasslands, pond construction and assisted migration through translocation will be increasingly important.


Introduction
Habitat degradation and fragmentation are key drivers of the current global biodiversity crisis (Sala et al. 2000;Foley et al. 2005). Habitat fragmentation, through habitat loss and habitat patch isolation, changes animal behaviour including dispersal (Janin et al. 2012) and movement patterns (Poessel et al. 2014), mortality rates (Pinto et al. 2018), population growth (Bascompte et al. 2002), population structure (Haag et al. 2010) and population viability (Newman et al. 2013). Regions experiencing habitat loss have a greater proportion of species in decline than regions of intact habitat (Donovan and Flather 2000), with the degree of landscape fragmentation key to population persistence (Clobert et al. 2001;Hanski 2001). Unconstrained animal movements are important in foraging, breeding, dispersal, (re)colonization and essential for responding to environmental change (Zeller et al. 2012). These factors are particularly important in metapopulations where a species distribution is disjunct and separated into discrete, dispersed populations. Factors limiting dispersal rates increase isolation and are thus likely to limit gene flow and increase genetic differentiation between populations, potentially accelerating inbreeding and genetic drift (Fahrig 2002;Johst et al. 2002).
The global amphibian crisis, whilst predominately driven by diseases (fungal Chytridiomycosis), has been exacerbated by habitat destruction and degradation such that 41% of species are threatened with extinction (IUCN 2020). Amphibians have poor dispersal abilities and are particularly vulnerable to habitat fragmentation (Bowne and Bowers 2004;Graeter et al. 2008) due to small effective population size (Funk et al. 1999), high site fidelity (Joly et al. 2003) and a complex terrestrial-freshwater life cycle which necessitates two distinct environments narrowing their habitat tolerances (Houlahan and Findlay 2003). Lack of connectivity between breeding sites and suitable terrestrial habitats can lead to high mortality during breeding migrations and low recruitment of dispersing individuals (Bulger 2003;Janin et al. 2012). Amphibian metapopulations are highly dynamic and susceptible to local extirpations and turnover, hence habitat connectivity is crucial for recolonization (Brown and Kodric-Brown 1977) and population persistence (Trenham et al. 2003). Pond-breeding amphibians are notably sensitive to isolation due to habitat fragmentation between breeding sites, which is a key issue in their conservation (Cushman 2006). Understanding dispersal, population genetic structure and the impact of landscape features on habitat connectivity are crucial for developing targeted and species-specific conservation management approaches at the landscape level. Moreover, such information can be used to predict the impact of proposed land use changes and/or infrastructural developments helping shape mitigation strategies such as the creation of dispersal corridors (Storfer 2007).
Monitoring dispersal through direct observations can be costly, difficult and time consuming (Broquet and Petit 2009), where it is possible at all. Recent development in geospatial information technology, molecular biology and the resolution of several statistical problems in spatial genetics such as nonindependence among samples (Prunier et al. 2015;Peterman 2018) facilitate the indirect study of movements using gene flow and Global Information Systems (GIS) (Balkenhol et al. 2009). Landscape genetics integrates population genetics and landscape ecology to detect genetic discontinuities to quantify the effects of geographical distance and landscape permeability on metapopulation genetic structure (Manel et al. 2003;Storfer et al. 2007). As an alternative to a null model where no spatial genetic structure is assumed, an isolation-by-distance (IBD) approach assumes that observed genetic differences between populations are a function of geographic distance i.e. isolation alone. IBD is expected for populations inhabiting continuous habitats where the main factor contributing to the genetic differentiation is the distance between populations and the dispersal capabilities of the studied species.
Habitat fragmentation creates dispersal barriers limiting gene flow. An isolation-by-resistance (IBR) approach assumes that observed genetic differences between populations are a function of landscape resistance to dispersing individuals (McRae 2006). Resistance surfaces represent the cost to a dispersing individual in crossing a landscape where low resistance values represent ease of movement and high resistance values represent restricted movement due to the presence of barriers along a permeability gradient (Zeller et al. 2012). IBD indicates a lack of population isolation due to habitat fragmentation while IBR indicates population isolation due to limited dispersal explained by habitat fragmentation and loss of connectivity (Kobayashi and Sota 2019).
This study aimed to investigate the role of geographic distance, climate, habitat, geographical features, and anthropogenic pressure in determining the spatial genetic structure of a species of conservation concern in Ireland, the Natterjack toad (Epidalea calamita (King et al. 2011). The species is widely distributed over much of southwest and central Europe, inhabiting environments with a wide range of climatic and habitat conditions including human-modified landscapes like agricultural land and quarries (Reyne et al. 2021a). However, the species is often associated with open habitats on sandy substrates with shallow, ephemeral ponds for breeding (Beebee 1979). Radiotracking data suggests that adult Natterjack toads can travel up to 4 km outside the breeding season with a maximum daily movement of approximately 300 m . Movement across terrain is facilitated by bare ground or where vegetation is low and spare like sand dunes and grasslands with moderate grazing, while habitats with dense vegetation like forest plantations and highly fragmented landscapes through urban development restrict movement (Cobert 2006;Maes et al. 2019).
The species is highly range-restricted in Ireland, found only in Co Kerry. Ireland has lost most of its freshwater ponds (Reid et al. 2012), with the loss of suitable breeding sites identified as the single most important driver of Natterjack toad population declines (Reyne et al. 2021b). The Natterjack toad is Red-listed as Endangered with a 23% decline in fecundity (number of eggs strings deposited annually) in Ireland over a 14-year period. However, at some breeding sites (Roscullen and Dooks Golf Course DGC) egg string production declined by over 90% since 2004. Fecundity has increased only at three breeding sites in the North of the species distribution (Magharees, Castlegregory Golf Course CGC and Lough Gill), but this is insufficient to offset the overall species decline in Ireland (Reyne et al 2021b). To halt the Natterjack toad decline, the National Parks and Wildlife Service (NPWS) initiated a 'Pond Creation Scheme' in 2008 installing over 100 artificial ponds on agricultural grasslands throughout the Natterjack toad range. Only artificial ponds in close proximity to natural breeding pools had been colonised by 2018 (Reyne et al. 2019(Reyne et al. , 2021a. Consequently, the NPWS initiated a Natterjack toad 'Head-start and Reintroduction Programme' in 2016 by collecting egg strings/tadpoles annually to raise toadlets in captivity before release back into the wild as part of assisted migration and translocation to newly created ponds. The success of the programme is yet to be determined, as breeding will occur 4 to 5 years post-release when toads reach sexual maturity and return to the ponds for breeding. The objectives of this study were to: i) estimate genetic structure and quantify pairwise genetic distances between remaining breeding sites, ii) quantify climatic and habitat landscape variability, and iii) relate genetic distance to a) geographic distance and b) landscape dispersal resistance explicitly testing isolationby-distance (IBD) and isolation-by-resistance (IBR) models. We hypothesised that due to the Natterjack toad's affiliation with coastal habitats in Ireland, inland climate and habitat factors would limit dispersal and explain landscape genetic structure better than geographic distance alone. Our results inform ongoing species conservation management particularly in respect to pond creation and assisted migration through translocation.

Study area and sampling
We investigated the gene flow between all known breeding sites (twelve discrete patches with clusters of breeding ponds within 500 m) (Fig. 1). Distances as the crow-flies between breeding sites ranged from 900 m to 24 km. Breeding ponds were in a variety of habitats: sand dunes, wet grasslands, golf courses and agricultural land mainly used for sheep grazing. Field surveys were conducted during the Natterjack toad breeding season (April to July) in 2017. We collected samples for DNA extraction consisting of welldeveloped egg strings, tadpoles or tissue samples from adult toads that had died naturally. Eggs and tadpoles were collected from forty natural and artificial breeding ponds, with only a single egg genotyped from each egg string, to avoid pseudoreplication by analysing siblings. The Natterjack toad 'Head-start and Reintroduction Programme' had no impact on this study, as at the time of data collection the translocated toadlets were too young to breed. Samples collected in the field were stored in 100% ethanol at ambient temperature until extraction. All samples were collected under license to capture protected wildlife animals (Licence No. C098/2016) issued by the National Parks & Wildlife Service (NPWS), Department of Culture, Heritage and the Gaeltacht, Republic of Ireland.

DNA extraction and genotyping
DNA extraction and genotyping followed Reyne et al. (2022). Genomic DNA was extracted following a high salt extraction protocol. Individuals were genotyped using thirteen highly polymorphic fluorescently labelled microsatellites. A summary of laboratory methods is given in Supplementary Information (Online Resource 1). We checked for null alleles with the R package PopGenReport v3.0.0 (Adamack and Gruber 2014). Approximately 10% of the samples (25 samples) were randomly selected and genotyped three times to calculate the rate of genotyping errors (allelic dropout and false alleles) using PEDANT v 1.0 software (Johnson 2007).

Genetic diversity
Analyses of genetic diversity were performed at the breeding site level. We estimated the number of alleles per locus, rarefied allelic richness, observed (H O ) and expected heterozygosity (H E ), and fixation index (Fis) using the hierfstat package (Goudet 2005) in R 3.6.2 (R Core Team 2019).

Population structure
We performed two clustering methods to infer population genetic structure: Bayesian and non-Bayesian approaches. We visualised the Natterjack toad genetic diversity in Ireland and estimated the number of clusters using STRU CTU RE v2.3.4. (Pritchard et al. 2010). K values between 1 and 10 were tested by 10 independent runs of STRU CTU RE with a burn-in of 250,000 iterations and run-length of 100,000 MCMC interactions for each run. STRU CTU RE output was visualised using the R package ggplot2 (Wickham 2016). Principal Coordinate Analysis (PCoA) is a multivariate technique that is not dependent on Hardy Weinberg or linkage equilibrium and allows major patterns within a multivariate dataset to be identified based on algorithms developed by Orloci (1978). We used two measures of genetic distances to build PCoA biplots (Gower 1966), F ST (Weir and Cockerham 1984) and corrected for sampling size G" ST (Meirmans and Hedrick 2011) calculated in FSTAT v2.9.4. (Goudet 1995) and GenAlEx 6.5 respectively. To investigate cryptic genetic patterns as a result of isolation, we performed a spatial Principal Components Analysis (sPCA; following Jombart et al. 2008). The method maximized the variance in individual allele frequencies while simultaneously accounting for spatial autocorrelation estimated using Moran's I (Moran 1948(Moran , 1950 with 9,999 randomized Monte-Carlo permutations to test for differences in allelic frequencies between neighbouring breeding sites. Data were analysed using two multivariate statistical tests (global and local tests) to assess if there were significant patterns at either scale. Presence of a global pattern indicates positive spatial autocorrelation i.e. neighbouring breeding sites tend to be similar, while presence of local structure indicates negative spatial autocorrelation i.e. neighbouring sites are dissimilar. Genetic structure was analysed hierarchically, initially using the entire dataset then using sample subsets based on genetic clusters detected by PCoA analysis. Data were analysed using the package ade4 for multivariate analysis (Dray et al. 2007), spdep for spatial autocorrelation (Bivand 2013(Bivand , 2018 and adegenet for sPCA and multivariate tests (Jombart 2008) in R 3.6.2 (R Core Team 2019).

Environmental parameters
The Irish Centre for High End Computing (ICHEC) provided five spatially explicit climate variables derived from the COSMO-CLM5 ensemble model (averaging five different climate models: CNRM-CM5, EC-EARTH, HadGEM2-ES, MIROC5 and MPI-ESM-LR) downscaled at a 4 km, the smallest available cell resolution for Ireland (Table 1, Figure S1). Variables were selected for their perceived relevance to Natterjack toad biology, specifically, temperature (air and surface) influencing activity levels and tadpole development and precipitation, humidity and surface evaporation influencing the permanency of breeding ponds. Climate variables (especially air and surface temperatures) accounted also for differences in the altitude across the study area.
A total of eleven land cover/habitat, geographical feature or anthropogenic pressure variables were captured (Table 1, Figure S1). Land cover/habitat was extracted from CORINE2018 (EEA 2018) and summarized at 4 km grid cell resolution to match the climate dataset. Individual CORINE land codes were aggregated and collapsed to derive simplified, ecologically relevant habitat classifications (Table 1). Simpson's Diversity Index was derived to quantify habitat diversity. Geographical features included distance to coast and riparian corridor density derived from an Ireland-specific GIS line shapefile of freshwater watercourses. Anthropogenic pressure was captured by the Human Influence Index (HII) downloaded from the Last of the Wild Project (Wildlife Conservation Society 2005) aggregating human population density, human land use and infrastructure (built-up areas, night-time lights), and human access (roads etc.). All surfaces were mapped at 4 km cell resolution to match the available climate data. Spatial data were extracted using ArcMap v 10.5 (ESRI, California, USA).

Landscape genetic analysis
We assessed three main models: (a) isolation-by-distance (IBD), where gene flow was a function of the distance between breeding sites; (b) isolation-by-resistance (IBR) where gene flow was a function of landscape resistance between breeding sites; and (c) a null model where gene flow was not characterised by spatial structuring between breeding sites. IBD analysis related pairwise genetic distances between breeding sites using linear mixed effects models with a maximum likelihood population effects (MLPE) parameterization (following Clarke et al. 2002) performed in the package ResistanceGA (Peterman 2018). Three measures of distance were used: (i) Euclidian distance which was the shortest distance between each pair of sites as the-crow-flies, (ii) geographic distance which was the distance between each pair of sites taking into account marine areas and including 3D topography i.e. slopes and (iii) historical distance which was geographic distance modified to account for a potential connection in the past during low tides and in presence of sand spits between the breeding sites at Inch sand dunes and those on the Iveragh Peninsula (Fig. 1c). Distances were calculated using ArcMap v10.8. Two measures of genetic distance were used (i) F ST and (ii) G" ST calculated as described  (322) and Inland marshes (411)  11 Open Bare rocks (332), Sparsely vegetated areas (333) and Burnt areas (334)  12 Coastal habitats Beaches, dunes, sand (331) and salt marshes (421)  above. The dependent variable in each MLPE model was a genetic distance matrix, breeding site was fitted as a Random Factor to account for non-independence of values in the pairwise matrix, and each measure of geographic distance was fixed as covariate in three separate models (one for each distance measure). IBR analysis related connectivity with climate, habitat, geographic features and anthropogenic pressure variables following the ResistanceGA approach (Peterman 2018;Peterman et al. 2014). This approach transforms landscape surfaces into resistance surfaces by providing the best fit to genetic data. Pairwise resistance distances were calculated between breeding sites by implementing genetic algorithms to optimise surface resistance based on pairwise genetic and resistance distances. Two IBR scenarios were tested: i) leastcost path and ii) circuit theory distances (following Kivimäki et al. 2014) using the costDistance and commuteDistance functions (Peterman 2018). Least-cost path distances estimated the optimal path with least resistance between two breeding sites, while circuit theory distance estimated commute-time distances performed in ResistanceGA equivalent to Circuitscape simultaneously considering all possible routes between the breeding sites (McRae 2006). IBR analyses were performed in two steps. The first step optimised the resistance of a single environmental surface by testing all 16 environmental parameters separately. Surface optimizations were performed twice to check for consistency. Model fit for each resistance surface was assessed based on corrected Akaike Information Criterion (AICc) values. The fit of the two measurements of genetic distances (F ST and G" ST ) was evaluated based on marginal R 2 values (Kimmig et al. 2020). The second step combined the most relevant environmental variables (resistance surfaces from singlepredictor ResistanceGA analysis with the greatest statistical support) into a multiple resistance surface. MLPE models were performed using the same structure as described above. Environmental variables were selected based on low AICc values from single surface optimisation. Analyses were performed using 1,000 bootstrap iterations to assess the relative support of each optimised surface and the robustness of model selection. Each surface optimisation was performed twice to check for convergence following recommendations by Peterman (2018).

Genetic diversity
We successfully genotyped 247 samples at 13 microsatellite loci. Sample size between breeding sites was uneven ranging between 7 and 43 samples ( Table 2). All microsatellite markers were polymorphic. Null alleles were observed for all loci but Bcalµ1. For these loci, null allele frequency estimates ranged between 0.01 and 0.13 (Table S2). The mean error rate, determined by repeat genotyping, in the 13 microsatellite loci was 5% for allelic dropout and 0.2% for false alleles (Table S2). The total number of alleles per breeding site ranged from 40 to 73 alleles with mean rarefied allelic richness ranging from 2.5 to 3.0. Expected heterozygosity was the lowest (H E = 0.52) for Roscullen and highest (H E = 0.60) for Yganavan, Inch and Quarry (Table 2).

Population structure
The Natterjack toad metapopulation in Ireland exhibited significant genetic structuring with clear differentiation between breeding sites (Table S3) with genetic structure reflecting clustering and proximity of breeding sites (Fig. 2). The analysis of 247 samples using structure algorithms suggested the best clustering solution was K = 2 (Delta K = 59.442; see Table S3). Clustering of individuals for K = 2 separated the north breeding sites (Magharees, CGC and Lough Gill) from the rest of the breeding sites in Ireland  Figure S3a). sPCA analysis based on the first and second principal component axes scores suggested three genetic clusters (Online Resource 2 Figure S3b). Hierarchical analysis of those three clusters identified additional structure. Positive spatial autocorrelation was detected among breeding sites at Inch sand dune and Iveragh Peninsula (global R = 0.013, p = 0.006) with three clusters: (i) Inch (ii) Dooks, Dooks Golf Club, Yganavan and Nambrackdarrig (iii) Quarry and Glenbeigh (Fig. 3c). No spatial autocorrelation was detected between the three breeding sites at the north of the Dingle Peninsula (global R = 0.032, p = 0.319; local R = 0.034, p = 0.108) or between the two breeding sites at Roscullen Island (global R = 0.0519, p = 0.215; local R = 0.045, p = 0.838).

Landscape genetics analysis
Genetic distance increased with geographic distance regardless of the genetic and distance measures used (Fig. 4). Mantel correlation coefficients whilst highest for Euclidian distance were similarly high for historical distance (0.809-0.792) suggesting existence of a possible past connection between Inch sand dune and the Iveragh Peninsula as current geographic distance had notably lower coefficients. Euclidian distance showed the best correlation with genetic distance (0.864-0.859) and was used in all subsequent analyses. Natterjack toad spatial genetic structure was not explained by IBD alone (Euclidean distance) but was also strongly influenced by all resistance surfaces using both leastcost path and circuit theory scenarios (Online Resource 2 Table S6 and Table S7). For the least-cost path scenario, marginal R 2 values were similar regardless of the genetic distance measure used (0.737-0.850) with relatively small differences in the predictive power between environmental resistance surfaces. Riparian density was consistently identified as the single most significant resistance surface to gene flow > 2 ΔAIC C units away from all other parameters (Online Resource 2 Table S5). For the circuit theory scenario, marginal R 2 values were also similar regardless of the genetic distance measure used (0.768-0.807), however, the two measures of genetic distance did not converge on the same model parameters. Euclidean distance was identified as the most significant predictor for gene flow when using F ST as the genetic distance measure though conifer ranked as the second-best model < 2 ΔAIC C units away from the single best model suggesting comparable predictive power. Distance to coast was identified as the single most significant predictor for gene flow > 2 ΔAIC C units away from all other parameters when using G" ST as the genetic distance measure (Online Resources Table S7). Data transformations were Fig. 2 Structure analysis of the Natterjack toad population in Ireland using microsatellite data. Each panel with graphs corresponds to a different number of genetic clusters (K). Individuals are shown with vertical lines, and the membership proportion is shown on the y-axis tested during single resistance surface optimizations (Fig. 5). The ranking of the sixteen single surface optimisation models for both scenarios and genetic distances remained identical between the two independent optimisation runs.
After single surface optimisation, the top four models were combined into a multiple resistance surface. None of the climate variables was identified as predictors of gene flow, thus were not selected for further analysis. Results of the multiple surface optimisation indicated that gene flow was better explained by a combination of habitats rather than by IBD alone (Euclidean distance) or by IBR using a single predictor (Table 3, Online Resources 2 Table S7). The difference between the IBD model (Euclidean distance) and the best multiple resistance surfaces IBR model was large regardless of the scenario or genetic distance measure (ΔAIC C > 119). Models using least-cost path performed considerably better in comparison to circuit theory distances based on models' AIC and identified more areas with high resistance representing restricted movement (Table 3, Fig. 6). Landscape genetic analysis of multiple resistance surfaces using both IBR scenarios (least-cost path and circuit theory) suggested that Natterjack toad gene flow between breeding sites in Ireland was positively influenced by coastal habitats (beaches, dunes, sand and salt marshes) and grassland whilst habitats with the greatest resistance to gene flow included coniferous forestry plantations, bog, marsh, moor and heath, scrub, anthropogenic presence (Human Influence Index) and rivers (riparian density) (Fig. 6).

Discussion
The Natterjack toad in Ireland, regionally Red-listed as Endangered, is not a single panmictic population. Our results suggest a high level of genetic differentiation and population structure with restricted dispersal and gene flow between One of the most important traits affecting population genetic structure is species dispersal capability, as this is essential in maintaining gene flow between breeding sites and populations. The Natterjack toad is a highly mobile amphibian with annual migrations to breeding sites up to 3-4 km , thus a high rate of gene flow might be expected between breeding sites within the toad's dispersal ability. However, significant genetic differentiation was found between breeding sites in Ireland separated by less than 3 km, so geographic distance and dispersal ability were not the principal factors limiting gene flow. This was confirmed by landscape genetic analysis, where resistance distances based on habitat explained Natterjack toad genetic structure and gene flow better than Euclidian distance alone. None of the climate variables were selected in the final models suggesting that restrictions to gene flow imposed by differences in climatic conditions may operate on a larger scale than those related to land cover and anthropogenic pressure (Reyne 2021a). Analysis suggested spatial genetic structuring with three main clusters separated by distance and unfavourable habitat: (i) Magharees, Castlegregory Golf Course (CGC) and Lough Gill; (ii) Killeen and Roscullen; (iii) Inch, Dooks, Dooks Golf Club (DGC), Nambrackdarrig, Yganavan, Quarry and Glenbeigh.
Results suggested a recent exchange of individuals between the three breeding sites at the north of the Dingle Peninsula (Magharees, Castlegregory Golf Course and Lough Gill). The main habitats in this area are sand dunes and grasslands both identified in the present study as facilitating dispersal and gene flow. Cox et al. (2017) demonstrated that dispersal of the Natterjack toad is facilitated by sand dunes and beaches, and the largest Natterjack toad population in Ireland is in the Magharees sand dune system, but breeding depends on the formation of ephemeral breeding pools in sand dune slacks (Reyne et al. 2019(Reyne et al. , 2021b. Access to permanent ponds at the nearby Castlegregory Golf Course between two measures of genetic distance (y-axes) and three measures of geographic distance (x-axes). Euclidian distance was the shortest distance between each pair of sites as the-crow-flies. Geographic distance was the distance between each pair of sites including 3D topography. Historical distance was geographic distance modified to account for the existence of a sand spit previously connecting the population at Inch sand dunes and those on the Iveragh Peninsula. Regression lines with 95% confidence intervals are based on fitted values of linear mixed-effects model with maximum likelihood population effects parameterization (MLPE) 1 3 (CGC) and suitable sites along the lakeshore of Lough Gill provide a valuable alternative in dry years, thus maintaining connectivity between the breeding sites is essential. Keeping sand dune and coastal grassland vegetation open in structure to support toad dispersal may require cattle or pony grazing and control of scrub encroachment, particularly of invasive sea buckthorn, Hippophae rhamnoides (Plassmann et al. 2010).
The three breeding sites at the north of the Dingle Peninsula were isolated from those of Castlemaine Harbour by the Slieve Mish Mountains dominated by peat bogs, marsh, heath and moor with patches of scrub and coniferous forestry plantations. Such habitats were identified by landscape genetic analyses as surfaces with high resistance to dispersal, which, along with the distance between breeding sites north and south of the Dingle peninsula, suggests the Slieve Mish Mountains are impassable. Peatland and their associated marshy wetlands are acidic, negatively affecting local distribution and abundance of amphibians (Freda 1986) whilst upland heath and moor tend to be drier habitats and at higher elevation both being suboptimal for amphibians. Frei et al. (2016) also found forest to negatively affect Natterjack toad gene flow and population size in Switzerland, yet woodland was preferred to pastures and agricultural fields by Natterjack toadlets in southern Belgium (Stevens et al. 2004(Stevens et al. , 2006a. Natterjack toad habitat selection varies notably across its range throughout Europe and regional idiosyncrasies in ecology remain poorly studied. Coniferous forestry in Ireland usually consists of tightly planted nonnative sitka spruce (Picea sitchensis) or similar tree species which produce dense stands of trees poorly penetrated by light, with little or no ground cover, few broadleaved species, and a highly acidic pine needle leaf litter supporting fewer invertebrates compared to native broad-leaved woodland (Pedley et al. 2014). Such conditions are likely to be avoided by toads. Additionally, roads extending east to west along the Dingle peninsula probably represented further barriers to gene flow relevant to the Human Influence Index in this landscape.
The genetic cluster at Killeen and Roscullen breed entirely in artificial ponds constructed on grassland < 200 m from the coast as part of the NPWS Pond Creation Scheme. Grasslands and proximity to the coast facilitated Natterjack toad dispersal and gene flow with the species known to inhabit and forage even on intensively managed agricultural lands most notably during summer Miaud and Sanuy 2005;Schweizer 2014). A further 20 artificial ponds that have been created in proximity (within 4 km) to Killeen and Roscullen have not yet been colonised naturally up to 10 years after their creation (Reyne et al. 2019(Reyne et al. , 2021b. Most adult toads have high breeding-site fidelity with dispersal usually by juveniles (Stevens et al. 2004). Toadlets may avoid agricultural areas (Stevens et al. 2004(Stevens et al. , 2006a, or suffer high rates of mortality due to its management practices there e.g. grass harvest by silage cutting which may represent an ecological trap for wildlife more generally (e.g. Reid et al. 2010). Thus, colonisation of new ponds at the limits of Natterjack dispersal capability may be slow. Other demographic parameters such as low local breeding success may also contribute to lower recruitment and dispersal. Certainly, Natterjack toadsat Killeen and Roscullen have exhibited substantial census size declines over the past 15 years (Reyne et al. 2019(Reyne et al. , 2021b. We propose including Killeen and Roscullen in the Natterjack toad NPWS 'Head-start and Reintroduction Programme' with occupied breeding sites used as a source population, with assisted migration and translocation to adjacent artificial ponds. This will reduce the distance between breeding ponds and increase the success of dispersal. More generally, any Natterjack toad translocations in Ireland should utilise the closest breeding sites as the source population to maintain genetic provenance and spatial genetic structuring; existing genetic diversity precludes the need for population admixture between genetic clusters. Our results suggest gene flow between Natterjack breeding site at Inch sand dunes and those of the Iveragh Peninsula. As these breeding sites are now separated by open sea, this gene flow presumably occurred historically when Table 3 Multiple surface optimisation for (a) least-cost path and (b) circuit theory distances Model evaluation metrics were produced using 1,000 bootstrap iterations: Avg.AIC C is an averaged AICc value; avg.weight is the averaged AIC C weights; avg.mR 2 is the averaged marginal R 2 dispersal between the sites might have been possible during low tides and in presence of sand spits since eroded by strong winds and storms (Orford et al. 1999). In this study, we estimated high genetic diversity at Inch sand dune. However, complete isolation, lack of long-term population data and dependence of breeding success on ephemeral ponds formation in sand dune slacks raise concerns for the future viability of this population. Close monitoring is required whilst ongoing inclusion in the NPWS 'Head-start and Reintroduction Programme', supported regionally by Dingle Oceanworld aquarium will ensure continued recruitment at the site even in dry years when all breeding sites evaporate before metamorphosis. Genetic connectivity was detected between breeding sites at Dooks, Dooks Golf Club (DGC), Nambrackdarrig and Yganavan where the main habitat was grassland interspersed by bog and scrub. Abandonment of agricultural grassland management results in rank vegetation and scrub encroachment which are threats to Natterjack toads in the area (Reyne et al. 2019(Reyne et al. , 2021b. Scrub was also identified as a barrier to dispersal thus local management should focus on preventing further deterioration in habitat quality. Construction of additional artificial ponds between existing breeding sites may be beneficial (Cox et al. 2017).
We detected potential reduced gene flow between Glenbeigh/Quarry and the rest of the Iveragh Peninsula breeding sites likely because of two landscape features: Glenbeigh town and the Caragh River and estuary. Increased traffic density due to the N70 road that crosses Glenbeigh town, part of the high traffic 'Ring of Kerry' tourist route and local urbanisation (captured in the Human Influence Index) may also present barriers to movement. Riparian density i.e. rivers and streams also had a negative impact on gene flow. Toads are poor swimmers and waterways such as the Caragh River and its associated saltwater tidal estuary are likely barriers to dispersal. Recent development in the area of Glenbeigh and Quarry breeding sites may increase population isolation, and whilst these sites had good genetic diversity, they also have a small number of breeding females (Reyne et al. 2019(Reyne et al. , 2021b and the highest fixation coeficientin Ireland raising concerns about future viability. Translocation of toadlets as part of the NPWS 'Head-start and Reintroduction Programme' sourced from Yganavan and Nambrackdarrig (breeding sites in the same genetic cluster with high numbers of egg strings and tadpoles) to Glenbeigh and Quarry may increase population numbers and mitigate any risk of future genetic erosion. Consideration should also be given to the creation of new artificial breeding ponds in this area.
While methods used in this study offer high potential to investigate the functional connectivity of the landscape, some technical aspects need careful consideration. A variety of methods exist to quantify genetic distances between populations, but at present, there is no consensus on which genetic matric to use to maximise model selection accuracy in landscape genetic studies. In this study, both genetic measurements (F ST and G" ST ) performed equally well with high marginal R 2 values likely as a result of a high degree of spatial structure. Most genetic matrices perform well when sample size and genetic structure are high (Shirk et al 2017). Thus, the use of several measures of genetic distance to evaluate model performance might be needed when the population genetic structure is low. Moreover, spatial resolution in landscape genetic study can be crucial for detecting the effect of landscape patterns on gene flow and assessing functional connectivity. Studies need to consider species ecology and habitat heterogeneity but also the availability of environmental data and computational time. We considered the 4 km grid cell size to be an adequate compromise considering the Natterjack toad dispersal capabilities and lack of high-quality climate data at finer resolution for the study area.
Understanding endangered amphibian movement and connectivity between breeding sites and metapopulations is key to their conservation. We show clear spatial structuring of the Natterjack toad in Ireland explained by Isolation-by-resistance with coniferous forestry plantations, bog, marsh, moor and heath, scrub, anthropogenic presence (roads) and rivers identified as barriers to gene flow while metapopulation connectivity was enhanced by coastal habitats (beaches, dunes, sand and salt marshes) and coastal grassland. Substantial population declines over the past 15 years necessitate increased conservation management efforts by the Government: principally artificial pond creation and assisted migration through translocation. Our results are invaluable in informing planned improvements in connectivity by suggesting sites where corridors of new ponds may be beneficial and the clusters of sites that can be used as source populations for translocation while maintaining genetic provenance. Moreover, we identified habitats with restricted dispersal and suggest site management measures to mitigate their effects, for example, prevention of scrub encroachment. Thus, this study is valuable in understanding how landscape affects dispersal, gene flow and habitat connectivity of a declining pond-breeding amphibian. The same approach could be used in other regions to direct and inform conservation managers.
in data collection and we thank the landowners and farmers who provided access to their land.
Author contributions MR: conducted field work, laboratory and statistical analysis. KD: provided guidelines on laboratory procedures including DNA extraction methods, selection of primers, PCR optimization and interpretation of results. JF and PN: prepared the climate data and contributed to the data analysis. JT, AA and ME: assisted with planning the fieldwork, verified the analytical methods and provided critical feedback on data interpretation. FM: contributed to the design of the study, substantially contributed to revising the manuscript. SH: supervised the molecular work and statistical analysis of genetic data, contributed to the study conception and designed. NR: was the Principal Investigator contributing throughout from conception, methodology, and data analysis to manuscript editing. The first draft of the manuscript was written by MR. All authors discussed the results, contributed to the draft, read the final manuscript and gave approval for publication.

Data availability
The dataset generated during this research will be deposited on DRYAD. The datasets generated during the current study will also be available from the corresponding author on reasonable request.

Code availability
No custom codes were used; software and used codes are cited in the text.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.