Lessons for conservation management: Monitoring temporal changes in genetic diversity of Cape mountain zebra (Equus zebra zebra)

The Cape mountain zebra (Equus zebra zebra) is a subspecies of mountain zebra endemic to South Africa. The Cape mountain zebra experienced near extinction in the early 1900’s and their numbers have since recovered to more than 4,800 individuals. However, there are still threats to their long-term persistence. A previous study reported that Cape mountain zebra had low genetic diversity in three relict populations and that urgent conservation management actions were needed to mitigate the risk of further loss. As these suggestions went largely unheeded, we undertook the present study, fifteen years later to determine the impact of management on genetic diversity in three key populations. Our results show a substantial loss of heterozygosity across the Cape mountain zebra populations studied. The most severe losses occurred at De Hoop Nature Reserve where expected heterozygosity reduced by 22.85% from 0.385 to 0.297. This is alarming, as the De Hoop Nature Reserve was previously identified as the most genetically diverse population owing to its founders originating from two of the three remaining relict stocks. Furthermore, we observed a complete loss of multiple private alleles from all populations, and a related reduction in genetic structure across the subspecies. These losses could lead to inbreeding depression and reduce the evolutionary potential of the Cape mountain zebra. We recommend immediate implementation of evidence-based genetic management and monitoring to prevent further losses, which could jeopardise the long term survival of Cape mountain zebra, especially in the face of habitat and climate change and emerging diseases.


Introduction
The Cape mountain zebra (CMZ, Equus zebra zebra) provides a useful case study to help understand and advance the usefulness of genetic tools and monitoring for biodiversity conservation. The Cape mountain zebra is a subspecies of mountain zebra that is endemic to South Africa. The subspecies is listed as Vulnerable on the International Union for Conservation of Nature (IUCN) Red List [1] and on Appendix II by the Convention on International Trade in Endangered Species (CITES; [2,3]). Historically, Cape mountain zebra had a widespread distribution in the mountainous Fynbos, Karoo, and grassland regions of the Western complex interactions between the aetiologic agent, the environment and the host genome. [13,14,15]. In addition, the two relict CMZ populations (GNR and KNR) further exhibited inflated genetic differentiation due to genetic drift and inbreeding effects resulting from lack of dispersal [7].
To manage and monitor the evolutionary potential of the CMZ metapopulation, a Biodiversity Management Plan (BMP) for this species in South Africa was developed by various stakeholders. The purpose of the BMP is to ensure the long term survival of CMZ in nature, using a strategy underpinned by specific goals and objectives aimed at addressing the threats faced by this subspecies, such as population fragmentation, disease, inbreeding, hybridization with plains zebra and habitat loss [6,9]. It was suggested by several authors [6,7,8,15,16,17] that mixing of aboriginal populations is required to further reduce genetic diversity losses, especially considering that only populations descending from MZNP stock are experiencing high population growth whereas KNR and GNR are not.
However, it was recommended that introductions into either KNR and/or GNR populations be avoided due to the observed population structure, and mixing should only be considered at alternative locations which would be monitored. These plans were hindered, however, by recent ecological studies that suggested CMZ numbers in the GNR were too low to risk removal of individuals to seed mixed populations [18,19], despite the high diversity reported for the mixed DHNR population. Currently, an estimated 40 mountain zebra are being removed from the MZNP for the purpose of re-establishing a breeding herd within the historical range as well as stocking private reserves with animals, both within and outside the natural distribution range. To date, Cape mountain zebra occur in more than 75 localities, including over 30 national parks. Population sizes are estimated to vary between 4 and 1,191 individuals, with the largest population found in the MZNP. The average annual population increase for the subspecies over the period 2009-2015 was 11% [8].
Although the CMZ metapopulation continues to grow, it is yet to be determined whether management strategies have affected genetic diversity over time. Thus, this study aims to investigate temporal changes in genetic diversity in three key CMZ subpopulations by comparing present day (2015-2016) levels with those of samples collected from 1999 to 2001 [7]. This represents a time span of approximately 15 years or approximately 1.5 generations [20,6]. Genetic diversity within populations can only be expected to increase through gene flow between relict stocks (e.g., as in DHNR) or new mutations, with the latter considered inconsequential in the timescale involved. Natural selection is unlikely to increase diversity except at individual loci under balancing selection. Given that our survey is only across a single generation, and there has been no gene flow between relict stocks during this time, we do not expect diversity to have increased over the study period.
However, diversity may be maintained in larger populations experiencing rapid growth. Therefore, we expect that within one generation, populations with larger size and high growth rates (e.g., MZNP) should maintain diversity, whereas populations with lower growth (e.g., DHNR), and smaller size (KNR), would be expected to show a loss in genetic diversity. Furthermore, changes in genetic diversity such as loss of heterozygosity and rare alleles may also have consequences for how populations are structured relative to each other. The loss of shared alleles from populations would be expected to inflate genetic structure (differentiation), as seen in Moodley and Harley [7], whereas a loss of private alleles would reduce genetic differentiation, making populations appear more similar. The results of this study will be used to inform management strategies employed by the CMZ BMP, by providing additional data on population diversity and differentiation.

Ethical approval and sample collection
Ethical approval was obtained from the Research Ethics and Scientific Committee (RESC) of the National Zoological Garden, South African National Biodiversity Institute (NZG SANBI, NZG/RES/P17/19), as well as the Animal Ethics Committee of University of the Free State (UFS-AED2017/0011). Permission was also obtained from the Department Agriculture Forestry and Fisheries of South Africa under Section 20 of the Animal Disease Act 1984 with (Ref: 12/11/1/1/8). The CMZ were chemically immobilised by helicopter. The dosages of sedation and reversal drugs as well as administration were carried out by a qualified, licensed veterinarian, registered with the South African Veterinary Council. Whole blood samples (5 ml) were collected by a qualified veterinarian from the MZNP (n = 75) and DHNR (n = 27) in 2016 and four samples were obtained from KNR in 2015. In addition, this study includes samples from MZNP (n = 12), DHNR (n = 15) and KNR (n = 9) collected between 1999 and 2001 [7]. Thus, a total of 142 samples were analysed. Samples were stored at -20˚C in the Biobank of the NZG, SANBI until used.

Molecular methods
We extracted DNA using the Quick-DNA Universal kit (Zymo Research) following the manufacturer's protocol for blood. We selected 14 cross-species microsatellite markers (AHT21, UCDEQ505, HTG3, HTG7, HTG9, HTG11, HTG14, HTG15, LEX20, LEX52, TKY273, VHL47, HMB1 and COR014) used in the study conducted by Moodley and Harley, (2005). Polymerase Chain Reaction (PCR) amplification was conducted in a 12.5 μl reaction volume consisting of Ampliqon Taq DNA Polymerase Master Mix RED, forward and reverse primers (0.5 μM each), and 50 ng genomic DNA template. The conditions for PCR amplification were as follows: 5 min at 95˚C denaturation, 30 cycles for 30 sec at 95˚C, 30 sec at 55-60˚C (depending on the marker amplified, S1 Table), and 30 sec at 72˚C, followed by extension at 72˚C for 40 min in a T100 Thermal Cycler (Bio-Rad Laboratories, Inc.). PCR products were run against a Genescan 500 LIZ internal size standard on an ABI 3130 Genetic Analyzer (Applied Biosystems Inc.). Samples were genotyped using GeneMapper v. 4.0 software (Applied Biosystems Inc.).

Genetic variation
MICRO-CHECKER software [21] was used to detect possible genotyping errors, allele dropout, and null alleles. Allelic richness (A r ), was estimated correcting for sample size through rarefaction using HP-RARE v. June-6-2006 [22]. Allele Frequencies, observed Heterozygosity (H o ) and expected heterozygosity (H e ) and number of private alleles per population was calculated using GenAlEx 6.5 [23,24]. To determine the significance of changes (H e or A r ) between the two time periods, a one tailed pairwise T-test (α = 0.05 and α = 0.1) was performed with a null hypothesis that no loss in diversity has occurred. Deviations from expected Hardy-Weinberg (HW) proportions were tested (Markov Chain length of 105 and 100,000 dememorization steps). We also tested for gametic disequilibrium between all pairs of loci using the exact test described by Guo and Thompson [25] in GenAlEx 6.5.

Bottleneck simulations
The programme Bottleneck version 1.2.02 [26] was used to detect evidence of recent population bottlenecks. This programme measures significant differences between the measured expected heterozygosity (H e ; i.e., gene diversity) and the theoretical expected heterozygosity assuming mutation-drift equilibrium (H eq ) given the sample size and observed number of alleles (A). This test identifies populations that have recently undergone a decline in effective population size (N e ), resulting in a heterozygosity excess and deficit of rare alleles [27]. It was suggested that, through testing for bottlenecks using heterozygosity excess, a population size reductions of~50 N e , occurring approximately 25-250 generations ago can be detected [27]. Previously recommended parameters for microsatellite data, using the two-phase mutation model, were used [28,27]. This model accommodates for a small proportion of multiple-step mutations, with most mutations being a single step change in allele length. We used two different mutation models including 95% and 80% single-step mutations (SMM) with a two-phase model variance of 12% for a total of 10,000 iterations [29]. Heterozygosity excess was tested for using the Wilcoxon sign-rank test (Significance at α = 0.05) [30]. Allelic frequency, mode-shift deviation from the L-shaped distribution was examined, which was to corroborate detection of a recent bottleneck event [31].

Effective population size (N e ) estimation
The two-sample temporal method [32,33,34] was applied to estimate the variance effective population size (N eV ). This method assumes temporal changes in allele frequencies are caused solely by genetic drift, based on the Wright Fisher model [35,36]. The standardized variance in allele frequency was calculated using two moment based F-statistic estimators, namely Fs [37] and Fc [34] under the model that assumes animals sampled will contribute to future generations. This was implemented using the program NeEstimator v2.1 [38]. This analysis was performed using population samples from MZNP which consisted of adults and foals. The KNR and DHNR were omitted from this analysis due to the small sample size and the recent admixture.

Genetic structure
Since changes in genetic diversity can also affect genetic structure, we determined the relative structure (differentiation) among our three sampled populations at the time periods 1999-2001 and 2015-2016 using three methods. First, we used the Principal Component Analysis (PCA), which is a multivariate method using K-means clustering, and implemented in the R package Adegenet version 2.1.1 [39]. We then used the model based Bayesian clustering algorithm in STRUCTURE version 2.3.4 [40], which determines the most probable number of populations and assigns individuals to their most likely population of origin.
We ran STRUCTURE with the following models: admixture model with both correlated and independent allele frequencies and no admixture model with correlated and independent allele frequencies. Each of the models were run without prior population information for ten replicates each with K = 1-10, with a run-length of 700,000 Markov Chain Monte Carlo repetitions, following a burn-in period of 200,000 iterations. The ten values for the estimated loglikelihood (ln(Pr(X|K)) were averaged across runs and posterior probabilities were calculated. The K with the greatest increase in posterior probability (ΔK, [41]) was identified as the optimum number of sub-populations using STRUCTURE HARVESTER [42]. The membership coefficient matrices (Q-matrices) of replicate runs for the optimum number of sub-populations was combined using CLUMPP version 1.1.2 [43] with the FullSearch algorithm and G 0 pairwise matrix similarity statistics. The results were visualized using DISTRUCT version 1.1 [44]. Lastly, we used an F st -based hierarchical analysis of molecular variance (AMOVA, [45]) to estimate how genetic diversity was partitioned between and within the MZNP and KNR populations for both time periods (Arlequin 3.5; [46]). We excluded the DHNR population for this particular analysis as it is descended from a mix of MZNP and KNR.

Genetic variation in different populations and time periods
Deviations from HW proportions, following Bonferroni correction [47], were only observed in two loci from two populations namely: COR014 (MZNP, DHNR, 2015-2016) and UCDEQ505 (DHNR, 2015-2016). One locus, COR014, showed evidence of null alleles in the MZNP population at both time periods. The locus, UCDEQ505 showed evidence of null alleles in one population at one time period (2015-2016) according to the MICRO-CHECKER results (S2 Table). Thus, all further analysis was performed with and without the marker COR14. The KNR 2015-2016 population was not tested for null-alleles due to insufficient sample size. Significant gametic disequilibrium was not observed between loci in any population. The markers, HTG15, LEX52 and HTG03 were monomorphic in all CMZ samples for both time periods and were omitted from further analysis. Analysis of microsatellite data identified low to moderate genetic diversity (1999)(2000)(2001): H e = 0.511 and 2015-2016: H e = 0.338) in CMZ populations from both time periods (S3 Table) compared to those reported for plains zebra (Equus quagga; H e ranged from 0.71 to 0.80) [48].

Genetic diversity and effective size (N e ): temporal changes within populations
Analysis per population indicated that the highest heterozygosity and A r was observed for the DHNR population at both temporal periods compared to the MZNP and KNR populations ( Table 1). Heterozygosity within each reserve (KNR, DHNR and MZNP) declined between temporal sampling periods (Table 1). A decline in genetic diversity was observed for DHNR in A r , which decreased from 2.  (Table 1). However, only the decline in H e in DHNR was observed to be statistically significant (p > 0.05). The number of private alleles for KNR, MZNP and DHNR in 1999-2001 was 0.143, 0.357 and 0.071 respectively, whereas the frequency of private alleles in these populations in 2015-2016 was 0.0 ( Fig  1). The temporal variance estimates of N e ranged from 1.7 to 6 when performing F c analysis and varied from 1.7 to 18.2 when performing F s analysis ( Table 2). Analysis where the locus COR014 was removed was similar for F c analysis but differed for F s (1.6 to 31.2, Table 2).

Bottleneck tests
The bottleneck tests were only carried out on the MZNP, for the two temporal periods. Similar results were obtained for both 80% and 95% SMM ( Table 1). The DHNR population was not tested as the heterozygote excess method assumes that no recent admixture has taken place in KNR, and the sample size was too low. No significant heterozygote excess was detected for MZNP (p > 0.10).  (Fig 3A), while the 2015-2016 DHNR populations were clustered mainly with KNR (qi = 0.9048, Fig 3B). Analysis of Molecular Variance (AMOVA) ( Table 3)

Table 1. Genetic diversity at the two time points for the Cape mountain zebra in the three reserves, Mountain zebra National Park (MZNP), De Hoop Nature Reserve (DHNR) and Kammanassie Nature Reserve (KNR). Population genetic diversity over for two periods spanning around 16 years (1999-2001 to 2015-2016) was evaluated based on allelic richness (A r ), observed heterozygosity (H o ) as well as unbiased expected Heterozygosity (H e ) and fixation (F IS ) across 11 polymorphic microsatellite loci.
Detection of recent bottlenecks were tested by identifying whether Mountain Zebra National Park had significant levels of heterozygous excess using the Two Phase Mutation model (Wilcoxon P TPM and P SMM ) and then confirmed by checking allelic mode shift at mutation-drift equilibrium (A DIST ). � and bold text = significant temporal change, p < 0.05. 1 indicates estimates of A r are corrected for sample size through rarefaction in HP-RARE using the smallest number of gene copies per population (MZNP = 14, DHNR = 24, KNR = 4) and 2 indicates were analysis was not applicable, since sample size was too small for KNR and in DHNR, the population is an admixture of MZNP and KNR.

Discussion
In this study, we provide a rare test of the genetic consequences of different conservation management strategies among populations of a large mammal. Management recommendations from the BMP included (1) deliberate mixing of relict populations to maintain and improve genetic diversity (excluding KNR and GNR), (2) re-enforcement of existing populations prioritised over the establishment of new populations, (3) translocation of animals to other protected areas, (3) acquisition of land adjacent to protected areas with CMZ, (4) alteration in fire management in the habitat preferred by CMZ to increase availability of palatable grasses; and (5) formation of conservancies with adjacent landowners. However, our data provides Table 2. Estimated effective population sizes for the Mountain zebra National Park (MZNP) using the temporal method. The two analyses used were the F c method (Nei and Tajima, 1981) and the F s method (Jorde and Ryman, 2007) using the Program NeEstmator ver. 2.1 (Do et al., 2014). P crit = the criterion for excluding rare alleles if the frequency of rare alleles and less than the P crit value they are excluded, GI = the Generation interval, N e = the estimated effective population size, Min and Max = confidence interval values and CV = the Coefficient of variation.

MZNP Analysis with marker COR14 Analysis without marker COR14
Fc Fs Fc Fs P crit = 0,05 P crit = 0,02 P crit = 0,05 P crit = 0,02 P crit = 0,05 P crit = 0,02 P crit = 0,05 P crit = 0,02  These results suggest that despite a high metapopulation growth rate, the CMZ has lost a significant proportion of its genetic diversity within a single generation. At the population level, all three reserves lost genetic diversity (A r and H e ), with a statistically significant loss being detected in DHNR. The loss of diversity in DHNR was compounded by decreases in H o and an increase in the inbreeding co-efficient (F IS ). In contrast, non-significant increases in H o and a reduction in F IS were found in the relict populations MZNP and KNR. This slight increase in observed heterozygosity may reflect a sampling effect in KNR, since only four samples were genotyped from this reserve in the 2015-2016 time period. It is more likely due to several loci within these two populations reaching HWE, where H o is not significantly different from H e . Both KNR and MZNP, therefore, did not display significant heterozygous excess.
The N e for the MZNP was very low suggesting that continued isolation of this population will result in further loss of genetic variation through drift. A small effective population size could assist in explaining how genetic diversity was lost even though the population demographic census size increased. The effective size can be far smaller than the census size due to a mating system that produces high variance in male and/or female reproductive success. However, bias in this case may be introduced by small sample size, limited time period between  collection of samples and number of microsatellite markers [33,49]. Bias is unlikely to explain entirely the low N e estimates which are well below N e = 50 and can lead to excessive inbreeding and inbreeding depression [50]. Similar low N e values have been reported in other species including red deer (Cervus elaphus) from Sardinia and Mesola (N e = 2 to 8) as a consequence of bottlenecks and near-extinction [51]. Of additional concern is the complete loss of private alleles from all three sampled population's. This suggests rare and low frequency alleles have likely been lost genome wide. While a loss in such population-specific alleles may not have an effect on overall heterozygosity (Norris et al., 1999), it represents a loss of alleles potentially important for future adaptation and loss the uniqueness of the KNR and MZNP stocks. This means that a proportion of the historical diversity of the Cape mountain zebra, conserved through different conservation strategies and present in 1999-2001 has already been lost in the present generation. The loss of unique alleles in MZNP and KNR could also affect an individual or populations ability to adapt and cope with future environmental change. Equally worrying is that the third stock population, inhabiting the Gamkaberg Nature Reserve (GNR), with a growth rate even lower than that of KNR, could also be similarly affected.
The erosion of private alleles may have affected population genetic structure (differentiation) of these populations. Model-based and model-free algorithms all document a decrease in genetic differentiation between the KNR and MZNP stocks, which are now more similar to each other than they were a generation ago. Furthermore, the allele frequencies for DHNR, which were intermediate between its MZNP and KNR founding stocks in 1999-2001, are now clearly more Kammanassie-like, with 90% of genotyped individuals assigned to the KNR cluster (Fig 3). In addition, these results are supported by a reduction in F ST (0.544 to 0.354). Such a substantial shift or homogenization in allele frequencies within a single generation underscores the erosive effect of random genetic drift, even in populations that are expanding demographically, and threatens to undo much of the population benefits of a strategy to restore gene flow.
Major changes in the number of private alleles could also be brought about by non-random mating, where a handful of males dominate most of the breeding opportunities. This idea is corroborated by empirical data showing that in 2005, the population was already male-biased, with a deficit of females resulting in an excess of non-breeding males with limited reproductive potential. Population growth at DHNR has declined from 6.6% in 1995-1999 to 4.5% in 1999-2005 [52]. Thus, DHNR, a population that previously benefitted from admixture, now requires urgent intervention to mitigate this loss. The lower growth rates of KNR and GNR has been attributed to lower abundance of palatable grass species in those reserves [16,53]. The conservative practice of managing these populations separately to protect their uniqueness has had the opposite effect of loss of alleles that made them unique in the first place. Given the evidence of genetic declines reported in this study, the erosive effects of genetic drift and non-random mating can only be rectified through new introductions. Further suggested management practices to facilitate population growth and promote increased genetic diversity include establishing studbooks for all newly founded mixed-stock populations, the use of fertility-control methods to ensure equity in mating opportunities among males and females [54] and ensuring range quality and hence overall body condition on new reserves identified as suitable for CMZ populations. We therefore advocate changes to the conservation management of these important populations, to try to arrest these worrying population trends, likely to cause further loss of diversity, evolutionary potential and the onset of fitness related problems.
Results and the approach from this study could help design and implement management and conservation strategies in other species with only a few small populations remaining. In addition to census size, genetic monitoring of multiple metrics (e.g., heterozygosity, allelic uniqueness, and effective size) can provide early detection of loss of diversity even when a population is large or growing.
Supporting information S1