Population genetic structure of the Natterjack toad (Epidalea calamita) in Ireland: implications for conservation management

Molecular methods can play a crucial role in species management and conservation. Despite the usefulness of genetic approaches, they are often not explicitly included as part of species recovery plans and conservation practises. The Natterjack toad (Epidalea calamita) is regionally Red-Listed as Endangered in Ireland. The species is declining and is now present at just seven sites within a highly restricted range. This study used 13 highly polymorphic microsatellite markers to analyse the population genetic diversity and structure. Genetic diversity was high with expected heterozygosity between 0.55 and 0.61 and allelic richness between 4.77 and 5.92. Effective population sizes were small (Ne < 100 individuals), but not abnormal for pond breeding amphibians. However, there was no evidence of historical or contemporary genetic bottlenecks or high levels of inbreeding. We identified a positive relationship between Ne and breeding pond surface area, suggesting that environmental factors are a key determinant of population size. Significant genetic structuring was detected throughout the species’ range, and we identified four genetic entities that should be considered in the species’ conservation strategies. Management should focus on preventing further population declines and future loss of genetic diversity overall and within genetic entities while maintaining adequate local effective population size through site-specific protection, human-mediated translocations and head-start programs. The apparent high levels of genetic variation give hope for the conservation of Ireland’s rarest amphibian if appropriately protected and managed.


Introduction
Given the current global 'amphibian crisis' (IUCN 2020), comprehensive understanding of declining and threatened species' genetic structure is essential for effective conservation strategies and population management (Emel and Storfer 2012). Small population size promotes loss of rare alleles, increased homozygosity and accumulation of detrimental recessive alleles (Westemeier et al. 1998;Madsen et al. 1999), resulting in local declines and elevated extinction risk (Frankel and Soule 1981;Shaffer 1990;Hedrick 2001;Emel and Storfer 2012). Amphibians are particularly at risk of reduced genetic diversity due to their sensitivity to environmental change, low dispersal capability, and strong natal philopatry Beebee 2005;Kraaijeveld-Smit et al. 2005).
Low genetic variability in amphibians has been linked to an increased presence of physical abnormalities (Hitchings and Beebee 1988) and negative impacts on oxygen consumption (Mitton et al. 1986), immune response (O'Brien and Evermann 1988;Altizer et al. 2003;Cabido et al. 2010Cabido et al. , 2011, ability to compete for resources (Rowe and Beebee 2005), clutch size (McAlpine 1993), hatching success (Blaustein et al. 1994) and growth rates (Rowe and Beebee 2004). For example, loss of genetic diversity in the Italian agile frog (Rana latastei) has been related to higher susceptibility to Ranavirus (Pearman and Garner 2005). Such findings are a concern given the major role disease plays in global amphibian declines (Daszak et al. 2003;Skerratt et al. 2007; Wake and Vredenburg, 2008). Thus, omitting genetics from design of species management plans may lead to inappropriate allocation of resources, ineffective conservation strategies, and further population declines (Frankham 2003;Cushman 2006;Noel et al. 2007).
Species genetic management is often overlooked in conservation programmes due to lack of resources or expertise (Taylor et al. 2017). Understanding of the genetic structure of endangered populations can inform a targeted management approach. Genetic analysis permits assignment of individuals to genetic clusters independent of the sampling location and identifies discrete units for management (Moritz 1994;Palsbøll et al. 2007;Schwartz et al. 2007). Below the species-level, conservation priorities focus on intraspecific diversity for evolutionary plasticity and adaptability in response to environmental change (Shaffer et al. 2015). For instance, genetic analysis of the Australian agamid lizard (Diporiphora nobbi) identified several Evolutionary Significant Units (ESUs), emphasising the importance of local populations of a widespread species in harbouring intraspecific genetic diversity (Driscoll and Hardy 2005). Conversely, where genetic diversity is critically low, conservation managers may need to actively consider genetic rescue by deliberately maximising assisted gene flow between distinct units (Beauclerc et al. 2010;Whiteley et al. 2015). An isolated population of adders (Vipera berus) in Sweden dramatically increased in number after the introduction of new genes to avoid severe inbreeding averting local extirpation (Madsen et al. 1999(Madsen et al. , 2004. Thus, the approach taken to genetic management is species-, population-and context-dependent. Effective population size (N e ) plays a crucial role in predicting a population's extinction risk and is often more valuable than estimating absolute census size (Beebee 2005). It is defined as the size of an idealized population that meets Hardy-Weinberg assumptions showing the same rate of loss of genetic diversity as the target population (Wright 1931). Effective population size can be derived from genetic data with no additional life-history information (Schwartz et al. 1998). Furthermore, molecular markers allow historical patterns of population decline or expansion to be inferred, informing contemporary management. There are many factors that can influence N e , including census size, sex ratio and past declines, breeding strategy, habitat carrying capacity and connectivity between populations (Waples and Gaggiotti 2006;Mills 2007). Wang et al. (2011) found that suitability of breeding habitats can contribute to variation in N e among populations of the California tiger salamander (Ambystoma californiense). Thus, understanding the relationship between demographic parameters, key habitat features, and effective population size is essential for a comprehensive conservation approach.
Our study assesses the impact of small population size and recent declines on the genetic structure of the Natterjack toad (Epidalea calamita), to inform species conservation management. In Ireland, the species is at its northwesternmost range edge margin (Gasc et al. 1997) and is highly range restricted, occurring in seven isolated breeding sites on the coast of Co. Kerry ( Fig. 1; Reyne et al. 2021a). One of those breeding sites (Caherdaniel) was established in the 1990s as a result of introduction of individuals from the Magharees to increase the species geographic range. Regardless of conservation efforts, the species is in decline and is regionally Red-Listed as Endangered in Ireland (King et al. 2011). The Natterjack appears native to Ireland and, despite reporting an apparent bottleneck, May and Beebee (2008) suggested its genetic diversity was "reassuringly high" despite a lack of gene flow from other populations. However, low levels of polymorphism in the selected microsatellite markers may have contributed to their low reported allelic richness (1.86-2.89). Since then, the Natterjack toad population in Ireland has continued to decline, with > 90% loss of spawning at some breeding sites (Reyne et al. 2019(Reyne et al. , 2021a. In an effort to improve habitat availability, over 100 new breeding ponds were created as part of an agrienvironment scheme since 2008 by the National Parks and Wildlife Service (NPWS). NPWS now manages a 'head-start and translocation' programme through which eggs strings and tadpoles are collected annually, raised in captivity and metamorph toadlets released back into the wild, supplementing existing populations as part of assisted migration and translocation to newly created ponds. Since 2016, 13 translocations have been performed to aid colonisation of artificial ponds. As yet, there are no data available to assess the head-start and translocation program success, as breeding will occur 4 to 5 years post release when toads have reached sexual maturity and return to the ponds for breeding (Beebee 1979). A primary goal of ex situ conservation is to preserve maximum intraspecific genetic variability to ensure long-term survival (Pelletier et al. 2009). However, inbreeding is a major problem in many ex situ breeding programs (Witzenberger and Hochkirch 2011;Santymire et al. 2019;Phillips et al. 2020) with success often determined by the appropriate selection of founders from source populations (Tzika et al. 2009). Selecting individuals and populations for ex situ breeding and translocations should, therefore, be informed by empirical genetic data ensuring provenance and/or appropriate management of genetic lineages (Witzenberger and Hochkirch 2011;IUCN/SSC 2013). Hence, the specific objectives of our study were to: (i) provide genetic characterisation of each breeding site, (ii) reconstruct parentage of offspring samples, (iii) estimate the effective and census population sizes, (iv) evaluate the impact of pond characteristics on effective population size, (v) detect any genetic bottleneck effect(s) due to historical or recent declines, (vi) quantify genetic differentiation among breeding sites, (vii) designated discrete conservation genetic entities, if appropriate, and (viii) provide management recommendations.

Field surveys and sample collection
We conducted fieldwork in 2017 during the Natterjack toad's breeding season (April-July). Breeding ponds were visited every 7 to 10 days throughout the duration of the breeding season. We recorded the total number of egg strings in each pond and measured the size (surface area) of each breeding site. For full details of sampling methods see Reyne et al. (2019). We collected minimum 30 DNA samples per breeding site from different, well-developed Natterjack toad egg strings (embryos with clearly defined head and tail) within the same pond or whole tadpoles from different breeding ponds to avoid analysing full siblings (see Fig. 1 for the location of the breeding sites and number of breeding ponds/areas per site). We did not sample any egg strings or tadpoles that were part of the headstart and translocation program i.e., translocated between sites. Samples were preserved in 100% ethanol at ambient temperature until extraction.

Genotyping and data validation
Genomic DNA was extracted following a high salt extraction protocol (Supporting Protocol S1). Thirteen highly polymorphic fluorescently labelled microsatellite markers were amplified (Table 1). PCR reactions were performed in two multiplexes, and forward primers were labelled with 6-FAM™, NED™ and VIC™ fluorescent dye (Applied Biosystems, Integrated DNA Technologies). Multiplex PCR reactions had a total volume of 10 µl, and contained 1 µl of genomic DNA, 5 µl Type-it Multiplex PCR Master Mix (Qiagen) and primer mix of labelled forward and reverse with 0.1-0.3 µM final concentrations. PCR conditions were: an initial denaturation of 95 °C for 15 min; followed by 35 cycles of 94 °C for 30 s, annealing temperature (Ta) of 55 °C or 58 °C for 90 s for multiplex one and two respectively, and 72 °C for 60 s, with a final extension at 60 °C for 30 min. Samples were randomized across breeding sites during the genotyping analysis to avoid bias and negative controls were used throughout. PCR products were diluted ten times with ddH 2 O, and 1 µl diluted PCR product, was added to 8.95 µl Hi-Di™ Formamide (Thermo Scientific) and 0.05 µl GeneS-can™ 600LIZ ladder (Applied Biosystems). Fragment analysis was performed on ABI 3730xl Genetic Analyzer (Applied Biosystems).
Alleles were scored with GeneMarker® V1.8. We rounded all genotype calls to an even or odd number based on the known motif of respective locus. Differences between assigned and actual allele size were between 0.1 and 1 bp. To calculate potential genotyping errors (Bonin et al. 2004;Pompanon et al. 2005), 28 samples (approx. 10%) were selected randomly using random number function in Microsoft Excel (2016), genotyped three times, and genotyping error rates (allelic dropout and false allele) were calculated using PEDANT v.1.0 software (Johnson and Haydon 2007). To calculate error rates from three repeated genotypes, we compared each repeat and averaged the error estimates, as recommended by Johnson (2007).

Genetic diversity
To quantify genetic diversity, we estimated allelic richness (A R ), observed heterozygosity (H O ), expected heterozygosity (H E ) and inbreeding coefficient (F IT ) for each locus and breeding site using FSTAT version 2.9.4 (Goudet 1995) and Genepop 4.2 (Raymond and Rousset 1995). We also calculated rarefied private allelic richness (A P ) in HP-rare (Kalinowski 2004(Kalinowski , 2005. Compliance with Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) was tested among loci. Statistical significance was adjusted for multiple comparisons using strict Bonferroni correction.

Polygyny
Polygyny is a common anuran breeding strategy. While multiple paternity of egg clutches as a result of polyandrous mating or clutch piracy have been observed in several amphibian species (e.g., VietIs et al. 2004;Knopp and Merilä 2009;Byrne and Roberts 2012), it is unlikely for the Natterjack toad due to species segregated-pair breeding behaviour (May et al. 2011). Polygyny was estimated following Ficetola et al. (2010). We used the computer software COLONY (Jones and Wang 2010) to estimate half-sibling groups among the sampled eggs and tadpoles. We assumed that all eggs and tadpoles had different mothers based on the sample collection method (i.e., samples were collected from different egg strings or different breeding ponds), but that one male could fertilise more than one egg string. The reconstructed sibships were used to assess the number of egg strings fertilized by each male. We identified all half-siblings and calculated the degree of polygyny at a breeding site level as the average number of egg strings fertilized by reproductive males. Eggs and tadpoles without half-siblings were treated as the offspring of males that fertilized only one egg string. Variance in reproductive success at a breeding site level was calculated as the variance in the number of egg strings fertilized by breeding males. We estimated Spearman's rank correlations to investigate relationships among polygyny, variance in reproductive success and number of egg strings. All statistical analyses in the study were performed in R 4.0.0. using the stats package (R Core Team 2019) unless stated overwise.

Population measures
We used egg string counts to monitor population size, a widely used method by previous Natterjack toad surveys in Great Britain (Smith and Skelcher 2019) and Ireland (Bécart et al. 2007;Sweeney et al. 2013;Reyne et al. 2019). Females typically deposit one egg string per year (Buckley and Beebee 2004); thus the number of egg strings can be used as a proxy for the female breeding population size. We  Rowe et al. (2000) plotted the cumulative number of egg strings at each pond visit for each population and the asymptote was calculated using a generalized additive model (GAM) which allowed 95% confidence intervals to be associated with counts. The lower confidence limit (LCL) accounted for potential double counting and the upper confidence limit (UCL) accounted for egg strings potentially being missed. We estimated effective population size (N e ) using LDNe 1.0 software (Waples and Do 2008). Any bias caused by a small sample size is corrected by implementing the biascorrection method of Waples (2006). We ran the model with three different critical values (Pcrit) of allele frequency. Pcrit was set at 0.05, 0.02 and 0.01 which means that all alleles with frequency lower than the critical value were excluded from analysis (Waples 2006). We used linear regression to evaluate the role of number of breeding females in explaining N e .

Breeding area
Total available area to breed at each breeding site was calculated as a sum of pond surface areas where breeding activity was recorded. We performed linear regressions between available breeding surface area and each of the following response variables: effective population size (N e ), mean observed heterozygosity (H O ), mean allelic richness (A R ) and polygyny to examine the relationship between breeding pond characteristics and genetic diversity.

Bottlenecks
BOTTLENECK v.1.2.0.2 was used to identify that recently had undergone a reduction in effective population size through testing for significant deviations from the mutationdrift equilibrium. The program was run under two mutational models for microsatellite data: two-phased (TPM) and stepwise mutation (SMM). The TPM model was run with 95% single-step mutations, 5% multiple-step mutations and variance among multiple steps of 12, as recommended by Piry et al. (1999). A Wilcoxon single-rank test was used to test for heterozygosity excess. Significant results of heterozygosity excess indicate evidence of a recent reduction in the effective population size (Piry et al. 1999). We also performed a two-way ANOVA between mean observed heterozygosity and mean allelic richness across markers and breeding sites to compare genetic diversity between stable and declining populations.

Population genetic structure
We performed discriminant analysis of principal components (DAPC), a multivariate analysis which is free from HWE assumptions and shows greater resolution of population structure when levels of genetic differentiation are low (Jombart et al. 2010). We used snapclust clustering algorithm to identify optimal number of genetic clusters (K) and assign individuals to panmictic populations (Beugin et al. 2018). The optimal number of clusters to retain in the analysis was identified using the Akaike Information Criterion (Akaike 1998) and the Kullback Information Criterion (Cavanaugh 1999) where lower values indicate better fit. Analyses were carried out in the adegent package in R following Jombat (2008) and Beugin et al. (2018).
F ST (Weir and Cockerham 1984) was used to investigate patterns of differentiation among populations and as an indirect measurement of historical gene flow. We tested for significance across these comparisons using 10,000 permutations. Statistical significance was adjusted for multiple comparisons using strict Bonferroni correction. Analysis was done using R package hierfstat (Goudet 2005) and FSTAT v2.9.4 (Goudet 1995). Genetic structure was also assessed using analysis of molecular variance (AMOVA, Excoffier et al. 1992) which calculates the proportion of the total genetic variance originating within and between breeding sites. Permutation tests were performed at three hierarchical levels: among breeding sites, among individuals within breeding sites and within individuals. Analyses were conducted in Arlequin v3.5.2.2. (Excoffier and Lischer 2010).

Results
In total, 316 Natterjack toad samples were collected from all breeding sites in Ireland and successfully genotyped. Sample sizes varied between 32 and 58 individuals ( Table 2). All 13 microsatellite loci amplified successfully. Genotyping error rate varied among loci ranging from 0 to 0.10 for allelic dropout and from 0 to 0.03 for false alleles (Table S1). All microsatellite markers were polymorphic with allelic  Table 2). The mean inbreeding coefficient across all populations was low (ranged 0.04-0.16; Table 2). The breeding sites and overall, the Natterjack toad population in Ireland were not in HWE (p = 0.01). The number of egg strings fertilized by individual males was between 1 and 8. The degree of polygyny at a population level varied between 2.13 and 3.73 ( Table 2). Variance of male breeding success differed strongly among populations and ranged from 1.98 to 13.22 (Table 2). There was no significant relationship between polygyny, variance of breeding success and number of egg strings i.e., breeding females (p > 0.05).
In total, we recorded 1457 egg strings. There was a significant sigmoidal accumulation of egg strings for all populations resulting in narrow 95% confidence intervals ( Fig.  S1; Table 3). The narrow 95% confidence intervals for most populations are indicative of reasonably precise estimates. Based on the number of egg strings, the most productive population was at the Magharees sand dune and the least productive at Inch (Table 3). Estimates of N e were low, ranging from 19 at Caherdaniel to 519 at Yganavan but typically < 100 individuals at most sites (Table 3). There was no significant relationship between the number of breeding females and N e regardless of the critical value of allele frequencies used to calculate N e (p > 0.05).
BOTTLENECK analysis suggested significant deviations from mutational drift equilibrium for some breeding sites (Table 4). However, there is no evidence for significant heterozygosity excess at any breeding site. In addition, there was no evidence of mode-shift in allele frequencies as all distributions were L-shaped suggesting no recent reduction in the effective population size. Additionally, a two-way ANOVA showed no significant difference in the allelic richness (F df=6,84 = 0.594, p = 0.735) and observed heterozygosity (F df=6,84 = 0.490, p = 0.814) among breeding sites.
Cluster analysis suggested an optimal number of K = 5 using AIC and K = 4 using KIC (Fig. S2). KIC approach indicated a better fit (sharp decrease in the values followed by a sharp increase), hence the number of clusters 1 3 retained in the analysis was four. DAPC indicated that the Magharees differed from the rest of the breeding sites with no overlap (Fig. 2A). The Caherdaniel was also distinct from all other breeding sites, although it had similar DAPC axis 1 scores (which carried most of the observed variation) to the Magharees breeding site from which it was derived. The Castlemaine Harbour sites (Glenbeigh, Yganavan, Inch, Dooks and Roscullen) clustered together but re-analysis after the exclusion of the Magharees and Caherdaniel suggested Roscullen was distinct from the rest whilst Glenbeigh had minimal overlap (Fig. 2B). Snapclust analysis revealed similar patterns where Magharees and Caherdaniel breeding sites were different from the rest of the breeding aggregations; however it failed to clearly differentiate among the Castlemaine Harbour breeding sites (Fig. S3). Based on our analysis, we identified four genetic entities: (1) Magharees, (2) Caherdaniel, (3) Roscullen and (4) the remainder of the Castlemaine Harbour sites i.e., Dooks, Glenbeigh, Yganavan and Inch. AMOVA analysis suggested most molecular variation (96.43%) was attributed to within individuals while variation between breeding sites accounted for 4.58% (Table 5). AMOVA detected significant (p < 0.001) differentiation among sites. The pairwise F ST values suggested higher levels of gene flow between Glenbeigh, Yganavan, Dooks and Inch with F ST ranging between 0.016 and 0.020 (Table 6). Higher F ST values (> 0.05) were observed between Magharees and the rest of the populations apart from Caherdaniel.

Discussion
The regionally Red-Listed Irish Natterjack toad population, whilst having undergone substantial declines in census population size (Reyne et al. 2021a), exhibited no evidence of genetic bottlenecks or inbreeding with relatively high genetic diversity (allelic richness and heterozygosity), despite low effective population sizes. Analysis of population structure suggested four distinct genetic entities that should be considered in species conservation programmes.
All seven Natterjack toad breeding sites surveyed in Ireland were polymorphic at the 13 microsatellite loci. Expected heterozygosity (0.55-0.61) and observed heterozygosity (0.49-0.59) were higher than previous estimates for the same populations in Ireland and higher than most estimates for Natterjack toad populations throughout Europe (Table S4). However, estimates of heterozygosity and allelic richness are not directly comparable due to the use of a different set of microsatellite markers. It is also important to emphasis, that neutral and selected genes can provide different information regarding species genetic diversity and adaptative capacity. In this study, we used microsatellite markers consequently focusing only on the putatively neutral genetic diversity. However, it is increasingly recognised that adaptive genetic diversity is a more important indicator of species ability to adapt to changing environments and population persistence in the wild. Hence, shifts towards genome-wide single nucleotide polymorphisms (SNPs) or a combined approach can Pond breeding amphibians are predisposed to lower levels of genetic variation compared to other taxa resulting from high amplitude fluctuations in population size (Alford and Richards 1999;Newman and Squire 2001). For example, the Natterjack toad typically breeds in ephemeral ponds where breeding success can depend on stochastic climatic One of the key factors determining population genetic diversity is effective population size, a metric that provides additional insights for management that complements census population size. Low N e values can potentially indicate loss of genetic variability within a population (Frankham 2003;Ficetola et al. 2007;Palstra and Ruzzante 2008). Our estimates of the effective population size were low for all seven remaining populations (mainly < 100 individuals) despite reasonably observed high genetic diversity. The N e for pond breeding amphibians is typically in tens rather than hundreds or thousands of individuals regardless of large census sizes . Similar estimates of the effective population size were reported for ten Natterjack toad populations in Europe (Faucher et al. 2017) and other amphibian species including the marsh frog (Rana redibunda, Zeisset and Beebee 2003), Italian agile frog (R. latastei, Ficetola et al. 2010) and common frog (Rana temporaria, Brede and Beebee 2006). This generally reflects the most common anuran breeding strategy of scramble competition where large males dominate breeding ponds and available females resulting in only a few males contributing to the gene pool of the next generation (Ficetola et al. 2010). Higher levels of polygyny and variance of male breeding success were recorded at Caherdaniel, Dooks and Glenbeigh where few, small breeding ponds are available. Variation in breeding habitat has previously been associated with differences in effective population size (Wang et al. 2011), and this is true for our study. Natterjack toad N e was significantly associated with the breeding surface area available. The highest N e was calculated for the Yganavan population where toads breed along the shore of a large lake, suggesting less competition among males (i.e., low polygyny). Populations with small available breeding habitat had smaller N e , but it did not result in lower genetic diversity. Similar results were found for the California tiger salamander (A. californiense, Wang et al. 2011) highlighting the need for a more complete understanding of the parameters influencing N e .
Despite small N e and substantial declines in egg string numbers (breeding females) since 2004 for populations at Roscullen, Yganavan, Dooks and Glenbeigh (Reyne et al. 2019;2021a), our analysis did not provide evidence of recent genetic bottlenecks. Dooks and Yganavan populations have declined by over 70% (Reyne et al. 2021a); however, migration between these sites may decrease the number of rare alleles, consequently masking any H E excess (Cornuet and Luikart 1996) or a lag between when population size reduction took place and an observable increase in the inbreeding coefficient suggesting that changes in the genetic diversity of the offspring are yet to be observed. The BOTTLENECK program has been reported to fail to detect known demographic bottlenecks (Whitehouse and Harley 2001) and   (Spencer et al. 2000;Whitehouse and Harley 2001). The Magharees had the largest Natterjack toad population in Ireland accounting for 90% of the recorded egg strings deposited annually and this number has remained largely stable over time (Reyne et al. 2021a). However, allelic richness for the Magharees was among the lowest (4.79) in Ireland, though there was no significant difference in allelic richness among breeding sites. This would support the assertion of no decline in the effective population size consistent with a lack of evidence of a genetic bottleneck. Similarly, the Australian bell frog (Litoria aurea) has undergone an 80% population decline over the past 40 years with no apparent genetic bottleneck with levels of allelic richness not significantly different for 19 out of 21 populations when compared to a large demographically stable population (Burns et al. 2004). The Natterjack toad in Ireland exhibits genetic structuring throughout its range with all sampled locations being significantly different. DAPC results imply that the Natterjack toad population in Ireland can be divided into four discrete management units: (1) Magharees, (2) Caherdaniel, (3) Roscullen, (4) Dooks, Glenbeigh, Yganavan and Inch. Amphibians exhibit strong site fidelity and limited dispersal and migration, resulting in low levels of gene flow among populations (Reading et al. 1991;Kusano et al. 1999;Pittman et al. 2008). Several studies suggested that highly structured populations are often typical of amphibians with distinction at scales less than 5 km (e.g., Shaffer and Breden 1989;Routman 1993;Driscoll 1998) questioning the notion of metapopulation dynamics (Marsh and Trenham 2001). For amphibian species with highly structured populations site-specific protection and human-mediated translocations/ reintroductions may be critical management tools to preserve intraspecific genetic diversity (Shaffer et al. 2000).
Conservation recommendations should focus on maintaining high genetic diversity as well as protection of the genetic integrity of identified genetic entities. This can be achieved by maintaining adequate effective population size (Storfer et al. 2007;Wang 2009) especially in small and fragmented populations (Wang et al. 2011). The NPWS is not currently creating new ponds but should future measures to incentivise farmers to create new potential breeding ponds, and these should be clustered within and around each identified management unit. Furthermore, our results suggest that having large breeding ponds or high numbers of small ponds in close proximity can be particularly valuable for ensuring a large N e associated with large breeding surface area and higher number of successfully mating male toads. The benefits of multiple pools in supporting Natterjack toad's metapopulation structure have been clearly demonstrated in Britain after decades of active species conservation (Denton et al. 1997;Beebee 2014). The ongoing NPWS 'head-start and translocation' programme should be cognisant of our identified genetic entities when releasing toadlets back into the wild respecting their genetic provenance, not translocating individuals between genetic entities and selecting the source populations from geographically proximate sites within the same entity. The only exception is Caherdaniel, as the population is likely distinct due to founder effects, subsequent gene drift and lack of gene flow after its establishment using translocated individuals from the Magharees. Thus, should further translocations be required in the future to maintain the population these can be drawn from the original source population. The lowest observed heterozygosity and highest inbreeding coefficient were recorded for the Glenbeigh population. Inbred Natterjack toad's tadpoles have slower growth rates and lower survival rates (Rowe and Beebee 2005). Thus, we have conservation concerns about the long-term survival of the Glenbeigh population. A genetic restoration program of a small and isolated Natterjack toad colony at Saltfleetby, Britain has led to increase in metamorph production and breeding success (Beebee 2014). Hence, consideration should be given to population supplementation at Glenbeigh with translocated individuals ideally from Yganavan as it is genetically the closest population, but alternatively from any population within the management unit. It is important to note that while genetically most relevant population is preferable as a source population, it becomes unfavourable if removal of those individuals might put the source population at risk. Hence, other factors like population size and trend, potential threats and pressures should be considered when selecting source populations.
There may be an aspiration that the Natterjack toad's range in Ireland could be enlarged beyond its current highly restricted range, for example, by introducing animals to sites in the southwest occupied historically or introducing animals to suitable habitat in the west more generally which has never been occupied (for example, assisted migration of populations northward to track climate change). In these cases, donor populations should either be geographically proximate or, if a new location is distinct from existing genetic entities and isolated by barriers to dispersal, the Magharees population can be considered as donors due to the high genetic diversity and large population size. An experimental approach could be taken by mixing individuals from each of the genetic entity to create artificially high genetic diversity buffering any new population against local extirpation due to small initial effective population size (Beauclerc et al. 2010;IUCN/SSC 2013). Monitoring such populations would be warranted.
Our findings have important conservation implications for the management of the Natterjack toad, which is regionally Red-Listed as Endangered in Ireland.

3
Populations appear to lack deficiency in allelic richness or heterozygosity despite low effective population sizes and there is no evidence of genetic bottlenecks despite declines in census size. We identified discrete genetic entities, which we urge species conservation programmes to consider when undertaking population supplementation, translocation or assisted migration. We also highlight the importance of long-term genetic monitoring especially of declining populations and at reintroduction sites. Apparent high levels of genetic variation gives hope for the conservation of Ireland's rarest amphibian.