Genetic variability, management, and conservation implications of the critically endangered Brazilian pitviper Bothrops insularis

Abstract Information on demographic, genetic, and environmental parameters of wild and captive animal populations has proven to be crucial to conservation programs and strategies. Genetic approaches in conservation programs of Brazilian snakes remain scarce despite their importance for critically endangered species, such as Bothrops insularis, the golden lancehead, which is endemic to Ilha da Queimada Grande, coast of São Paulo State, Brazil. This study aims to (a) characterize the genetic diversity of ex situ and in situ populations of B. insularis using heterologous microsatellites; (b) investigate genetic structure among and within these populations; and (c) provide data for the conservation program of the species. Twelve informative microsatellites obtained from three species of the B. neuwiedi group were used to access genetic diversity indexes of ex situ and in situ populations. Low‐to‐medium genetic diversity parameters were found. Both populations showed low—albeit significant—values of system of mating inbreeding coefficient, whereas only the in situ population showed a significant value of pedigree inbreeding coefficient. Significant values of genetic differentiation indexes suggest a small differentiation between the two populations. Discriminant analysis of principal components (DAPC) recovered five clusters. No geographic relationship was found in the island, suggesting the occurrence of gene flow. Also, our data allowed the establishment of six preferential breeding couples, aiming to minimize inbreeding and elucidate uncertain parental relationships in the captive population. In a conservation perspective, continuous monitoring of both populations is demanded: it involves the incorporation of new individuals from the island into the captive population to avoid inbreeding and to achieve the recommended allelic similarity between the two populations. At last, we recommend that the genetic data support researches as a base to maintain a viable and healthy captive population, highly genetically similar to the in situ one, which is crucial for considering a reintroduction process into the island.


| INTRODUC TI ON
Genetic diversity is one of the three classes of biodiversity recognized as global conservation priorities and plays a decisive role in conservation efforts (AZA, 2020; IUCN -International Union for Conservation of Nature, 2019; Sodhi & Ehrlich, 2010). Genetic diversity data on wild and captive populations have shown to be useful for evaluating the consequences of fragmentation and habitat loss, elucidating gene flow and population genetic structure, and defining important evolutionary areas (see Arruda et al., 2017;Gallego-Gárcia et al, 2018;Madsen et al., 2000;McAlileyet al., 2016;Michaelidis et al., 2015;Monzón-Argüello et al., 2015;Wallis, 2019). However, a recent review on 30 years of conservation genetics in New Zealand, Wallis (2019) highlighted that genetic approaches have some limitations related to the definition of minimum viable populations, eco-sourcing, inference of gene flow, and species boundaries. To address these limitations, genetic approaches should grief with other knowledges (for example, demography and reproductive traits) in an integrative conservation effort (Allendorf et al., 2013).
Both ex situ and in situ conservation programs focus on retaining genetic diversity for a minimum period-usually 100 years (Frankham et al., 2008). However, demographic variations associated with captive adaptation, anthropogenic impacts, or the occurrence of inbreeding lead to a decrease in genetic diversity. In several cases, this may carry populations toward an extinction vortex and imminent extinction, as occurred to the isle royal wolves (Canis lupus) (Frankham, 2005;Frankham, 2010;Hedrick et al., 2019;West et al., 2018). Based on this, the management of wild and captive populations has focused on the maintenance of a viable population and its genetic diversity through the avoidance of inbreeding and captive adaptation, and a continuous monitoring of demographic events and genetic diversity (IUCN, 2002;Shafer et al., 2015).
Brazil is considered as a megadiverse South American country due to its high levels of species richness and endemism, with biodiversity hotspot regions (e.g., Cerrado, Rainforest; ICMBio/MMA, 2018a). The last census of Brazilian fauna identified 34 snake species classified as vulnerable, endangered, or critically endangered (ICMBio/MMA, 2018b), most of which are incorporated into national conservation plans. Other than improvement of Brazilian snakes conservation (Navega-Gonçalves & Porto, 2016), only one study focuses on genetic data using Random Amplified Polymorphic DNA (RAPD) in Bothrops moojeni populations (see Dutra et al., 2008).
Even though the access in the island is restricted to the marine and authorized scientific researchers, natural and deliberative bushfires, and anthropogenic actions such as capture of specimens and biopiracy have already been reported (Duarte et al., 1995;Guimarães et al., 2014;Martins et al., 2008). These threats associated with restrict geographic distribution and evidence of populational decline lead this species to be classified as critically endangered (according to the parameters B1ab (iii) + 2ab (iii)) by the International Union for The majority of the studies with B. insularis focused predominantly on ecology (Guimarães et al., 2014;Marques et al., 2012;Martins et al., 2008), reproduction (Amorim et al., 2019;Marques et al., 2013;Silva et al., 2015), and phylogeography (Grazziotin et al., 2006). Populational genetics and molecular parameters of the species, however, remain scarce, even though accumulation of genetic information and maintenance of genetically viable ex situ populations of the species were a component of the aims included in a national plan specific for insular herpetofauna (Bataus et al., 2011), which may allow the design of future reintroduction plans.
Based on this brief historic, due to this species vulnerability, ex situ conservation programs have been designed in Brazilian scientific centers, such as the Instituto Butantan. The ex situ population housed at the Laboratório de Ecologia e Evolução, Instituto Butantan, was established in 2009/2010 with 20 founders (IBAMA no. 25.650-1), aiming to develop a healthy population which could be used for either scientific researches or future reintroduction (Kasperoviczus & Almeida-Santos, 2012). Up to now, breeding couples were designed based on the sperm health and viability for male's selection (Silva et al., 2015), and x-ray analysis to select females in vitellogenesis. Although mating was reported, no study book was maintained for this population.
In this context, knowing that the access of genetic data of wild and captive populations is crucial for conservation programs, our study aims to: (a) characterize the genetic diversity of ex situ and in situ populations of B. insularis using heterologous microsatellites; (b) investigate the existence of genetic structures among and within populations; and (c) provide data for the conservation program of the species.

| Sample collection
A total of 80 samples from representatives of B. insularis ( Figure 1) were used in this study: 49 specimens belonging to an ex situ population maintained at the Laboratório de Ecologia e Evolução (LEEV), Instituto Butantan (São Paulo State, Brazil; sampling area 1; Figure 2a; Appendix S1), 31 specimens belonging to the in situ population from Ilha da Queimada Grande (São Paulo State, Brazil; Figure 2b; Appendix S1), 12 of them were taken in sampling area 2, and 19 in sampling area 3 (Figure 2b). We are using sampling area terminology because of the low accuracy of GPS data in the island, so the centroid point from each area was used as the geographic coordinate of the island representatives.

| Genetic analyses
Details and the script used for the analysis performed in the software R v.3.5.1 (R Core Team, 2018) are available in the Appendix S3. The following sections will focus on reporting the parameters, packages, and software used. Critical probability values (p ≤ .05) for all tests described herein were adjusted with a sequential Bonferroni correction (Rice, 1989).

| Genetic diversity indexes
The number of alleles (A) and the number of private alleles (PA) were calculated using POPPR. Allelic richness was obtained using HIERFSTAT v.0.44-22 package (Goudet, 2005). Additionally, the observed and expected heterozygosity indexes were inferred through the basic package implemented in R.

| Inbreeding indexes
The literature presents two different inbreeding coefficients that possess distinct biological meanings. The inbreeding coefficient (F is ) evaluates the mating system in the population, and the inbreeding coefficient of the genogram or per kinship (f is ) assesses the individual probability of exhibiting an allele that is identical by descent (IBD) due to its parental relationships (Templeton, 2011).
Since both biological meanings offer different insights for conservation, we investigated the inbreeding coefficient for the mating system (F is ) using the FSTAT v.2.9.3.2 software (Goudet, 2002), and the average individual inbreeding coefficient per kinship (f is ) in each population using ADEGENET. Significance of each index was obtained using bootstrap, with 1,000 permutations.

| Bottleneck in the in situ population
The Bottleneck software (Piry et al., 1999) was used to test for a recent demographic reduction in the in situ population. The analysis was performed with 1,000 iterations under two models: the infinite allele model (IAM) and the stepwise mutation model (SMM). Sign up and Wilcoxon tests were used to infer recent bottleneck signals, through significant heterozygosity excess.

| Populational genetic structure
Offspring representatives from the ex situ population were excluded of the genetic differentiation and structure analyses in order to avoid bias.
Populational differentiation coefficients are sorted in three classes: (a) the classical Wright's F st , (b) the standardized analogues F st , and (c) D est . Though none of these indexes are considered an ideal summary statistic, combining them may increase accuracy in attempts to elucidate demographic and genetic population structure (Meirmans & Hedrick, 2011). Wright's F st was calculated using FSTAT. Hendrick's G st ' and Jost's D est indexes were obtained using MMOD (Winter, 2012). AMOVA was performed with POPPR package. Significances of each index and AMOVA were obtained using bootstrap with 1,000 permutations.
For the DAPC analysis, the number of clusters (K) was estimated from one to ten using the lowest value of the Bayesian information criterion (BIC) retaining 30 principal components (PCs) using ADEGENET. We increased the numbers of possible K to 10 to improve the visualization of the plot obtained for BIC values. Subsequently, we conducted an initial DAPC analysis with parameters of 25 PCs, three discriminant analyses, and assuming K number of genetic clusters obtained in the step above. Once bias became possible as a result of the excess of PCs used, the correct number of PCs (i.e., six) was accessed, and the results were validated via the a-scored method.

| Genetic and breeding management of the ex situ B. insularis population
Kinship coefficient values were inferred for 37 captive individuals, 34 of which were alive at the final stage of this study, through the ML-Relate software (Kalinowski et al., 2006). Based on these val-

| RE SULTS
After data quality, linkage disequilibrium and null alleles tests have been performed, a total of 21 loci and four individuals had to be excluded from the study (Appendix S4) because heterologous amplification was not successful in eleven loci; excessive mistyping and missing data (>30%) were found in three loci and four individuals; evidence of linkage disequilibrium and null alleles were recovered in two and five loci, respectively. Therefore, populational analyses were conducted on the basis of a genetic matrix (Appendix S5) con-

| Genetic diversity
The final 11 microsatellites showed similar values of genetic diversity indexes between the populations (Table 1)

| Inbreeding
We obtained negative and significant average values of F is for both Table 1). However, both populations showed positive f is values, which were only significant in the in situ population (f is ex situ = 0.05/CI = -0.42 < x < −0.02 and f is in

| Bottleneck
Under the models tested, the results obtained showed no evidence of a recent bottleneck in the in situ population after Bonferroni sequential correction (IAM: 0.03; SMM: 0.55). However, it is worth mentioning that under the IAM model, nine out of 11 loci showed heterozygosity excess.
However, some loci show moderate to high values of genetic differentiation indexes (for instance, Bmat_070)-it seems to be related to both the existence of private alleles (Table 1)

| Genetic and breeding management of the ex situ B. insularis population
We recovered a mean kinship value of 0.20 for the captive population based on the weighted average of the mean relatedness values from the 34 alive individuals (Appendix S6). Furthermore, the genetic data allowed the two cases of uncertain paternity and uncertain maternity to be resolved. Due the decease of the putative father ID0001 and absence of any tissue for accessing any genetic information, we were not able to test his paternity. Then, paternity analyses were only conduct using the possible father ID011F, revealing relatedness values (R) equal to zero among this possible father and the three offspring (ID06FF, ID08FF, and ID09FF; Table 3; Appendix S7), excluding this possible parentage. Maternity analyses revealed a founder female as mother for each offspring (Table 3; Appendix S7): mother ID0009 to the offspring ID011F, mother ID0010 to the offspring ID013F, mother ID0010 to the offspring ID014F, mother ID0006 to the offspring ID015F, and mother ID0013 to the offspring ID021F. Individuals ID06FF, ID013F, and ID021F were excluded from Appendix S6, given that they are dead, therefore unable to contribute as relatives to the next generations. After addressing these uncertain parental relationships, a pedigree of the ex situ population was drawn, assembling these data with the previous breeding information (Appendix S7).
Aiming to improve the genetic contribution of each founder, and to avoid inbreeding, six preferential couples are suggested (Table 4).  (b) the uniqueness of ecological and reproductive traits, such as its diet based on birds as adults and the hemipenis in females (Hoge et al., 1953;Marques et al., 2012;Martins et al., 2002), its unique venom which toxicity changes ontogenetically and is higher upon birds than mammals (Zelanis et al., 2008), and the possibility of using this species as a model for evolutionary and ecological studies (Duarte et al., 1995).

| Population genetics of B. insularis
Heterologous amplification is a tool widely used in studies of conser-    B. mattogrossensis, and B. pauloensis;Machado, 2015) and other continental species (Duan et al., 2017;King, 2009). Low genetic diversity was also recovered in other critically endangered snakes (Jaeger et al., 2016;King, 2009;Meister et al., 2012;Wang et al., 2015); heterozygosity values found herein were similar to the ones compiled by King (2009) (Piry et al., 1999). However, reduction in the number of alleles is faster than the decrease in heterozygosity rates, so that in some cases observed heterozygosity rates may be higher than the expected ones, which are based on the allele frequencies (King, 2009;Luikart & Cornuet, 1998). It might be considered that populational growth after a bottleneck event may increase heterozygosity rates (Allendorf et al., 2013).

TA B L E 3 Parentage analyses of offspring with uncertain paternity and unknown mother; p-values obtained through the ML Related analyses
Demographic history records of B. insularis indicate that the in situ population declined in the past 20 years due to bushfires and anthropogenic actions. The last populational census recovered a high density and stability in the in situ population, with records of growth rates (Guimarães et al., 2014). Therefore, even though our results did not recover evidences of a recent bottleneck, we may hypothesize that the heterozygosity rates obtained herein might be explained by the existence of bottleneck effect during the last two decades (or the founder effect in the ex situ population) associated with demographic growth in both populations during the last years.
We recovered a significant negative value of inbreeding mating system coefficient for both ex situ and in situ populations. These results contradict the tendency for positive inbreeding values observed in snakes (see King, 2009). Though the analyses indicate that outbreeding is the mating system in both populations, the in situ population showed a significant positive value of kinship inbreeding coefficient, suggesting the occurrence of endogamy, which is the mating of related individuals (Allendorf et al., 2013). Endogamy was also reported in another endemic insular snake, Gloydius shedaoensis (Wang et al., 2015). In their study, the authors discuss the possibility that this species evolved dispersal strategies to avoid endogamy.
Likewise, we hypothesize that the outbreeding mating system could also have evolved in B. insularis as a response to endogamy.
Furthermore, the existence of inbreeding may reduce the genetic fitness in the population, leading, in some cases, to a decline in fecundity and survival rates, as well as sexual abnormalities.
Reproductive studies in B. insularis have been reported hemipenis in different developmental stages in females (Hoge et al., 1953;Kasperoviczus, 2009), small fecundity of this species when compared to its mainland relative B. jararaca (Marques et al., 2013), and high level of mutations in males sperm which lead to a reduction in the number of the viable ones (Silva et al., 2015). Therefore, even though we used neutral molecular markers, the data generated herein may be considered an initial step to associate molecular and reproductive data on B. insularis.
The genetic divergence between populations could be explained

| Conservation implications
Although   Frankham, 2010;Witzenberger & Hochkirch, 2011). Additionally, the so-called sampling address of the founder individuals must be selected in order to minimize inbreeding and outbreeding effects in the ex situ population (Újvari et al., 2002), taking into account that it should retain 90%-95% genetic similarity to its wild counterpart (Castellanos-Morales et al., 2016;Frankham et al., 2008). We observed a mere 49% of genetic similarity between the in situ and ex situ populations, explained by the low but significative divergence between them due to distinct allelic frequencies. Therefore, we strongly suggest that new individuals from the island should be incorporated into the ex situ population; besides, a breeding protocol (based on genetic and reproductive information) should also be established, in order to choose the ideal couples which may reduce such biases, maintaining the genetic similarity between the populations, as well as avoiding loss of genetic diversity and inbreeding. We also recommend that these data should be associated with ecology and reproductive information guiding future research on this species, so that the genetic data will be used as a base to maintain a healthy and viable captive population, with highly genetically similar to the in situ one, which is crucial for future reintroductions.
Several reproduction and management strategies are proposed to maintain genetic variability and to avoid inbreeding and captive adaptations in the ex situ population. Fernández and Caballero (2001) pointed out that strategies focusing on the reduction of mean kinship among individuals, dampening of the founder effect, or a combination of both strategies are widely used in ex situ conservation programs. Another important issue that must be accounted for is the lack of parentage information for some individuals, which may bias management strategies due to underestimation of kinship and inbreeding within populations (Jiménez-Mena et al., 2015). Our results allowed us to refute a potential relationship between a putative father and three descendants and identify four maternal founders related to five descendants, whose mothers were undetermined.
The average kinship value suggests that the ex situ population is composed of half-siblings and representatives with small parentage degree levels (e.g., cousins, uncles, nephews, grandparents), which can be explained by the high genetic contribution of few founders (see Appendix S6). Based on the kinship values and the strategies cited above, we suggest six preferential crosses to the ex situ population (Table 4).
Importantly, a significant inbreeding coefficient was found in the in situ population. Studies have shown that inbreeding is correlated with the fixation of deleterious genes and, when reaching inbreeding depression, results in extinction (Allendorf et al., 2013;Brook et al., 2002;Frankham et al., 2008). Therefore, continuous monitoring of genetic parameters and delimitation of strategies to reduce inbreeding is critical to guarantee the successful conservation of the in situ population, as well as to access genes which may be related to the fitness of this species.
Due to the importance of this species, the Laboratório de Ecologia e Evolução of the Instituto Butantan is developing an interdisciplinary conservation program aiming to implement successful It is also worth mentioning that all the information generated herein should provide the base of the development and implementation of the conservation center of the Instituto Butantan' Zoo Park.

| CON CLUS IONS
Bothrops insularis showed levels of genetic diversity similar to those Lastly, we are grateful to anonymous reviewers for their valuable comments and suggestions on the manuscript.

CO N FLI C T O F I NTE R E S T
The authors have knowledge about the content within this manuscript and declare the absence of conflict interest.