Genetic and demographic vulnerability of adder populations: Results of a genetic study in mainland Britain

Genetic factors are often overlooked in conservation planning, despite their importance in small isolated populations. We used mitochondrial and microsatellite markers to investigate population genetics of the adder (Vipera berus) in southern Britain, where numbers are declining. We found no evidence for loss of heterozygosity in any of the populations studied. Genetic diversity was comparable across sites, in line with published levels for mainland Europe. However, further analysis revealed a striking level of relatedness. Genetic networks constructed from inferred first degree relationships suggested a high proportion of individuals to be related at a level equivalent to that of half-siblings, with rare inferred full-sib dyads. These patterns of relatedness can be attributed to the high philopatry and low vagility of adders, which creates high local relatedness, in combination with the polyandrous breeding system in the adder, which may offset the risk of inbreeding in closed populations. We suggest that reliance on standard genetic indicators of inbreeding and diversity may underestimate demographic and genetic factors that make adder populations vulnerable to extirpation. We stress the importance of an integrated genetic and demographic approach in the conservation of adders, and other taxa of similar ecology.


Introduction
Population genetics remain overlooked in conservation planning [1], although genetic factors may lead to population extinction even after other threats have been addressed [2]. Loss of genetic diversity and inbreeding depression represent the primary genetic threats [3], with the potential to contribute to an extinction vortex [4].
Tingley et al [5] have stressed the importance of addressing the optimal genetic management of small isolated reptile populations. The adder Vipera berus (Linnaeus, 1758) is a terrestrial snake with an extremely wide geographic range [6], which accounts for its IUCN  status as "Least Concern", although with a decreasing population trend [7,8]. Adders, like many other temperate snakes, are viviparous with low fecundity, low vagility and high philopatry [9,10], life-history traits that render them vulnerable to local extinction [11,12]. The negative outlook for adder populations is exacerbated by snakes being among the least popular terrestrial vertebrates, more likely to be targets of intentional killing than of conservation management [13]. The potential risk of inbreeding depression in adders is highlighted by an isolated adder population in Sweden, in which a decline in numbers was associated with stillbirths and deformities, and a reduction in genetic diversity, all of which responded to the introduction of adult males from a large outbred population [14][15][16]. Ú jvári et al [17] have similarly reported low juvenile survival and birth deformities with reduced genetic diversity and increased homozygosity in fragmented populations of the congeneric Hungarian meadow viper, V. ursinii rakosiensis. Small population size is an important factor in loss of genetic diversity, exacerbated by bottleneck events. This has led to the concept of a minimum viable population size [18,19], based on the inverse relationship between the effective population size (N e ) [20] and the rate of erosion of genetic variation by drift, which is supported by studies of wild populations [21,22]. N e estimates tend to be low in relation to census population size in natural populations [23,24], influenced by demographic fluctuation and life-history traits [25,26]. Both small population size and genetic erosion render populations more susceptible to stochastic environmental and demographic adverse events, such as climate change or disease [19,27,28]. Small populations isolated by habitat fragmentation are also at increased risk of inbreeding. However, several important questions regarding genetic variation and inbreeding depression in natural populations remain largely unanswered. In particular, it is unclear to what extent mating between close relatives and loss of genetic diversity contribute to population decline and extinction in the wild, and thus to how results of genetic studies should influence their conservation management [29].
We report the results of the UK Adder Genetic Project (UKAGP), a study into the genetic status of lowland adder populations in southern mainland Britain, where national distribution surveys have indicated a decline in comparison with historic records [30][31][32][33]. A national questionnaire-based investigation survey of adder populations showed that declines were more likely to be reported in small sites with fewer than ten adders, whether based on systematic surveys or anecdotal evidence [31]. The subsequent Make the Adder Count (MTAC) initiative, based on peak springtime adder counts over sequential years, further underscored the increased risk of decline in small populations [34], flagging threats of habitat loss, public disturbance and predation, especially by cats and birds. To these threats should be added the potential risk of disease caused by the release of captive non-native snakes onto adder habitat, especially in view of the recent finding of the causative agent for snake fungal disease (Ophidiomyces ophiodiicola) in UK adders in the wild [35]. This consistently emerging picture of habitat fragmentation and local decline forms the background for our study, in which we have used a combination of mitochondrial DNA (mtDNA) and microsatellite markers to investigate the potential role of genetic factors in their decline of adders in mainland Britain.
The aim of this study was to document population genetic structure and differentiation, and to estimate indicators of inbreeding and genetic diversity in lowland adder populations in southern mainland Britain. To assess the ecological and conservation significance of our study, we interpret our results in comparison with published studies of adder populations in mainland Europe, and with reference to the size, and thus the likely risk status, of the study populations.

UKAGP study sites and samples
Ethics and animal welfare: the project was reviewed prospectively and approved by the ZSL Ethics Committee. Sampling was undertaken in concordance with ZSL ethical guidelines. Cloacal swabs were collected without anaesthetic from adult snakes, and buccal swabs from juveniles [36], by ecologists experienced in snake handling. All snakes were released at the site of capture. Handling of adult females was avoided after mid-May, to minimise disturbance to gravid snakes. No study animal was subjected to euthanasia.
Samples were collected from 220 adders at 16 sites in southern mainland Britain between March and May 2011 (Table 1). No permits were required, as there were no restrictions on site access, and the adder in the UK has no specific protection status other than against deliberate injury or killing, or collection for trade (Wildlife and Countryside Act 1981, as amended 1991). For each site, adders were caught over a one to two-day period. Dates and optimal weather conditions for sampling were determined according to local ecological expertise.
DNA was extracted using a DNeasy blood and tissue kit (Qiagen), following the manufacturer's protocol for swabs.

Mitochondrial DNA (mtDNA) sequencing
We designed primers from a 918-bp mtDNA control region (CR) sequence and a 1043-bp mtDNA cytochrome b sequence (Cytb) [37], based on high levels of variability across 40 European adder haplotypes, and on consistent flanking region stability. The primer sequences selected for this study are shown in Table 2A. PCR for both loci was performed in HotStarTaq-Plus (Qiagen), with an annealing temperature of 55˚C. The same primers were used for sequencing reactions. PCR products were cleaned using a QIAquick PCR purification kit (Qiagen, UK), and sequenced using ABI BigDye1 chemistry and 3130XL sequencer, following manufacturers' protocols. We generated alignments of concatenated Cytb/CR consensus sequences using MEGA 6 [38]. Individual haplotypes were identified using the haplotype function of pegas package v0.10 [39], implemented in R v3.4.0 [40]. A haplotype network [41] was constructed using the median-joining method [42] in NETWORK v4.6.1 (www.fluxusengineering.com). Phylogenetic analysis was carried out as described in Supporting Information (S1 Table,

Microsatellite genotyping
In preliminary studies we tested published microsatellite primers that had been developed for adders (V. berus) [43,44], selecting five polymorphic loci that we found to amplify consistently, with a minimum of stutter bands. To increase the number of loci, we also evaluated 15 congeneric microsatellite markers which had been developed for meadow vipers (V. ursinii) [45], selecting three that successfully amplified and demonstrated polymorphism in adder samples (results not shown). PCR was performed in 10 μl volumes with 20-100 ng DNA, 5 μl mastermix (HotStarTaq Plus or Multiplex; Qiagen), 5 μmol/L unlabelled reverse primer and 5 μmol/L fluorophore-labelled forward primer (Applied Biosystems). Amplification was performed in simplex with initial denaturation 95˚C 5 min, 60 cycles of 94˚C 60 sec, 57-59˚C 60 sec, 72˚C 60 sec, and final extension 72˚C 7 min, optimized in preliminary studies for each primer pair. Primer sequences and locus-specific PCR conditions are summarized in Table 3. Amplified products were resolved by capillary electrophoresis on a 3130xl Genetic Analyser with a LIZ-500 size standard (Applied Biosystems). Alleles were scored and binned manually, using PeakScanner 1.0 software (Applied Biosystems).

Microsatellite data analysis
Quality control. Replicates and template negative controls were included in all plates to confirm reproducibility of results. Results were analysed for genotyping errors and null alleles in Micro-Checker v 2.2.0.3 [46]. We used FSTAT v2.9.3.2 [47] and pegas [39], implemented in R, to test for Hardy-Weinberg equilibrium (HWE), and to exclude linkage disequilibrium.
Measures of population structure and differentiation. In this study, we use the term population to refer to all individuals sampled at a single study site in the sampling time period. Pairwise F ST values between populations were estimated in FSTAT.
To test for isolation by distance [48] we applied the mantel.rtest function of ade4 v1.7-13 [49], implemented in R, with 999 repetitions, using pairwise F ST to estimate genetic distance. Geographic distance was estimated at https://andrew.hedges.name/experiments/haversine (accessed 23 April 2108) using the Haversine great circle method, a measure of the shortest distance between two points on a sphere [50].
To investigate genetic differentiation, we used STRUCTURE v2.3 [51,52], using correlated allele frequencies and admixture models, with or without the locprior option [53]. The initial alpha was set at 1/n, where n is the number of sample locations, to allow for variation in sample sizes between populations [54]. We used burn-in of 10 5 , followed by 10 6 iterations for 10 independent replicate runs for values of K from 1 to the number of populations being studied. Results were uploaded to StructureHarvester [55] to derive mean log likelihood and delta-K as a function of K, detecting hierarchical levels of structure [56]. Results across replicate runs were permuted using the greedy function of CLUMPP [57] to derive proportional assignments to each cluster for supported values of K. We also studied population structure using discriminant analysis of principal components (DAPC) [58] in adegenet version 2.0.1 [59] implemented in R. DAPC is a multivariate method to identify clusters of genetically related individuals, which is not based on a predefined model, and makes no assumptions of HWE. The find.clusters function was applied to determine the optimal number of clusters (k) in each population, independent of the number of sampling sites. The dapc function was then applied, using the α-score function to determine the optimum number of principal components to retain in each analysis. Probabilities of assignment of individuals to each of the different clusters were visualised using the compoplot function of adegenet [59].
Microsatellite summary statistics: Baseline indicators of genetic diversity and inbreeding. For each population, we estimated F-statistics [60] and allele richness in FSTAT. Confidence intervals for F IS , a measure of intrapopulation heterozygote deficiency due to inbreeding, were calculated using the boot.ppfis function of hierfstat package v0.04-22 in R [61]. Mean allele richness was determined using a rarefaction method [62] in PopGenReport [63], implemented in R. We used FSTAT to compare populations with respect to allele richness and F-statistics, using 1000 permutations.
Detection of population bottlenecks. We used BOTTLENECK v 1.2.02 to test for significant heterozygosity excess, applying a one-tailed Wilcoxon test with 1000 iterations, using the two-phase model (TPM) (90% stepwise mutations, variance 10) [64]. A mode-shift test for distortion of the allele frequency distribution [65] was also implemented in BOTTLENECK. We derived M-ratio statistics [66] using the mRatio function of the strataG package v 2.0.2 [67], implemented in R. Effective population size (N e ). Two single sample methods were used for estimation of N e . The linkage disequilibrium method [68,69] was implemented in NeEstimator ver 2.1 [70], assuming random mating, deriving confidence intervals by jack-knifing (1000 iterations). We also used the sibship assignment method [71], which estimates the current effective breeding size of the population, implemented in COLONY 2.0.6.3 [72], using the same input parameters as detailed below for sibship and parentage analysis. Confidence intervals were obtained by bootstrapping.
Further investigation of breeding between relatives. We applied the inbreeding function of adegenet, version 2.0.1 [59] in R, to derive genetic estimates of the pedigree inbreeding coefficient FPED, which denotes the probability that both alleles at a single locus are identical by descent from a single ancestor [73]. Pairwise relatedness (Rxy) was estimated using a maximum likelihood method in ML-Relate [74]. To calibrate Rxy values with first degree family relationships in our data, we simulated genotypes for pairs of individuals with defined relationships (100 pairs for each category of unrelated, half-sibling, full-sibling and parent-offspring), using the familysim function of the package related v1.0 [75], implemented in R, based on allele frequency data of observed datasets. Pairwise Rxy between each pair of simulated genotypes of defined relationships was measured using ML-Relate. Means and confidence intervals were derived in R. Significance of within-population relatedness was tested using the grouprel function of related v1.0.
For sibship and parentage analysis we used a full-likelihood method, implemented in COL-ONY 2.0.6.4 [72], assuming both male and female polygamy. In the absence of known pedigree structure, all individuals for each study population were treated as a single offspring group. We used default settings for sibship priors, including small sibship size, with the aim of reducing false sibship assignments [76]. The outputs of three independent replicate runs, using independent seeds for random number generation, were examined to confirm convergence to the same configuration and log likelihood. The best maximum likelihood cluster configuration was used to infer half-and full sib dyads and inferred parentage.
Interpretation of results. We compared our results for F-statistics and allele richness from UKAGP study populations with summary statistics from two published studies of adder populations in mainland Europe [44,77], and from a site containing a very large (n >500) population of lowland adders in northern Belgium ( [78], Mergeay & Bauwens unpublished data). Direct statistical comparison was precluded by only partially overlapping microsatellite panels between the studies (2/8 loci of our study were in common with Ursenbacher et al [44,77]; 3/8 loci in common with Bauwens et al [78]).
In the UK, the MTAC survey demonstrated opposite average population trends between sites with small and large adder populations, the threshold being a mean normalised peak count of 10 adders, below which there was significant decline over time [34]. This provides an approach of demonstrated demographic relevance with which to classify and compare populations according to their likely risk of decline. At eight sites in our study, more than 10 individuals had been sampled on a single visit (range [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. The number of adders sampled at the other eight sites was lower (range 4-10) (Table 1), despite equivalent or higher sampling effort by an experienced ecologist familiar with the sites. We therefore applied an equivalent threshold to categorise UKAGP study populations as large (presumed lower risk) (count >10, n = 8), or small (presumed higher risk) (count �10, n = 8) (Table 1), based on the number of adders sampled, as an approximation of the MTAC criteria. We compared first-line summary statistics, FPED and R xy between small and large populations defined in this way, using a Wilcoxon rank sum test implemented in R.
Although exceeding the MTAC peak count threshold, the springtime counts for the UKAGP populations in the "large" category were still relatively low, with a maximum of 27. In the absence of a very large well-characterised UK population for comparison, we analysed 50 genotypes from the Belgian site [78], focusing primarily on indices of relatedness and inbreeding. This 1570 ha site ("Groot Schietveld", N 51˚20-22'-E 4˚32-37'), has been used as a military exercise zone since 1893, and is separated from neighbouring adder populations by a minimum of 18km of unsuitable agricultural habitat. The site is transected by a road, constructed in 1875 as a narrow cobble road, then transformed in 1982 to its present state as a two-lane asphalt trunk road, flanked by two asphalt cycle tracks (total width ca. 20 m). In capture-mark-recapture studies, the road was shown to constitute an effective hard barrier to adder migration (Claus & Bauwens, unpublished data). The sample of 50 Groot Schietveld genotypes (GS50) was selected randomly from the total dataset of microsatellite genotypes of individual adders, which had been sampled from 14 locations within the Groot Schietveld site between 2011 and 2013, and genotyped for 9 microsatellite loci as detailed in Bauwens et al 2018 [78]. The sample included 18 individuals from the north east (NE) segment of the site, and 32 from the south west (SW), relative to the transecting road. No mtDNA sequence data were available for the GS50 sample.

mtDNA haplotypes
Eight different Cytb/CR haplotypes were identified in 97 individuals across the 16 sites (Fig 1). Six haplotypes were unique to single sites. Three sites had more than one Cytb/CR haplotype. All 53 individual samples sequenced from a cluster of ten sites in the West Midlands/South West had the same Cytb/CR haplotype (WMids Haplogroup). As predicted, the UKAGP haplotypes broadly clustered with the mainland Europe Northern phylogenetic clade of V.

Microsatellites
Quality control. There was no evidence for null alleles, allele drop out or linkage disequilibrium between loci. No two samples had identical genotypes. Results from duplicate samples confirmed consistency of genotyping. Samples that failed to amplify at a minimum of 6/8 microsatellite loci were excluded from analysis. In all, 186 samples (84.5%) were retained for downstream analysis. There was no significant divergence from HWE. Two small populations TM (n = 5) and CC (n = 4) had missing data for more than one individual at a single locus (CA71 and Vu4 respectively). For these, genetic diversity was calculated with the omission of the relevant locus, which had little impact on summary statistics for the other populations (S2 Table).
Genetic substructure of WMids Haplogroup. The WMids Haplogroup includes ten geographically neighbouring sites, separated by up to 100km (Fig 1, Table 4). All individuals tested in these sites shared the same mtDNA haplotype, consistent with their origin from a common ancestor. We therefore investigated this group for evidence of more recent differentiation, using microsatellite markers. Although this haplotype was found two in other study sites, they were not included in this analysis, as both showed evidence for additional haplotypes, and were from less intensively sampled regions, with larger distances between sites.
Genetic differentiation within the WMids Haplogroup is reflected in the pairwise F ST matrix between these sites (Table 4). Mantel testing for isolation by distance was negative (r = 0.0042; simulated p value = 0.538). In STRUCTURE there was support for hierarchical clustering [56], both with and without applying the locprior option (Fig 2; S2 Fig). DAPC analysis also showed clustering within the WMids Haplogroup, optimal at four clusters (Fig 2). To evaluate concordance between STRUCTURE and DAPC results we also applied the dapc function for three and six clusters, the optimum values for K in STRUCTURE using the locprior option. The population proportional assignments to each cluster by each of the two methods was very similar at the higher hierarchical level of clustering (K = 3), but more divergent for K = 6 (S3 Fig).
Genetic diversity and inbreeding (F IS ). F IS values did not differ significantly from zero in any of the UK study populations (Table 5). This is in line with the findings of the mainland Europe study, where only 2/16 sites had been reported to have raised levels of F IS , in both cases attributable to high homozygosity at a single locus [44]. Estimates of genetic diversity (H S or H E ) were at similar levels across the study sites, and broadly equivalent to published results from European populations [44,77], including the large Belgium lowland population [78] ( Table 6). Allele richness was in a similar range to that of the large lowland Belgian population, and of the Belgian, Netherlands and northern France lowland populations in the studies of Ursenbacher et al [44,77], although with the caveat of only partially overlapping microsatellite panels. We found mean allelic richness was lower (p<0.05) in small populations. This was the only statistic for which small and large populations differed significantly, other than for size (Table 7).
Effective population size N e . The single sample LDNe method [69] failed to deliver plausible results, which may reflect small sample size and high levels of relatedness [79]. Results derived using the sibship assignment method are shown in Table 8. The small populations again gave very wide confidence intervals. Results for GS50 SW and NE by the sibship assignment method are also shown, generating lower results than expected for the very large number of adders on the site, discussed further below. Unfortunately, the combination of small sample size and high relatedness thus prevented us from deriving reliable estimates for effective population size, a significant disadvantage in the study of wild populations.  Population bottlenecks. As small sample size may give rise to false positives in bottleneck tests, we tested different sample sizes randomly selected from the simulated population of 100 unrelated pairs. A sample of 5 genotypes generated a positive result for both heterozygosity excess and modal shift. Results for both tests were negative for simulated samples of 10 or 20 (Table 9). In two of the UKAGP study populations with a sample size of �10, both heterozygosity excess and allele frequency modal shift tests generated positive bottleneck results (EH, WF). Two (BM, CH) was positive for modal shift only. MRatio results did not discriminate between any of the simulated or study populations, irrespective of size, or results of other bottleneck tests (Table 9).
Further investigation of breeding between relatives. Estimated inbreeding coefficient FPED results showed little variation between populations, irrespective of size ( Table 7). The mean population FPED was 0.233 (95% CI 0.218-0.248), although some populations had individual outliers with FPED >0.50 (Fig 3). We found a similar pattern of FPED results for the GS50 sample, again with occasional outliers FPED >0. 50. FPED results did not differentiate between simulated populations of 100 pairs of defined relationship, whether unrelated, half-or full sibling (Fig 3).
Intra-population mean Rxy estimates ranged from 0.135 to 0.377 (mean 0.220, 95% CI 0.188-0.252), with no significant difference between large and small populations (Fig 3;  Table 7). All UKAGP study populations showed significant within-population relatedness (S3 Table). For the GS50 sample, the mean Rxy between individual samples from the same side of the transecting road (intra-SW, intra-NE) was significantly higher than for pairwise Rxy across the road (SW-NE) (Fig 4). As expected, and in contrast to FPED, mean Rxy differed significantly between simulated pairs of defined relationship (Fig 3).
Rxy frequency distribution curves of simulated pairs show clear differences between the defined relationships, with a dominant density peak at Rxy = 0 in the unrelated pairs, and a progressive right shift of the curve for the simulated half-sib and full-sib pairs (Fig 5). The Rxy frequency distribution curves for the individual UKAGP study populations showed variable right shift of the curve, in association with and blunting or loss of the unrelated peak at Rxy = 0. An equivalent pattern was apparent in the GS50 SW and NE populations (Fig 5). The GS50 populations were analysed separately, as mean Rxy is affected by genetic structure within a sample.
Family structure and parentage analysis in COLONY. Fig 5 illustrates networks of individuals linked by first degree relationships (inferred full or half-sibship) in representative UKAGP study populations, based on the best maximum likelihood configurations in COL-ONY, shown with their respective Rxy frequency distribution curves. Rxy distribution curves for simulated unrelated, half-sib and full-sib pairs are shown for comparison. These networks are characterised by extensive linkage at the inferred half-sib level in all populations. In some study populations, especially those with a right shift of the Rxy frequency distribution curve, networks show dominant inferred half-sibships, some very large, sharing the same inferred parent.
When the GS50 sample was analysed as a single group in COLONY, the patterns of inferred parentage differed between individuals from the SW and NE sampling sites (Fig 4), providing further evidence that the dividing road acts as a barrier to gene flow. The network of COL-ONY-derived sibships in the GS50 NE and SW samples also showed a loose pedigree linked at the half-sib level, with occasional larger inferred shared-parent sibships, and rare inferred fullsib dyads (Fig 5).

Concordance between COLONY and DAPC cluster membership.
We found no evidence for intra-population substructure or admixture on STRUCTURE analysis of the individual UKAGP populations (not shown), nor in the GS50 SW or NE samples, consistent with the low F ST results for individual sites. By contrast, DAPC analysis revealed clustering within the all individual study populations, including GS50 (S4 Fig), despite the low F ST results. We hypothesised that patterns of relatedness might underlie this within-population clustering. We

Discussion
The aim of our study was to investigate the genetic status of lowland adders in the UK, in response to concerns about declining numbers, especially affecting small, fragmented populations. We initially adopted a standard panel of microsatellite-based summary statistics, including genetic diversity and the standard F IS measure of inbreeding, to allow comparison with published studies of adders in mainland Europe. We also applied the MTAC criterion of a threshold peak count to categorise study sites into small or large, predicted to be at high or low risk of decline respectively [34]. This initial panel of genetic tests generated a similar pattern of results across all the UKAGP study sites, irrespective of size, although there was a modest decrease in mean allele richness in small populations relative to large. This is likely to have been influenced by the inevitably small sample size, illustrating the difficulty inherent in analysing the unavoidably small sample sizes of the most vulnerable populations. The interpretation of genetic results requires a biologically relevant comparator, especially for single time-point samples. Estimates of allele richness in the UKAGP study populations were similar to those of the large Belgian population [78], and of populations of lowland adders in Belgium, NE France, and the Netherlands [77], despite the only partially overlapping panels of microsatellite markers used in the different studies. The mtDNA haplotypes of our study populations places them within the Northern phylogenetic clade of European adders [37]. They are therefore likely to have been part of the same post-glacial recolonization process as their counterparts in north-eastern France, in which a central-marginal effect on genetic diversity has been described [77].

PLOS ONE
Notably, F IS levels were not significantly above zero in any of the UKAGP study populations, irrespective of size, indicating retention of heterozygosity, despite the clear risk of inbreeding depression in fragmented, isolated populations. Our findings are again consistent with published studies of remnant adder populations of varying sizes in mainland Europe [44,77]. In small isolated populations of species with high philopatry and low vagility, like the adder, mating opportunities will inevitably be biased towards relatives [80,81], and total avoidance of consanguineous mating could rapidly lead to population extinction [82]. We therefore investigated other measures of consanguinity, a very valuable comparison for the effect of size being provided by the Belgian adder population, isolated but thriving, well-studied, and very much larger than the UKAGP study populations [78].

Identity by descent: Genetic legacy of population demographic history
In all the UK study populations, genetic estimates of the mean inbreeding coefficient FPED were at a level consistent with the degree of genetic identity expected in the offspring of a oneoff mating between half siblings [73]. However, the nature of the defined relationship of simulated pairs of genotypes, whether unrelated, full-sibling or half-sibling, did not appreciably influence FPED estimates. This is in keeping with FPED estimates being a measure of cumulative identity by descent, the embedded genetic legacy of long-standing consanguineous breeding. In the WMids Haplogroup, for example, the shared mtDNA haplotype provides evidence for a historic common ancestor, which could date back to as early as post-glacial colonisation. The Belgian site similarly represents a relict population, probably isolated for more than a century [78].

Identity by state: Reflection of contemporary relatedness
By contrast pairwise relatedness in simulated pairs of genotypes was significantly and predictably influenced by the defined relationship, providing a very useful template against which to interpret Rxy results for the study populations. We found a range of patterns of relatedness, with variable loss or blunting of the modal unrelated peak seen in the simulated unrelated pairs of genotypes, and right shift of the Rxy distribution curve seen in the half-or full-sibling simulated pairs. Pairwise Rxy thus generated the most informative results in our study, with the potential caveat that results may be influenced by cryptic genetic differentiation within the sample. Repeat studies will be necessary to determine trends and timescales in changing patterns of relatedness.
COLONY results, like Rxy, provide a snapshot estimate of contemporary relatedness between individuals within the sample, but in the format of best maximum likelihood combinations of inferred genotypes of sibling and parent-offspring dyads. This approach generated dramatic networks of inferred half-sibships in our study populations. However, in the absence of pedigree data to inform COLONY analysis, loose networks of inferred half sibships may simply reflect identity by state, rather than true first-degree relationships, especially for inferred half sibships in larger populations [76]. While we sought to minimise this phenomenon by applying stringent parameters for sibship assignment in COLONY, it is likely to have been exacerbated by the presence of a high level of background relatedness, as well as by the relatively limited number of genetic markers. This high level of relatedness is also the likely explanation for the unexpectedly low estimates of Ne by the sibship method for the GS50 populations.
We found the membership of inferred dominant large sibships in small populations to be concordant with the assignment of individuals to clusters in DAPC, in keeping with DAPC clusters reflecting allele frequency patterns driven by a polygynandrous mating system in a small population, rather than discrete panmictic subpopulations. An equivalent phenomenon of clustering in DAPC, but not STRUCTURE, has also been described in the Prairie rattlesnake (Crotalus viridis) [83]. This is an important consideration when DAPC is used to investigate cryptic genetic structure within consanguineous populations.

Inbreeding depression: Protective effect of polyandry
As a single year snapshot, our study is only a starting point, and inevitably limited by population size, and thus sample size, especially in the study of the most vulnerable populations. We have nevertheless exposed a previously undocumented degree of consanguinity in wild adder populations, despite their showing no increase in homozygosity to suggest inbreeding depression. In models of inbreeding, a system of half-sib mating is more likely to maintain heterozygosity than one of maximum avoidance of inbreeding [84,85]. Polyandry, which is widespread in taxa of live-bearing snakes [86], including the adder [87][88][89], may thus represent a protective mechanism against inbreeding [82,90]. Interestingly, Mourier et al (2013) described a very similar network pattern of extensive pairwise relatedness between individuals of the sicklefin lemon shark (Negaprion acutidens) in French Polynesia [91]. Despite the biological differences, there are clear similarities in the reproductive ecology between the taxa, the lemon shark also being viviparous, of low female fecundity, with a polyandrous mating system and limited distribution [91].

Inbreeding depression: Size matters
While polyandry may delay, it will not prevent the eventual loss of heterozygosity in isolated populations, where movement of adders is prevented by loss of connectivity between patches of fragmented habitat. It is therefore interesting to compare the large, thriving Belgian site, with high levels of relatedness but no loss of heterozygosity, with the well-documented Swedish population of adders with unequivocal evidence for inbreeding depression [14][15][16]. Both sites have been isolated for more than a century by agricultural landscape, are situated some 20km from neighbouring adder populations [15,78], and both would fulfil the MTAC definition of a large population [34]. However, at 1570ha, the Belgian site is significantly larger than its 20ha Swedish counterpart, which provides the likely explanation for the difference in inbreeding. As the Belgian site has more recently been subjected to asymmetric fragmentation by the truncating road, it will be especially interesting to monitor the genetic status of the smaller NE fragment in comparison with its larger SW counterpart unless gene flow can be restored across the truncating road.

Safety in numbers: Demographic vulnerability of smaller adder populations
"Unfortunately, the best way to find tipping points so far has been to cross them-a dangerous proposition" [92]. Our findings suggest that genetic factors are unlikely to be the direct cause of the observed decline in small populations of adders. Instead, small populations may already be "doomed to extinction by demographic factors before genetic effects act strongly" [93], representing the "living dead" [94,95], where continuation of a population or metapopulation becomes demographically impossible. For example, the reproductive ecology of adders renders small populations profoundly vulnerable to stochastic sex bias [96][97][98]. In adders, males are the actively mate-seeking sex [99,100], with only a short interval of female sexual receptiveness, and thus limited time available for mating, while female adders have low lifetime fecundity, with high fitness costs of reproduction [101,9,10], suggesting particular vulnerability of small adder populations to a limiting number of females. Conversely, a relative reduction in males would be predicted to impact on any protective effect of polyandry against inbreeding. Sampling over a limited period precluded an accurate field assessment of sex ratios in our study, but we are addressing this question of breeding sex ratios in ongoing work.

Future studies
We are using radio-telemetry of adult snakes to inform habitat management, especially with respect to connectivity [100], in combination with ongoing genetic monitoring of study populations. This will allow us to investigate the reproductive success of potential mating connections, generating very interesting data with respect to the breeding system in this secretive species, including assortative mate choice, overt or cryptic. Pedigree information will also enhance the interpretation of results in COLONY, including estimates of the effective numbers of breeding adults.
We plan to use genomic sequencing to increase the number of informative markers available, especially important for genetic diversity and pedigree studies. In addition, the development of a genomic SNP panel will help to increase consistency and comparability across sites and laboratories, currently hampered by the limited numbers of microsatellite markers, and potential inconsistencies in their application [102][103][104]. This will be especially important in decision making and post-release monitoring in any future adder translocations, whether for reasons of mitigation or conservation [105][106][107]. A panel of genomic SNPs will also allow the investigation of heterozygosity affecting different loci [104], facilitating the study of inbreeding at the whole genome level. In addition, the pattern of runs of homozygosity [108][109][110][111][112] will help to elucidate the demographic history of this fascinating species, as well as the identification of potential targets of selection [113].

Conclusions
Our results suggest that the most immediate threat to small adder populations is demographic rather than genetic. For larger populations high levels of relatedness indicate that genetic factors are likely to represent a real threat, albeit less imminent, but also less visible and thus more insidious [22]. Continuing monitoring will be essential to determine the urgency and nature of intervention. Our study thus underscores the need for a systematic, evidence-driven approach in conservation planning for adder populations, whether healthy or declining, integrating population genetics and traditional ecology [33]. In this the "true cost of loss and degradation of habitat" [13] should not be neglected, including public engagement to reduce persecution by changing the public perception of snakes [114]. Attempts at genetic or demographic rescue may be similarly doomed unless such underlying factors can be addressed [115].