Genetic purging in captive endangered ungulates with extremely low effective population sizes

Inbreeding threatens the survival of small populations by producing inbreeding depression, but also exposes recessive deleterious effects in homozygosis allowing for genetic purging. Using inbreeding-purging theory, we analyze early survival in four pedigreed captive breeding programs of endangered ungulates where population growth was prioritized so that most adult females were allowed to contribute offspring according to their fitness. We find evidence that purging can substantially reduce inbreeding depression in Gazella cuvieri (with effective population size Ne = 14) and Nanger dama (Ne = 11). No purging is detected in Ammotragus lervia (Ne = 4), in agreement with the notion that drift overcomes purging under fast inbreeding, nor in G. dorcas (Ne = 39) where, due to the larger population size, purging is slower and detection is expected to require more generations. Thus, although smaller populations are always expected to show smaller fitness (as well as less adaptive potential) than larger ones due to higher homozygosis and deleterious fixation, our results show that a substantial fraction of their inbreeding load and inbreeding depression can be purged when breeding contributions are governed by natural selection. Since management strategies intended to maximize the ratio from the effective to the actual population size tend to reduce purging, the search for a compromise between these strategies and purging could be beneficial in the long term. This could be achieved either by allowing some level of random mating and some role of natural selection in determining breeding contributions, or by undertaking reintroductions into the wild at the earliest opportunity.

Dorcas gazelle (G. dorcas) is one of the smallest African gazelles. Males and females are similar in size (mean body weight 15.9 kg for males and 15.7 for females; López and Abáigar 2013). Both sexes have horns, being bigger in males. Females are fertile at the age of 11-12 months and males at eight months (López and Abáigar 2013). The gestation period is 24-25 weeks (Alados 1982). Female reproductive senescence is around 11 years. There is only one calf per birth (Abáigar 2002). When conditions are harsh, dorcas gazelles live in pairs, but when conditions are more favourable, they occur in family herds with one adult male, several females and young individuals. During the breeding season, adult males tend to be territorial (Beudels et al. 2005). The founder population arrived La Hoya between 1970 and 1975 from Western Sahara, and consisted of 37 males and 37 females (Cano 1988, Abáigar 1993. Dama gazelle (Nanger dama) is a big desert/semidesert gazelle (body mass 40-75 Kg; Barbosa and Espeso 2005) sexually dimorphic in size. Adult males are 25% heavier than adult females. Males and females have horns, these being thicker and longer in males. Females are fertile at the age of 13-15 months (Alados and Escós 1991) and males about 16 (Espeso 2018). The gestation period is about 26 weeks. Female reproductive senescence is around 13-14 years. (Barbosa & Espeso 2005). Herds occur singly (males or females), or in mixed groups of 10 -15, composed of a dominant adult male, several adult females, and young individuals. It is highly nomadic ranging widely in order to obtain sufficient nutrition. These gazelles used to undertake large seasonal migrations, moving north into the Sahara desert during the rainy season, and retreating south into the Sahel during the dry season. Founders of the captive population in Almería has been controversial. Following the criteria of "unknown parents", Cano (1991) considered that 16 individuals (3 males and 13 females) were the founders; coinciding with those animals included in the first published studbook of the species (Barbosa and Espeso 2005). However, Ruiz-López et al. (2009) using molecular information found that the founding population of N. dama mhorr was much smaller. Finally, the history of its captive population at La Hoya was reconstructed based on information taken from Cano (1991) and references therein, Valverde (2004), and information available in the historical files of the Experimental Field Station "La Hoya" (Moreno et al. 2011). The current studbook derived from that reconstruction (Espeso 2018) assumes that five individuals (1 male and 4 females) should be considered as founders of the captive population of this gazelles at La Hoya, as many of the animals arriving at La Hoya from Western Sahara had been living in captivity for years and were actually relatives.
At La Hoya, once young adult males reach the age of 7-9 months, male-male interaction is common, and must be avoided to prevent injury. This social behaviour strongly determines the species management as enclosure must allow healthy animals to display the full range of normal behaviour (Escós 1992). Thus, at La Hoya animals are distributed in three different types of enclosures: a) breeding enclosure, generally occupied by a breeding group composed by one male and a group of 5 to 8 females their young (young males are removed when 7-9 months old to prevent injury by the adult male); b) bachelor enclosure, where young celibate males are kept forming bachelor groups; c) individual male enclosure where solitary males that have been used as breeders are kept (see Moreno & Espeso, 2008

Appendix S2. Individuals filtered
Original pedigree files were downloaded from the studbooks available online at http://www.eeza.csic.es/es/programadecria.aspx, but not all individuals were made available for analysis. Several filtering criteria were introduced, aimed to avoid sources of bias for inbreeding and fitness estimates. Table S1 contains the original number of individuals in the studbooks, as well as the number of individuals removed in each filtering step, and the number of individuals retained in the processed pedigrees. Table S1. Number of individuals during pedigree processing: Number of individuals in the studbook (Ntotal), number of individuals removed due to unknown ancestry (NUNK), number of individuals removed as consequence of accidental death (Naccident), number of individuals at the end of the pedigree, removed for female productivity analysis (Nrecent years), final number of individuals remaining for survival (NWS) and productivity (NWP) analysis. First, individuals with unknown ancestry were removed. This includes individuals recorded in the studbooks with unknown (UNK) as well as "multiple" (MULT) ancestry, meaning that several individuals were candidate parents, but accurate ancestry could not be established. Descendants from these individuals were also removed. This filtering step was aimed to avoid downward biases of inbreeding that could affect analysis of inbreeding depression and purging. The higher number of individuals excluded for G. cuvieri compared to other species (NUNK = 696) is mainly due to the higher degree of transference of this species to international institutions participating in its Ex situ Conservation Programme.

Species
Then, individuals that died younger than 15 days of life as consequence of environmental accidents (e.g. trauma), were also removed, to avoid a downward bias of fitness due circumstances non-connected to genetics. This included individuals with identities 1057, 1381, 1384, 1386 and 1426 in G. cuvieri, and 1654 in G. dorcas. No individuals were removed for this reason for A. lervia and N. dama.
When considering lifetime female productivity as fitness (WP in the main text), females at the end of the pedigree were also removed. The criteria was defined as a function of each species' longevity, so that more long-lived species had more years of data removed from the pedigree. This is a required filter to accurately estimate inbreeding depression and purging for traits such as productivity, since in the more recent years, when individuals are more inbred, yet-alive individuals will have their productivity estimated underscored, and only complete records would be available for individuals that died earlier or died more prematurely. We selected the 90 % percentile of longevity as threshold, and therefore removed the last 16 years of data for A. lervia, 11 years for N. dama, and 9 for the two Gazella species.
It must be noted that we only considered individuals belonging to La Hoya Experimental Field Station (see main text), and individuals from other locations had their fitness values set as NA (non available). Male productivity was always set to NA. As a last filtering step, individuals with NA values of fitness were also removed from the pedigrees, as long as they were not ancestors of individuals with non-NA values of fitness. Since this filtering is conditional to fitness, the number of individuals removed for survival and productivity analysis differs. The motivation of this last filtering step was to remove individuals not required for the remaining analyses, and to keep a shorter pedigree that improve the performance of downstream analyses.

Appendix S3. Complementary demographic parameters
The estimates of Ne reported in the main text were similar to those computed solving for N the classical non-overlapping generations expression = 1 − (1 − 1 2 ) , where F is the average inbreeding for individuals in the last cohort and t is their average number of equivalent complete generations minus one, to account for the absence of self-fertilization: Ne = 3.87 (0.05), 14.10 (0.16), 39.58 (1.40) and 11.18 (0.11) for A. lervia, G. cuvieri, G. dorcas and N. dama, respectively (bootstrap errors parenthetically). This supports the approach proposed by Gutierrez et al. (2008Gutierrez et al. ( , 2009, as well as the use of the classical expression with overlapping generations. The effective number of ancestors (Nea) defined as the minimum number of ancestors, founders or not, required to account for the IBD genetic diversity of the TP. The number of equivalent founder genomes (Ng), which also includes other sources of random loss of alleles such as genetic drift (Caballero and Toro 2000). Then, the existence of recent bottlenecks in the pedigree can be inferred when < 1, and drift can be expected when Ng < Nea. These parameters were estimated following the methods described by Boichard et al (1997). For the estimation of Ng, Monte-Carlo simulation was run for 10 4 iterations. They are reported in the table below. The effective number of equivalent founder genomes (Ng) approximates ½ Nf (Nf is the number of founders) for all species except for G. dorcas, for which it was lower ( Figure S1).  Figure S1. Partition of the number of founders (Nf) for different species (in %). The total height of the bar represents the total Nf. From top to bottom, stacked bars indicate the successive reduction of Nf after subtracting i) the component due to unbalanced founder contributions (Nf-Nef, in light grey), ii) the component due to unbalanced ancestor contributions (Nef-Nae, in grey), the component due to drift (Nae-Ng), and the number of equivalent founder genomes left (Ng).

Appendix S3. Quantitative genetic analysis for survival
Genetic variability for the fitness traits was also evaluated in the full pedigree through the computation of the additive variance (h 2 ), using the R package MCMCglmm (Hadfield 2010), as in Gauzere et al. (2020). In the corresponding animal model, inbreeding and mother identity were included as potential random factors, and the most likely solution was selected using the deviance information criterion (DIC, Spiegelhalter et al. 2002). The dominance matrix was computed using the R package nadiv (Wolak 2012), and was also included in the model. Noninformative priors were chosen for all variables involved. MCMC was run for 3 x 10 5 iterations with a burn-in period of 5 x 10 4 iterations, and storing one in every 250 samples. Output were always check for low autocorrelation (|r| < 0.1), and effective sample sizes close to 10 3 for VA.
After running the animal model, VA and h 2 were computed on the observed data scale using the QGglmm package (Villemereuil et al. 016). The mean-scaled additive variance (IA) was also computed from the previous analyses by dividing the estimated VA over the squared trait mean (Houle 1992, Hansen et al. 2011). An environment with R version 3.6.1 (R Core Team 2019) was used for all analyses.   Table S7. Inbreeding-purging parameters estimated for productivity WP: purging coefficient (d), Pvalue for the estimate of the purging coefficient, rate of inbreeding depression δ (i.e., estimated inbreeding load for Wp ascribed to the individual's genotype), maternal rate of inbreeding depression (i.e., estimate of the maternal inbreeding load δM , ascribed to the maternal component of the trait), sex effect (S, referred to males), period of management effect (POM, referred to the period when veterinary care becomes regular) and initial fitness (W0, numerical estimate of the average fitness of the initial non-inbred population). NA values indicate that the factor was not selected in the model giving the best fit. The last column gives the average fitness for non-inbred individuals with noninbred ancestors and its empirical standard error [W0 F=Fa=0 (±SE)].  Figure S3. Evolution of the expected fitness W predicted under the Inbreeding-Purging model for a panmictic population (including random self-fertilization) using different purging coefficients (d=0.05, 0.2, 0.35 and 0.5) relative to the corresponding fitness expectation in the absence of purging (i.e., using d=0 to compute gi), plotted against generation number for different effective population sizes: red, blue and green lines for Ne = 4, 10 and 40, respectively. Vertical lines correspond to the tm=√(2Ne) values with the same color code (rounded to the closest natural number).