Subspecies hybridization as a potential conservation tool in species reintroductions

Abstract Reintroductions are a powerful tool for the recovery of endangered species. However, their long‐term success is strongly influenced by the genetic diversity of the reintroduced population. The chances of population persistence can be improved by enhancing the population's adaptive ability through the mixing of individuals from different sources. However, where source populations are too diverse the reintroduced population could also suffer from outbreeding depression or unsuccessful admixture due to behavioural or genetic barriers. For the reintroduction of Asiatic wild ass Equus hemionus ssp. in Israel, a breeding core was created from individuals of two different subspecies (E. h. onager & E. h. kulan). Today the population comprises approximately 300 individuals and displays no signs of outbreeding depression. The aim of this study was a population genomic evaluation of this conservation reintroduction protocol. We used maximum likelihood methods and genetic clustering analyses to investigate subspecies admixture and test for spatial autocorrelation based on subspecies ancestry. Further, we analysed heterozygosity and effective population sizes in the breeding core prior to release and the current wild population. We discovered high levels of subspecies admixture in the breeding core and wild population, consistent with a significant heterozygote excess in the breeding core. Furthermore, we found no signs of spatial autocorrelation associated with subspecies ancestry in the wild population. Inbreeding and variance effective population size estimates were low. Our results indicate no genetic or behavioural barriers to admixture between the subspecies and suggest that their hybridization has led to greater genetic diversity in the reintroduced population. The study provides rare empirical evidence of the successful application of subspecies hybridization in a reintroduction. It supports use of intraspecific hybridization as a tool to increase genetic diversity in conservation translocations.


| INTRODUC TI ON
Reintroductions are an important and powerful tool in the conservation and recovery of endangered species. The long-term goal of a viable self-sustaining population is strongly dependent on the genetic makeup (Seddon & Armstrong, 2016), yet evolutionary genetic considerations are often neglected in applied conservation management (Mace & Purvis, 2008). Small numbers of founders and genetic bottlenecks experienced during the establishment phase cause an increased risk of inbreeding and reduction in genetic diversity due to random genetic drift (Frankham et al., 2002).
Consequently, reintroduced populations often display low levels of genetic diversity and reduced evolutionary potential which can impact the populations' chances of long-term persistence (Jamieson, 2010). The capture of sufficient genetic diversity is therefore critical for reintroduction success and should be prioritized when selecting founders (Seddon & Armstrong, 2016). A major practical difficulty is that for many endangered species, potential source populations are small and the number of individuals available for translocation is often limited. In these cases, it has been suggested that genetic diversity could be increased by mixing individuals from different source populations (Jahner et al., 2019;McLennan et al., 2020;Neuwald & Templeton, 2013). Indeed, Biebach and Keller (2012) compared 40 reintroduced populations of Alpine ibex (Capra ibex) and reported that the degree of admixture in founder individuals had a greater impact on heterozygosity in the established population than the number of founders.
The mixing of source populations prior to a reintroduction is not without risks, with the potential for outbreeding depression having been highlighted (Edmands, 2007;Huff et al., 2011). If founders stem from different geographical or ecological regions, individuals may have developed local adaptations. Population admixture is expected to break up these coadaptations, resulting in reduced fitness of the hybrid descendants, especially in later generations (Tallmon et al., 2004;Templeton et al., 1986). Additionally, genetic incompatibility or behavioural differences between the source populations may prevent successful interbreeding of founder individuals (Gottsberger & Mayer, 2019). Complete or partial admixture barriers have been reported between different subspecies and populations of the same species (Soland-Reckeweg et al., 2009). For example, female brown boobies (Sula leucogaster) actively selected against males of a different color morph, thereby preventing hybridization of different genetic clusters within the same species (López-Rull et al., 2016). A barrier to admixture could impact the success of a reintroduction by effectively creating two cryptic populations of smaller size, leading to increased extinction risk.
The application of intraspecific hybridization as a conservation management tool, while much debated, remains highly controversial (Allendorf et al., 2001). This is also due to a general lack of empirical studies on immediate and long-term effects on the translocated populations. The reintroduction of the Asiatic wild ass in Israel provides a rare study system to investigate the impact of this potentially powerful yet controversial strategy on reintroduction outcome. In 1968, a captive breeding core was established with 11 individuals from two different subspecies, the Iranian onager (Equus hemionus onager) and the Turkmen kulan (E. h. kulan) (Saltz et al., 2000). Breeding was unmanaged, and individuals of the different subspecies were allowed to interbreed (Saltz & Rubenstein, 1995). Between 1982 and1987, 28 individuals (14 females, 14 males) and between 1992 and 1993 an additional 10 individuals (7F, 3M) were released into the Negev desert. The reintroduced population has since rapidly increased in size and expanded its range across a large geographical area in Southern Israel (Gueta et al., 2014). No signs of outbreeding depression in the population have been observed. In fact, the reproductive rate of the population is higher than in other reintroduced Asiatic wild ass populations, the reproductive success of females is high, and the reintroduction is currently considered a success Saltz & Rubenstein, 1995). Nonetheless, a major question is whether this success followed admixture of the two subspecies (potentially increasing genetic diversity) or whether genetic or behavioural barriers led to distinct cryptic populations. Differences in habitat selection between parental subspecies or between hybrids and parental subspecies could prevent complete admixture. Indeed, a previous analysis of subspecies-specific mitochondrial haplotypes indicated spatial structuring of the population (Gueta et al., 2014).
The aim of this study was to investigate the long-term genetic consequences of sourcing founders from two different subspecies in a conservation reintroduction. We first explored levels of admixture in the historic captive breeding core and in the current wild population. Then, we investigated population genomic parameters including effective populations size estimates and heterozygosity levels, to evaluate the genetic status of the populations before and after release.

| DNA sample collection
In 1991, three generations after the establishment of the captive population, whole blood samples were collected from 25 individuals in the breeding core (hereafter founder population; Gueta et al., 2014). Between 2011 and 2017, a total of 33 blood and tissue samples of the reintroduced population (hereafter wild population) were collected opportunistically during veterinary treatments, fitting of radio collars and from animals killed in traffic accidents from across the population's range ( Figure 2; Table S1). In addition, we obtained tissue and whole blood samples of the E. h. onager (N = 6) and E. h. kulan (N = 15) subspecies from captive populations in European zoos, collected opportunistically during veterinary treatments or from dead individuals (Table S1). Whole blood samples were stored in EDTA tubes (BD Vacutainer K2EDTA 18.0 mg; Thermo Fisher Scientific; Vacuette K3EDTA 3 mg, Greiner Bio-One), and tissue samples were stored either untreated in paper bags or in screw-cap tubes in 70% ethanol. All samples were stored frozen (at −20°C or −80°C).

| DNA extraction, ddRADseq library preparation and sequencing
We used commercial silica spin column-based extraction kits (Qiagen DNeasy Blood and Tissue Kit, Thermo Fisher Scientific GeneJET Genomic DNA Purification Kit) to purify DNA, following the manufacturers protocol. Samples were digested with the high-fidelity versions of the restriction enzymes EcoRI (R3101S; New England Biolabs) and Sbfl1 (R3642L; NEB). The ddRADseq library preparation followed a protocol adapted from Peterson et al. (2012). Libraries were pooled in equimolar amounts (100 nM) and sequenced on one lane of an Illumina HiSeq4000 flowcell. Sequencing produced over 400 million raw paired end reads with a mean read length of 300 bp.
We assessed the quality of raw reads with the FastQC tool (Andrews, 2010) and a mean Phred+33 quality score >30 was recorded for all bases. Reads were de-multiplexed, and barcodes and Illumina adapters were trimmed using the process_radtags script in the Stacks 2.0 pipeline (Catchen et al., 2013). We used seven replicate samples and a parameter optimization approach adapted from Paris et al. (2017) to assemble loci de novo (optimal parameters: -m3 -N0 -M4 -n4). Subsequent alignment of de novo assembled loci to a reference genome resulted in significant reduction of assembled loci (Table S2).
Since mean individual coverage was exceptionally high (114×) and SNP error rates calculated from seven replicate pairs of samples were low for de novo assembly (mean ± SD = 0.011 ± 0.003; Table   S2), we decided to continue with the de novo assembled loci. Called SNPs were filtered in VCFtools (minimum mean individual coverage ≥35×, minimum minor allele count ≥3, SNPs present in minimum of 80% of individuals, Hardy-Weinberg equilibrium outliers for α = 0.05; Danecek et al., 2011). Eight samples were removed due to low sequencing quality. The final data set contained 69 unique individuals and 4231 SNPs.

| Admixture of the subspecies
To investigate levels of admixture in the Israeli populations, we calculated hybrid indices for each individual. First, we analysed the zoo samples (sample sizes: kulan n = 9, onagers n = 5) to identify diagnostic SNPs, which are fixed for alternative alleles in the two subspecies and calculated hybrid indices as the proportion of alleles inherited from either subspecies (M. Kardos, pers. comm.). Out of the total of 4231 SNP positions, two were not detected in any of the onagers and were removed from the analysis. Of the remaining 4229 SNPs, only 13 (0.3%) were found to be fixed for different alleles in the two subspecies populations.
We calculated individual hybrid indices both from diagnostic SNPs and using maximum likelihood (ML) methods implemented in the R package introgress (Gompert & Buerkle, 2010) on the entire data set. introgress provides estimates of parental allele frequencies and calculates hybrid indices, accounting for uncertainty in allele ancestry when parental populations share alleles. The zoo onager and kulan samples were specified as parental populations with allele ancestry information being used to calculate individual hybrid indices (as the proportion of kulan ancestry) for other samples. We further investigated admixture levels using pairwise fixation indices (F ST ) and calculated for all pairs of populations in the hierfstat R package (Goudet, 2005). We tested Weir and Cockerham's pairwise F ST values for significant deviation from 0 using the function boot.ppfst and 10,000 bootstrap permutations. Additionally, we investigated individual admixture levels using the Bayesian cluster algorithm implemented in the program Structure (Pritchard et al., 2000). To investigate the admixture between two ancestral populations, we ran Structure with the admixture model and correlated allele frequencies, and we predefined the number of genetic clusters (K) as 2. We set the alternative ancestry prior (POPALPHAS = 1) which allows for distinct α-values for each ancestral population, as recommended when sampling is unbalanced (Wang, 2017). Runs were performed with 1 × 10 6 iterations of the Markov Chain Monte Carlo (MCMC) chain preceded by 1 × 10 5 burn-in iterations and 10 repetitions.
We further investigated the relationship between the sampled populations using multivariate methods in the adegenet R package (Jombart, 2008). We first used a principal component analysis (PCA) to explore the data set without a priori groupings. We then performed a discriminant function analysis of the first 15 principal components (referred to as DAPC), which accounted for 52% of the total variance. The DAPC fits orthogonal discriminant functions that maximize between group relative to within-group variation for a priori defined groups (Jombart et al., 2010). We defined four groups consistent with the sampled populations.
To test for a potential barrier to admixture caused by different habitat preferences between the subspecies, we investigated spatial autocorrelation between individuals of the wild population based on hybrid indices in ArcGIS 10.0 (ESRI, 2011). We set a fixed distance threshold of 3400 m between sample locations with row standardization to ensure that each individual had at least one neighbour.

| Heterozygosity, inbreeding and effective population sizes
To avoid any bias due to obvious linkage disequilibrium between SNPs, we used a reduced data set with only 1 SNP per locus (total N = 1738) for this analysis. We calculated expected (H e ) and observed heterozygosity (H o ) and the mean inbreeding coefficient across loci for each population using the pegas R package (Paradis et al., 2010). We expressed individual heterozygosity as the proportion of heterozygote markers for each individual. We estimated variance (N ev ) and inbreeding (N ef ) effective population sizes using the program NeEstimator V2 (Do et al., 2014). N ev refers to the size of an ideal population which displays the same sampling variance in allele frequencies as the focal population (Nei & Tajima, 1981). We estimated N ev using the temporal method (Waples, 1989) with the founder and the wild population used as two samples of the same population at generation 0 and generation 3 (based on generation time of 7.5 years, Ransom et al., 2016). Standardized variance in allele frequencies were computed using the method described by Pollak (1983). N ef describes the size of an ideal population with the same probability of alleles being identical by decent as the population in question. We used the linkage disequilibrium method to estimate N ef for the founder and the wild population (Waples & Do, 2008). For both estimates of N e , we set the lowest allele frequency to 0.02 to avoid bias in estimates and we calculated jackknifed 95% confidence intervals (CIs), as recommend for larger numbers of markers (Do et al., 2014).  Figure S1). Nonetheless, these differences had only a minor impact on the biological interpretation of the results, as both methods indicated high admixture levels. Since the identified diagnostic SNPs were few with unknown locations in the genome, they could be clumped and poorly represent genome-wide admixture levels.

| Admixture of the subspecies
Therefore, we used hybrid indices estimated from maximum likelihood methods in subsequent analyses. The founder and the wild population displayed high levels of individual admixture ( Figure S2).
These values indicated a very even subspecies ancestry in the founder population which shifted after release to a slightly higher proportion of onager ancestry in the wild population. This was consistent with pairwise F ST values, which were all significant at p < 0.05.
The founder and wild population both displayed a stronger genetic differentiation from the kulans than from the onagers (Table 1A).
In the Bayesian clustering approach, the two inferred ancestral populations did not coincide with the subspecies samples. While kulans were clearly differentiated from the other populations, onagers displayed greater levels of admixture, comparable to the Israeli populations ( Figure 1a). This is consistent with the results of the multivariate analyses. The PCA clearly separated the kulans from the other samples along PC1 (11.03% of total variance), while there was some overlap between the onagers, founder and wild populations along PC2 (6.26%; Figure 1b

| Heterozygosity, inbreeding and effective population sizes
The observed heterozygosity was significantly higher than expected in all population (Table 1B) Table 1B).
Effective population size estimates for the wild and the founder population were small. For the wild population, the temporal method estimated a variance effective population size (95% CIs) of N ev = 17.1   ulation, yet with a slight bias towards increased onager ancestry.

| D ISCUSS I ON
No spatial autocorrelation arising from subspecies ancestry was detected in the wild population. These results indicated that there were no genetic or behavioural barriers to subspecies hybridization and that no cryptic subpopulations have formed due to a lack of admixture. This may be due to weak genetic differentiation between the subspecies; we found that only 0.03% of the analysed SNPs were fixed for opposite alleles in the onager and kulan samples, which is consistent with previous findings demonstrating only low genetic divergence between onagers and kulans (Bennett et al., 2017). Species such as the Asiatic wild ass, which show low levels of divergence between subpopulations or subspecies may therefore present good candidates for hybridization as part of genetic management.
Further supporting the potential benefits of hybridization as a conservation tool, we discovered that the hybrid founder population showed significantly higher observed heterozygosity than the zoo onager and kulan populations. There was a 10.85% reduction in observed heterozygosity from the founder to the wild population, although this was still above levels found in the zoo onager and kulan populations and was significantly greater than the expected heterozygosity for the wild population. A loss of heterozygosity is commonly observed in reintroduced populations and likely the result of a genetic bottleneck due to a subset of just 38 individuals from the founder population being released into the wild (Broders et al., 1999;Grossen et al., 2018;Mock et al., 2004). This could be further intensified by strong genetic drift due to extreme polygyny in the species Renan et al., 2018). Similarly, reintroduced populations often display very small effective population sizes as found in the Israeli wild ass population (Manlick et al., 2018;Murphy et al., 2018). However, if the effective population sizes remain small over an extended period, this can seriously threaten the populations' long-term persistence (Franklin, 1980). The zoo populations of kulan and onager also showed significant heterozygote excess. This could be due to managed breeding by zoos as both subspecies are part of the European Endangered Species Programs of the European Association of Zoos and Aquaria, which aims to maximize retained genetic diversity in the captive populations (EAZA, 2019). Heterozygote excess could also be enhanced by population substructure, which was previously reported for the captive onager population, whereby heterozygosity could be enhanced if individuals from different clusters were interbred (Nielsen et al., 2007). Stochastic effects causing random changes in allele frequencies between the sexes can also result in an observed heterozygote excess, commonly observed in small populations (Templeton, 2018).
We found that the Bayesian clustering algorithm and the multi- Overall, this study provides valuable insights that are useful for the management of remaining onager and kulan populations. Both subspecies are currently classified as endangered by the IUCN and a major conservation concern is the extreme fragmentation and small sizes of remaining wild populations (Bennett et al., 2017;Kaczensky et al., 2018). Furthermore, on the basis of numerous phylogenetic analyses (Geigl & Grange, 2012;Orlando et al., 2009;Vilstrup et al., 2013), some authors advocate a revision of the separate subspecies status of E. h. onager and E. h. kulan and joint management of the remaining populations (Bennett et al., 2017;Oakenfull et al., 2000).
The results presented here demonstrate that genetic or behavioural barriers are unlikely to compromise a mixed stock management approach in in situ and ex situ conservation. and adaptive potential (Biebach & Keller, 2012;Mueller et al., 2020;Tordoff & Redig, 2001 (Frankham et al., 2011). The resulting reduced reproductive fitness and loss of local adaptations is often cited as the reason to avoid mixing as it could jeopardize successful population establishment (Edmands, 2007;Templeton et al., 1986). Additionally, there may be concerns of compromising the genetic uniqueness of rare populations. While the risk of outbreeding depression is real and deserves careful consideration, it has likely been overstated in many cases (Frankham, 2015), whereas inbreeding depression is a more imminent threat to many small and isolated populations (Frankham et al., 2011). Moreover, genetic differentiation between populations due to recent isolation and genetic drift must not be mistaken for local adaptations (Templeton, 1996 where interbreeding and individual fitness of hybrid offspring are monitored in a captive breeding facility for at least two generations, as fitness effects may be obscured in the F1 generation (Edmands, 2007). If no fitness reductions are recorded, the final (third) stage involves release of admixed individuals into the wild.
While there is increasing evidence for the potential benefits of mixed source reintroductions, long-term genetic data of successfully admixed reintroduced populations are scarce. Yet such case studies are important for evaluating this strategy as a conservation tool and advancing general conservation translocation guidelines. Our study therefore represents an advance through provision of valuable empirical data on a successful and complete admixture between individuals from two different subspecies.

ACK N OWLED G EM ENTS
We would like to thank the Israeli Nature and Parks authorities, University. This is publication 1091 of the Mitrani Department of Desert Ecology.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in NCBI at https://www.ncbi.nlm.nih.gov/, SRA accession: PRJNA650264 (Zecherle et al., 2020).