Visitors   Views   Downloads
PEER-REVIEWED

Introduction

There is general agreement among scientists that we are in the midst of the sixth great mass extinction. Anthropogenic pressure is affecting natural environments both directly or indirectly (Wake & Vredenburg, 2008). Amphibians are amongst the most severely affected; being the vertebrate group catalogued globally at risk (Hoffmann et al., 2010), effective in situ conservation efforts, such as habitat and ecosystem protection or restoration, are urgent and should be prioritized. Ex situ measures, such as captive breeding, should be viewed as a last resort in species recovery (Snyder et al., 1996). Captive breeding programs are typically associated with many limitations: high economic costs, adaptation to captivity (Araki, Cooper & Blouin, 2007; Williams & Hoffman, 2009), poor success in reintroductions (Harding, Griffiths & Pavajeau, 2015; McCleery, Hostetler & Oli, 2014) and more (see below) have been reported. The need for such programs is justified only when it is essential for species’ survival, and should be implemented only after a careful evaluation of costs and benefits of all conservation alternatives (Dolman et al., 2015; McGowan, Traylor-Holzer & Leus, 2016). Often, primary threats such as habitat loss, disease, or overexploitation lead to small isolated populations, which in turn become highly susceptible to additional stochastic threats that can lead to population decline and, eventually, to extinction (often referred to as the ‘extinction vortex’; IUCN/SSC, 2014). Without excluding other in situ actions, captive breeding can play a crucial role in the recovery of some species for which effective alternatives are unavailable in the short term (Frankham, 1995; Griffiths & Pavajeau, 2008; Ralls & Ballou, 1986; Snyder et al., 1996; McGowan, Traylor-Holzer & Leus, 2016).

When a captive breeding program is recommended, the number of initial founders and the correct selection of individuals to maximize genetic diversity are key factors, as they will determine the genetic diversity present in the captive population (Frankham, Ballou & Briscoe, 2010; Mace, Pemberton & Stanley, 1992; Senner, 1980; Willis & Willis, 2010). Captive populations can often only represent a small subset of a natural population due to spatial or financial limitations of the organizations that manage them (Willoughby et al., 2015). Moreover, captive populations are frequently established when the wild populations of the species are at risk or have already suffered significant reductions in population size (Conde et al., 2011; Lyles & May, 1987; Traylor-Holzer, 2011; Williams & Hoffman, 2009), which limits the number of indivduals that can be collected to establish captive programs. For these reasons, captive populations may represent only a small fraction of the genetic diversity of the wild population (Ballou & Lacy, 1995; Mace, 1986; Ralls & Ballou, 1992), limiting the genetic diversity of the ex situ founders. If this is the case, the captive populations may experience a loss of genetic diversity and the accumulation of mildly deleterious alleles, which could lead to inbreeding depression, thereby jeopardizing the long-term viability of the captive populations (Lacy, 1997; Willis, 1993). Another problem associated with ex situ breeding programs is genetic adaptation to the captive environment. Captive environments differ from wild ones, and therefore the genetic variants favoured in captivity differ somewhat from those favoured in natural environments. Crucially, the genetic variants selected for captive conditions may be overwhelmingly disadvantageous in the natural environment, reducing reproductive fitness of captive individuals when released into wild environments (Frankham, 2008). It is then imperative that a population’s genetic adaptation to captivity be considered when managing breeding stocks (Ballou & Lacy, 1995; Frankham, 1995; Frankham, 2008; Princée, 1995; Williams & Hoffman, 2009). Carefully planned breeding programs are therefore required to minimize these potential threats.

Two main aspects must be considered when establishing a breeding program. First, there is broad agreement among wildlife managers that the retention of the 97.5% of the genetic diversity from wild populations should be a sufficient goal for the long-term viability of the founder populations (Edwards et al., 2014). However, in some cases, these values cannot be achieved without compromising the wild populations (Lacy, 1994; Frankham, Ballou & Briscoe, 2010). Under these circumstances, a genetic evaluation and an efficient selection of the founder stock is necessary to ensure the conservation of the highest levels of genetic variability possible. Second, it is estimated that a “self-sustaining” ex situ population must keep a minimum of 90% of the founder’s population genetic variability for the complete duration of the ex situ program; equivalent to the time required for habitat recovery in nature (Ballou et al., 2006; Soulé et al., 1986). Therefore, the capture of the maximum genetic diversity from the wild populations and the long-term conservation of a relatively high level of this genetic diversity within the captive populations are the two main challenges for a successful ex situ breeding program.

Genetic diversity is measured in terms of allelic diversity and heterozygosity (Ballou & Foose, 1996; Ballou et al., 2010). Allelic diversity may represent the evolutionary potential of a population and is an important value to consider for a population’s long-term viability; it plays an important role in the adaptation to environmental change (Ballou et al., 2010; Briscoe et al., 1992; Frankham, 2003; Princée, 1995). Allelic diversity may be affected by population bottlenecks, such as those that can arise as a result of the establishment of breeding stocks (Amos & Balmford, 2001). Heterozygosity can be characterized either as the observed or expected (i.e., according to expectations of Hardy-Weinberg Equilibrium; HDW) heterozygosity. The former is defined by the proportion of genetic loci for which the average individual in a population is heterozygous. The latter is the probability that two homologous genes randomly drawn from the population are distinct alleles (i.e., the mean heterozygosity that would exist in a population if it were in Hardy-Weinberg equilibrium: Höglund, 2009; Lacy, 1994). One of the key questions to be addressed when establishing a breeding stock is “how many individuals are needed in order to obtain and maintain the required levels of genetic diversity?” (Hale, Burg & Steeves, 2012). Previous studies have addressed this question by calculating the number of captive individuals needed to adequately represent the genetic diversity within wild populations using rarefaction algorithms (Hale, Burg & Steeves, 2012; Kalinowski, 2005). Overall, the number of founders that start a breeding program should be a balance between the genetic diversity captured from the wild and its sustainability, the number of individuals to be extracted from the endangered wild population without compromising their populations, and the capacity of the program to maintain the number of breeders and descendants (McGowan, Traylor-Holzer & Leus, 2016).

The Montseny brook newt (Calotriton arnoldi) is one of the most threatened vertebrates in Europe, and is classified as Critically Endangered by the International Union for Conservation of Nature (IUCN) (Carranza & Martínez-Solano, 2009). Its current distribution is limited to an area of only 8 km2 in the Montseny Natural Park, NE Iberian Peninsula, Spain. Furthermore, although it is found in seven closely located brooks, its populations are divided into two main sectors along the Tordera river valley, which are separated by unsuitable habitat: the Eastern sector—three brooks (A1–A3); and the western sector—four brooks (B1–B4) (Fig. 1; Valbuena-Ureña, Amat & Carranza, 2013). With an estimated census population size of less than 1,500 mature individuals (Amat & Carranza, 2005), the major threats affecting this primarily water–dependent species are habitat degradation and disturbance, manifested as over-extraction of water for commercial purposes, deforestation, and the existence of tracks and roads that disrupt the continuity of the brooks where it occurs (Amat et al., 2014). Emergent diseases, such as Batrachochytrium dendrobatidis and Batrachochytrium salamandrivorans, have not yet affected this species, but should not be overlooked as they comprise an important threat to many amphibians, including newts and salamandrids (Martel et al., 2014; Obon et al., 2013; Spitzen-van der Sluijs et al., 2016).

Location of the Calotriton arnoldi distribution range and the captive breeding facility.

Figure 1: Location of the Calotriton arnoldi distribution range and the captive breeding facility.

Localities of C. arnoldi represented do not correspond to the exact geographic locations intentionally due to conservation concerns. Source of the Montseny brook newt drawing: Toni Llobet.

Although the species conservation plan is not yet approved, some in situ management measures have been developed by the Barcelona Provincial Council, including habitat restoration and field monitoring (Amat et al., 2014). The IUCN guidelines recommend that captive breeding programs start when a species still numbers in the thousands. Considering the low number of individuals estimated for this species, the declining trend of its populations and the drying out of mountain streams, the need to start a captive breeding program immediately was fully justified. Thus, a provisional captive program was designed in 2007 by the Catalonian Government, in order to delve into the study and protection of this species. This work aims to maintain a genetic reserve that is representative of the wild populations, and to improve in its captive management while waiting for the results of detailed genetic studies (Catalan Government, 2010).

Following the recommendation of Amphibian Ark, and based on the precautionary principle, 20 individuals were captured in order to serve as founders for the captive population. As a result of preliminary data on the genetic and morphological variability between the two sectors (Carranza & Amat, 2005), founders from each sector were kept separately as two distinct breeding lines. Currently, a studbook for this species is being maintained using the Single Population Analysis and Record Keeping System (SPARKS), developed by the International Species Information System (ISIS) and the PMx program. However, it does not include genetic information about captive individuals. After ten years of operation, this program has successfully reared more than 1,700 individuals (Carbonell et al., 2016). Despite the viability of the breeding program, no analyses have yet been performed to determine the levels of genetic diversity among captive individuals, and to evaluate if the founders initially sampled are enough to satisfy the genetic goals of the captive breeding program.

The aim of this study is to use microsatellite loci to analyze the effectiveness of the current provisional breeding program of the Montseny brook newt (established with limited genetic information) and use this data to design the definitive captive stock in order to maximize its potential benefits to conservation. To do this, we evaluated whether the levels of genetic diversity (allelic diversity and heterozygosity) observed in the initial captive stock are comparable to the genetic diversity observed in the wild populations, and if this diversity is sufficient to be maintained over time. We also estimated the optimal numbers of individuals needed to ensure that the captive stocks accurately match the overall genetic diversity of the wild populations over time.

Materials and Methods

Sample collection and genotyping

Tissue samples were taken from 42 individuals (19 founders and 23 descendants) from the captive breeding facility (Table 1), and compared to results previously published from the wild populations (Valbuena-Ureña et al., 2017). Of the 19 founders, nine were from the eastern sector and belonged to population A1, and 10 were from the western sector and belonged to populations B1 and B2 (Fig. 1; see Valbuena-Ureña, Amat & Carranza (2013) for population codification). All F1 descendants were genotyped and are comprised of 17 F1 individuals from the eastern captive line, and six F1 individuals from the western captive line, all born in 2007. These descendants were interesting, as they originated from founder females that were already gravid when sampled from the wild. Therefore, they could exhibit genotypes from unsampled wild males. Due to spatial limitations in the facilities and readjustments of breeding methodology, only F1 experimental matings were designed; F2 is not yet available for genotyping. In total, 26 captive individuals from the eastern breeding line were analysed and compared to 77 genotyped wild animals from this sector (Valbuena-Ureña et al., 2017). Similarly, 16 captive individuals from the western breeding line were genotyped and compared to 83 wild animals from the western sector. The collection of all samples was conducted under the licenses required by the Catalan Government, with the permit numbers SF/298 and SF/469.

Table 1:
Genetic diversities of the captive breeding lines compared to wild population sectors of the Montseny brook newt (Calotriton arnoldi).
Values represent averages across 24 loci.
Group N A Ar PA HO HE FIS
Eastern
Founders 9 3.417 3.420 46 0.574 0.501 −0.088
F1 17 3.750 3.330 55 0.522 0.513 0.012
Total captive 26 4.083 3.380 57 0.540 0.519 −0.022
Wild* 77 4.167 4.112 75 0.544 0.538 0.090
Western
Founders 10 2.917 2.650 28 0.413 0.421 0.072
F1 6 2.458 2.460 24 0.507 0.387 −0.227
Total captive 16 3.083 2.590 30 0.448 0.418 −0.039
Wild* 83 2.724 2.610 38 0.359 0.352 0.184
DOI: 10.7717/peerj.3447/table-1

Notes:

N

sample size

A

number of alleles per locus

Ar

allelic richness

PA

number of private alleles

HO

observed heterozygosity

HE

expected heterozygosity

FIS

inbreeding coefficient

Values in bold indicate statistical significance after Bonferroni correction.

Wild* values were obtained from Valbuena-Ureña et al. (2017).

Tissue samples consisted of tail tips or fingers preserved in absolute ethanol. Genomic DNA was extracted using QiagenTM (Valencia, CA, USA) DNeasy Blood and Tissue Kit following the manufacturer’s protocol. For consistency, all 42 individuals were genotyped for the same 24 polymorphic microsatellite loci used to examine the wild populations (Drechsler et al., 2013; Valbuena-Ureña, Steinfartz & Carranza, 2014; Valbuena-Ureña et al., 2017). Details of the PCR-multiplex combinations of loci and specific amplification procedures are given in Valbuena-Ureña et al. (2017).

Data analysis

MICRO-CHECKER (Van Oosterhout et al., 2004) was used to check for potential scoring errors, large allele dropout and the presence of null alleles. Pairwise linkage disequilibrium between loci and deviations from Hardy-Weinberg equilibrium in each grouping and for each locus were checked using the software GENEPOP version 4.2.1 (Rousset, 2008). Genetic diversity was measured for each grouping as the mean number of alleles (A), observed (HO), and expected heterozygosity (HE) using FSTAT 2.9.3.2 (Goudet, 1995). Formal analyses of raw allelic richness values were not performed, because the number of alleles detected at a locus can be influenced by sample sizes. Therefore a rarefied estimate of allelic richness (Ar) was obtained with HP-RARE (Kalinowski, 2005). The nonparametric Wilcoxon signed-rank test (Wilcoxon, 1945) was used to test for differences in diversity levels between captive and wild groups using locus-specific values of A, Ar and HE as paired replicates. In analyses based on allelic richness, separate rarefactions were performed for each sector grouping, to account for the smallest sample size from a particular group. The observed number of private alleles (PA), defined here as alleles found in a single population throughout the study region, for each locus and each grouping was calculated with GDA (Lewis & Zaykin, 2000), and the ratio of presence compared to wild PA was computed. Further, overall and intrapopulation subdivision coefficients (FIS) for all markers and groupings were calculated using FSTAT 2.9.3.2, which calculates the estimator f of Weir & Cockerham (1984) for each marker, as well as a multilocus estimate. Genetic differentiation between the captive and the wild populations was derived from two parameters, FST and RST (Slatkin, 1995), values were obtained using ARLEQUIN 3.5.1.2 (Excoffier & Lischer, 2010).

Relatedness among individuals, effective population size, and bottleneck detection

Pairwise values of genetic relatedness (rxy) were calculated using two estimators: the rqg89 estimator (Queller & Goodnight, 1989) and the rlr99 estimator (Lynch & Ritland , 1999). These estimated indices were compared following the Blouin et al. (1996) criteria implemented in iREL (Gonçalves da Silva & Russello, 2009), using the population allele frequencies estimated for each sample. The cut-off values were applied following the method described in Blouin et al. (1996), in order to classify individuals into relationship categories: unrelated (UR), half-siblings (HS), full siblings (FS) and parent–offspring (PO). Using this method, there is a greater than zero probability of misclassifying individuals if observed values of relatedness fall outside theoretically expected values (Russello & Amato, 2004). To minimize this error, the cut-off values specific for our samples were calculated using the Monte Carlo simulation procedure implemented in Blouin et al. (1996). In total, 10,000 dyads in each of the four relatedness categories were randomly generated using genotypes and allele frequencies specific for each eastern and western captive population as input data; rqg89 and rlr99 estimates for each simulated dyads were then computed. Expected misclassification rates were computed using the cut-off values described in Blouin et al. (1996) (the midpoints between the means of the distributions of pairwise relatedness estimates of each simulation category). Then, the index at which the overlap in the empirical distribution of unrelated and half-sibling dyads was minimized was selected. The relationship category compatible with the observed rxy value was then determined for each individual pair. The observed distribution of pairwise relatedness among individuals for each population was plotted against those of the four simulated relatedness categories.

The effective population size (Ne) for each breeding line based on one-sample Ne estimator method was calculated with ONeSAMP (Tallmon et al., 2008). ONeSAMP employs approximate Bayesian computation. The analyses were run specifying a minimum Ne value of 2 and a maximum of 100. To detect evidence of a recent bottleneck, we tested for significant heterozygosity excess or deficit with the program BOTTLENECK, version 1.2.02 (Piry, Luikart & Cornuet, 1999). This approach is effective for detecting bottlenecks that have occurred recently, within the past 0.2–4.0Ne generations (Luikart & Cornuet, 1998). The stepwise mutation (SMM) and two-phase mixed (TPM) models were tested with default parameters. We also sought signatures of genetic bottlenecks by calculating Garza and Williamson’s (Garza & Williamson, 2001) M-ratio for each breeding line using ARLEQUIN. The M-ratio is the mean ratio of the number of alleles to the range of allele size. When alleles are lost from a population, the number of alleles decreases at a faster rate than the allelic size range, so small (<0.68) M-ratio values are indicative of populations that have gone through a genetic bottleneck at some time in the past.

Random subsampling of empirical datasets

In order to estimate how many captive individuals would be needed for each population to accurately match the level of the genetic diversity present in wild populations, comparisons of genetic diversity estimates between real and distinct sample-sized simulated populations were conducted following the methodology from Hale, Burg & Steeves (2012). A simulated dataset consisting of 100 replicates for each cluster was constructed, with the following sample sizes: 5, 10, 15, 20, 25, 30, 35, 40, 45 and 50 individuals. Each replicate contained a random subset of individuals from the empirical dataset, created using a macro in Excel. This was designed to assign a random number (between 1 and 10,000) to each individual from the empirical dataset, and then sort the dataset by this number. The macro selects the first five (or 10, 15, etc. depending on the sample size category) to a new worksheet, 100 times, resulting in 100 simulated ‘populations’ that are independent, i.e., random subsamples of the empirical dataset, at each sample size. Empirical datasets comprise all individuals genotyped, including wild and captive specimens. Sampling was done without replacement, so no individual was present more than once in the same replicate (as in a real population genetic dataset). As replicates were independent from each other, the same individual could be present in more than one replicate of the simulated dataset at each sample size. The mean expected heterozygosity (HE) and the mean percentage of alleles—including all alleles, common alleles (frequency >0.05) and rare alleles (frequency <0.05)—detected at each sample size (mean of the 100 random replicates per size) were then compared to each cluster (eastern sector, western sector, cluster A1A2 and cluster B1B2B4). The mean pairwise FST between the 100 random replicates at each sample size and the empirical dataset was also calculated for each dataset. In order to further examine variations between the simulated and empirical estimations, ANOVA with Tukey HSD post hoc tests were performed separately on data from each cluster using Statistica v.7.

Results

A total of 140 alleles were identified in the 42 captive individuals analyzed, compared to the 170 alleles identified in the 160 individuals analyzed from the wild populations by Valbuena-Ureña et al. (2017). Regarding the presence of alleles by sector, 98 alleles were identified in the eastern captive breeding line. Ninety-five of these were present in the wild populations previously examined (making up 72% of alleles from the wild population); whereas three alleles were new. For the western captive breeding line, 74 alleles were found, of which one was new (not previously identified in the wild) and 73 were also present in the corresponding wild populations (78% of the wild alleles). Captive and wild allele frequencies are given in Table S1.

All 24 microsatellites were polymorphic for the eastern captive breeding line, but Us3 and Us7 were monomorphic throughout the western captive line (Table S2). Only Ca8 from the eastern sector was detected as a null allele. However, this does not suggest a locus-specific problem with null alleles, as no locus was affected by null alleles in all samples or in all populations. One locus (Calarn 52354) in the eastern captive line deviated from Hardy-Weinberg equilibrium (HWE), while the others did not show significant departures from HWE after applying Bonferroni correction (p > 0.00104). Linkage disequilibrium was only found between one pair of loci (Ca07 and Calarn52354 in the eastern captive group) after applying a Bonferroni correction (p > 0.00018).

Generally, higher allelic richness values were observed in both the eastern captive line and wild populations, in comparison to their western counterparts (Table 1). When comparing genetic diversity measures of captive individuals to wild ones (Valbuena-Ureña et al., 2017), similar values of mean number of alleles (A) are found in the eastern sector (p = 0.5135), while captive individuals from the western sector showed slightly higher values of A than wild specimens (p = 0.0012). Regarding differences among groups whilst taking into account sample size, rarefied allelic richness was lower in the captive individuals than in the wild ones. The difference was not statistically significant in the eastern sector (p = 0.0675), though it was in the western sector (p = 0.0007). These differences are almost certainly biologically significant, as it is impossible to retain genetic diversity that did not permeate into the captive population. Also, similar values of allelic richness were found among founders and descendants. The number of alleles per locus ranged from one to seven (Table S2). Each breeding line was found to contain private alleles (PA) in only one sector. A total of 57 PA were found in the eastern captive line (two of them not found in the wild populations) compared to 75 present in the wild populations (i.e., 76% PA representation in eastern captive line). In contrast, 81.5% PA representation was found in the western captive line (31of 38 PAs previously found in the wild). Regarding F1, 11 and two new PAs were found in the eastern and western captive breeding lines, respectively.

The mean expected heterozygosities (HE) were 0.519 and 0.418 in the eastern and western captive breeding lines, respectively. On a locus by locus basis, HE of the captive eastern line was similar to the wild populations (p = 0.4750). No differences in HE were observed between captive and wild individuals of the western populations (p = 0.8074), although a tendency towards smaller average heterozygosity was observed in the wild populations. Overall, FIS did not show values significantly different from zero for each breeding line after applying Bonferroni correction. Pairwise FST and RST values revealed significant divergences among captive and wild counterparts. Values of FST and RST between the captive and wild populations was 0.061 and 0.036 ( p < 0.001) for the eastern sector, and 0.044 and 0.032 (p < 0.01) for the western sector.

Relatedness among individuals, effective population size, and bottleneck detection

Comparison of the relatedness indices of rqg89 and rlr99 suggests that the latter performed better at discriminating between unrelated and half-sibling dyads, with fewer unrelated individuals being misclassified as half-sibling dyads following Blouin et al. (1996) (Table S3). All further analyses were therefore based upon the index of Lynch & Ritland (1999). The observed distribution of pairwise relatedness among individuals for each captive population was log normal, with a peak around zero (Fig. 2). Using the empirical distribution derived from 10,000 simulated dyad for each relatedness category, a dyad with a relatedness value of ≥0.1881 and 0.2295 for the eastern and western breeding line respectively, was determined to have a ≤ probability of 0.05 of being unrelated (Table S4). The proportion of dyads classified as unrelated for the eastern and western founder population is 94.44% and 88.89%, respectively (Fig. 3).

Observed distributions of relatedness (pink) for the (A) A1A2 dataset and (B) B1B2B4 dataset, plotted against expected distributions for 10,000 simulated pairs for each of the following relationship categories: nonrelated (green); half-siblings (yellow); full-siblings (red); parent–offspring (blue).

Figure 2: Observed distributions of relatedness (pink) for the (A) A1A2 dataset and (B) B1B2B4 dataset, plotted against expected distributions for 10,000 simulated pairs for each of the following relationship categories: nonrelated (green); half-siblings (yellow); full-siblings (red); parent–offspring (blue).

Proportion of individuals classified as unrelated (dark grey), first order relatives (full siblings and parent–offspring, white), and half-siblings (light grey) for each of the two Calotriton arnoldi breeding line (founders, offspring, total captive).

Figure 3: Proportion of individuals classified as unrelated (dark grey), first order relatives (full siblings and parent–offspring, white), and half-siblings (light grey) for each of the two Calotriton arnoldi breeding line (founders, offspring, total captive).

The estimated Ne for the eastern captive breeding line was 38.96 (median, 95% CI [34.27–52.25]), and Ne = 17.62 (median, 95% CI [15.30–24.05]) for the western captive breeding line (Table 2). Founders Ne estimates from both breeding lines were roughly 10, while effective population sizes estimated from the F1 were 17.03 and 7.97 for eastern and western captive founder stocks, respectively. No captive population showed any signs of a recent bottleneck under the TPM (one tailed test for heterozygote excess) or SMM models (Wilcoxon one-tailed test for excess of heterozygosity; P > 0.05 for all tests). In contrast, M-ratio tests suggest that both breeding lines experienced bottlenecks, with a Garza-Williamson M-ratio of 0.27 and 0.30 for the eastern and western breeding lines, substantially lower than the threshold of 0.68 considered indicative of a bottleneck.

Table 2:
Estimates of effective population size (Ne) for each breeding line calculated with ONeSAMP, with estimations of the upper and lower 95% CI.
Group Ne 95% Cis
Eastern
Founder 10.20 9.20 12.18
F1 17.03 15.07 22.12
Total captive 38.96 34.27 52.25
Western
Founder 10.93 9.38 15.07
F1 7.97 6.93 10.22
Total captive 17.62 15.30 24.05
DOI: 10.7717/peerj.3447/table-2

Size of the captive lines needed to accurately represent the genetic diversity in wild populations

The mean percentages of allele detection were significantly different across the range of sample sizes (P < 0.001, in all four datasets) when considering all alleles, rare alleles, and common alleles (Table S5). In the two former cases, the Tukey HSD posthoc comparisons indicated significant differences among all sample size groups tested. The differences found when considering 20 individuals were not significant when only the common alleles were analysed. This is also the minimum number of individuals needed to detect at least the 99% of the common alleles (Fig. 4). The mean expected heterozygosities were far more precise (lower standard deviation) and accurate with increasing sampling size (ANOVA P < 0.001, in all four datasets), much of the increase occurring when the sample size grew from five to 20 individuals (Tukey HSD posthoc comparisons, Table S5; Fig. 5). Increasing sample size above 20 individuals seemed to have little impact on the mean HE and its standard deviation for all four datasets. As above, the genetic distance between simulated and empirical datasets decreased as the sample size was increased, but again the decrease was not linear. ANOVAs indicated that Pairwise FST between each replicate and the empirical dataset were significantly different between sample sizes (P < 0.001, in all four datasets, Table S5). The difference was mainly from samples 5 to 20–25. No significant differences in pairwise FST was detected from 25 individuals upwards (Tukey HSD posthoc comparisons, Table S5), suggesting that there is little to be gained from sampling more than 25 individuals per captive population. In all four datasets, the 0.5% of genetic differentiation from wild populations is achieved with up to 20 individuals (Fig. 6), showing that the captive stock matches 99.5% of the wild genetic composition.

Mean percentage of all alleles (red), alleles at a frequency >0.05 (green) and alleles at frequency <0.05 (blue) detected at each sample size (of the 100 random replicates per size).

Figure 4: Mean percentage of all alleles (red), alleles at a frequency >0.05 (green) and alleles at frequency <0.05 (blue) detected at each sample size (of the 100 random replicates per size).

The vertical grey lines show the sample size at which at least the 99% of the common alleles (frequency >0.05) were detected. (A) eastern sector; (B) western sector; (C) cluster A1A2; (D) cluster B1B2B4.
Mean and ± one standard deviation of expected heterozygosity (HE) for the 100 random replicates at each sample size.

Figure 5: Mean and ± one standard deviation of expected heterozygosity (HE) for the 100 random replicates at each sample size.

(A) eastern sector; (B) western sector; (C) cluster A1A2; (D) cluster B1B2B4.
Mean pairwise FST and ± one standard deviation comparison between the 100 random replicates at each sample size and the empirical dataset.

Figure 6: Mean pairwise FST and ± one standard deviation comparison between the 100 random replicates at each sample size and the empirical dataset.

(A) Eastern sector; (B) western sector; (C) cluster A1A2; (D) cluster B1B2B4.

Discussion

The main aims of the initial ex situ breeding program of Calotriton arnoldi when first established in 2007 were to maintain reservoir populations that reflected the genetic diversity found in the wild populations, to improve its captive management (produce and select the appropriate specimens for reintroduction), and to continue gathering knowledge on the biology of this species. The success of the program will depend largely on the genetic variability of the captive stock (Lacy, 1994). Decisions based on erroneous or incomplete genetic data may lead to a decrease of genetic diversity, jeopardizing the long-term survival of the captive stock and potentially negatively affecting the wild populations by the reintroduction of genetically depauperate captive specimens (Frankham, 1995; Hedrick & Kalinowski, 2000; Saccheri et al., 1998; Spielman, Brook & Frankham, 2004). This is the first study to evaluate the initial captive breeding program of the Critically Endangered Montseny brook newt. As such, these insights are crucial in ensuring the genetic quality of the breeding stock of the Montseny brook newt, and to successfully establish a long-term self-sustaining captive population. Furthermore, the methodology and conclusions reached herein can be used as an example for the genetic management of other ex situ conservation programs.

Here we discuss our two initial questions which we considered when establishing the ex situ program: (Q1) do we have a good genetic representation of neutral diversity from the wild populations? If not, (Q2) how large should the captive population be to reach this target? It is crucial to detect the point at which the genetic benefits gained by adding additional individuals to the captive populations is outweighed by the increased costs of sampling and maintaining those extra founders.

Q1: estimation of the genetic diversity present in the captive stock

The levels of the genetic diversity measured in the captive dataset indicate that the individuals that form the current captive program are good, but not optimal, representatives of the genetic diversity present in the wild. Genetic theory indicates that the founder population should capture 97.5% of the genetic diversity present in the wild populations (Leus, Traylor-Holzer & Lacy, 2011). The current (initial) captive stock of Calotriton arnoldi captures roughly 93–95% of the total genetic diversity previously observed in the wild. Moreover, the percentage of unrelatedness among dyads from both breeding lines does not exceed 95%. Therefore, in order to increase genetic diversity, new genetic material should be incorporated by introducing new unrelated individuals or their sperm. Since this species present sperm storage (Alonso, 2013) and multiple paternity (Valbuena-Ureña et al., 2017), an ideal strategy to maximize the gene pool of the captive population while minimizing the impact on wild populations would be to incorporate wild-caught gravid females. Thus, descendants may retain the genetic diversity from individuals not kept in captivity.

Another issue to consider is that the genetic diversity of the captive stock should be maintained throughout the entire duration of the breeding program. A drop of the effective population size is detected on the eastern captive line (Ne = 39) when compared to the eastern wild populations (Ne = 86; Valbuena-Ureña et al., 2017), as well as in the western captive line (Ne = 18) compared to the western wild populations (Ne = 80). Although it has been shown that this species is able to survive in the wild with very low values of Ne, the possibility of inbreeding depression in captivity should be considered, since the strategies used by this species to maintain high levels of genetic diversity in the wild despite the relatively low number of individuals are not known, or may not be possible in captivity. Genetic drift erodes genetic variation as Ne decreases (Allentoft & O’Brien, 2010; Montgomery et al., 2000), elevating the probability of fixation of deleterious alleles, and reducing the effectiveness of selection. This may lead to a reduction of overall fitness by limiting adaptive responses (Hare et al., 2011). Animal models suggest that, in captivity, a Ne > 50 is needed to avoid the risk of inbreeding depression (Frankham, Ballou & Briscoe, 2010; Traill et al., 2010). Thus, efforts to enhance the Ne of the captive population of C. arnoldi captive should be undertaken.

In keeping with the wild situation, a higher allelic richness observed in the eastern dataset when compared to the western one (Valbuena-Ureña, Amat & Carranza, 2013; Valbuena-Ureña et al., 2017). Further, the high levels of genetic differentiation observed between the eastern and western wild and captive populations suggests that the maintenance of two separated breeding lines is needed. Although outcrossing the two management units may help to increase the genetic diversity (Frankham, Ballou & Briscoe, 2010), we do not recommend it at this point for two main reasons: (a) the importance of maintaining the two genetically and morphologically differentiated groups; and (b) the potential threat to the existing populations due to outbreeding depression (OD) (Allendorf & Luikart, 2007). Existing conservation guidelines might focus on maintaining high levels of genetic diversity within the two distinct breeding lines. In the case that captive-bred individuals are introduced into brooks currently not occupied by this species, it is important to consider the eastern or western geographical position of the brook, in order to avoid possible future interbreeding as a result habitat reconnection following floods in the rivers of each sector.

Q2: increasing the captive breeding population

Increasing the sample size clearly results in a better representation of alleles and expected heterozygosities between the simulated and the empirical datasets. The genetic distance between these two datasets also decreases when sample size is increased. The probability that a particular allele will be included among the founders is determined by its frequency in the population. Although the presence of rare alleles does not contribute much to the global genetic diversity of a certain species, they a enhance populations distinctiveness and adaptive potential. However, as seen in Fig. 4, a very large number of founders (>40 founders by sector/cluster) would be required to represent 80% of the rare alleles. In order to avoid extra costs associated with maintaining oversized captive stocks, it is crucial to detect the point at which increasing the sample size will not result in a significant benefit in terms of allelic representation (e.g., retaining all common alleles and the maximum number of rare alleles, or maximizing the mean expected heterozygosity of the founders; Lacy, 1994). This sample size is defined as the minimum captive stock size to meet the genetic goals.

Our results suggest that a minimum of 20–25 individuals for each captive breeding line should be kept in order to maintain good representation of the genetic diversity. Our results are in accordance with the Amphibian Ark guidelines (Schad, 2008) and Lacy (1994), which recommend breeding at least 20 unrelated founders by sector, with an equal sex ratio.

The expected heterozygosity (HE) of the eastern captive stock is significantly lower than that of the wild populations (p = 0.0082), indicating that the eastern captive stock only partially represents the gene pool of the eastern wild populations. Captive individuals of this breeding line originated from a single population (A1). As suggested by previous analyses of both mitochondrial and nuclear data, the eastern populations are grouped into two clusters, population A3 being the most differentiated and genetically diverse, and populations A1 and A2 clustering together in the microsatellite analysis (Valbuena-Ureña et al., 2017). Therefore, we recommended to introduce further individuals captured mainly from population A3.

In the western breeding line, in order to reduce the relatedness detected (Fig. 3), and to avoid the possibility of inbreeding, an increase of breeding individuals from this sector is also recommended. According to microsatellite loci data, the western sector clusters into two groups, the main one formed by populations B1, B2 and B4, and the other by population B3, which is separated into a distinct cluster (Valbuena-Ureña et al., 2017). Whereas the overall HE of the western wild populations is 0.352, the heterozygosity increases up to 0.423 when population B3 is excluded. The low genetic diversity (HE = 0.197) detected in this latter population biases the mean HE value of the western wild populations. Currently, the western captive line comprises no founders from B3, but as a result of its low genetic diversity, it is unadvisable to include breeders from this site. Therefore, if additional founders are incorporated into this breeding line, individuals from B1, B2 and B4 should be prioritized.

A need for periodic re-evaluation

Similar values of rarefied allelic richness are found among founders and descendants, which indicates that the descendants inherited almost all the microsatellite alleles from the founders. Eleven and two new PAs were found in the descendant group of the eastern and western sectors respectively; these were not detected in the founders. This could be explained by the fact that some founder females were already gravid when they were caught in the wild, so the genotype of the males was unknown. Also, some females may have fertilized their eggs with sperm from wild males stored in their reproductive apparatus, using a strategy suspected within C. arnoldi (Alonso, 2013) and present in many urodel species (Kühnel, Reinhard & Kupfer, 2010). Therefore, it is highly recommended that the F2, F3, etc. descendants are genotyped in order to ensure that genetic diversity is maintained over time.

Although no signs of a recent bottleneck are detected by BOTTLENECK (more effective at detecting very recent bottlenecks of low magnitude), M-ratio results (best suited for detecting more severe and older bottlenecks) suggested that a past bottleneck occurred in both breeding lines. Thus, it is crucial to keep monitoring for population bottlenecks in managed species since a reduction in Ne may enhance the rate of inbreeding, loss of genetic variation, and fixation of deleterious alleles considerably; thereby increasing the risk of population extinction (Luikart et al., 1998). As the captive breeding program of the Montseny brook newt was established very recently (less than 10 years ago) it is essential that signs of reduced genetic variation or a founder effect are monitored. Therefore, the captive stock and subsequent cohorts should be scrutinised over the long-term in order to effectively preserve genetic variation.

The current breeding protocol widely used by many captive breeding programs to minimize loss of genetic diversity involves pairing individuals according to the mean kinship value (MK) (Willoughby et al., 2015), which is based on the number and degree of relatedness of relatives that an individual has within the captive stock. The use of molecular data has proved useful to investigate the relatedness among individuals, to allow for the design of optimal mating groups, as well as allowing the identification of the natural population from which the founders originated. This underlines the importance of using molecular markers to evaluate genetic management of captive breeding programs. Therefore, further analyses on the resolution of the pedigree are needed to avoid possible mating of closely related individuals, allow deterministic parental assignment of offspring, and design optimal mating groups for maximizing diversity and avoiding genetic erosion in captivity.

Finally, a further problem associated with small founder populations is their significant rapid adaptation to captivity (Araki, Cooper & Blouin, 2007; Frankham, 2008; Gilligan & Frankham, 2003; Griffiths & Pavajeau, 2008; Witzenberger & Hochkirch, 2011). A possible solution would be to occasionally renew the captive populations with previously genotyped individuals from the wild, thereby reducing adaptation to captivity to a minimum (Williams & Hoffman, 2009; Montgomery et al., 2010).

Supplemental Information

Microsatellite allele frequencies in the wild and captive populations of the Montseny brook newt (Calotriton arnoldi)

Sample sizes at each locus are also provided. Data from wild specimens were obtained from Valbuena-Ureña et al. (2017).

DOI: 10.7717/peerj.3447/supp-1

Estimates of genetic parameters for each breeding line and locus

N, sample size; A, number of alleles per locus; Ar, allelic richness; PA, number of private alleles; HO, observed heterozygosity; HE, expected heterozygosity; FIS, inbreeding coefficient. Values in bold indicate statistical significance after Bonferroni correction.

DOI: 10.7717/peerj.3447/supp-2

Rate of missclassification for the related for the relatedness indices of Queller and Goodnight and Lynch and Ritland

Rate of misclassification of unrelated (un) to half-siblings (hs) and first-order relatives (fs, full-siblings and po, parent-offsprings) for the related for the relatedness indices of Queller & Goodnight (1989, rqg89) and Lynch & Ritland (1999, rlr99) based on the midpoint between the distributions of relatedness values calculated from 10,000 simulated dyads of each relatedness category using the population allele frequencies estimated for each sample, as described in Blouin et al. (1996).

DOI: 10.7717/peerj.3447/supp-3

Cut-off values (midpoints between the means of the distributions of pairwise relatedness estimates; Blouin et al., 1996) of each simulated relationship category

Un, unrelated; hs, half-sibling; fs, full-siblings; po, parent-offsprings; rel, related, including hs, fs and po. The relationship category compatible with the observed rlr99 value was and then determined for each individual-pair.

DOI: 10.7717/peerj.3447/supp-4

ANOVA and Tukey HSD posthoc comparisons among the simulated (mean of the 100 random replicates per each sample size) and empirical estimation

ANOVA and Tukey HSD posthoc comparisons among the simulated (mean of the 100 random replicates per each sample size) and empirical estimations of mean percentage of common alleles (frequency >0.05), mean expected heterozygosity (HE) and the mean pairwise FST by each captive dataset (eastern sector, western sector, cluster A1A2 and cluster B1B2B4).

DOI: 10.7717/peerj.3447/supp-5