Population genomics analyses of European ibex species show lower diversity and higher inbreeding in reintroduced populations

Abstract Restoration of lost species ranges to their native distribution is key for the survival of endangered species. However, reintroductions often fail and long‐term genetic consequences are poorly understood. Alpine ibex (Capra ibex) are wild goats that recovered from <100 individuals to ~50,000 within a century by population reintroductions. We analyzed the population genomic consequences of the Alpine ibex reintroduction strategy. We genotyped 101,822 genomewide single nucleotide polymorphism loci in 173 Alpine ibex, the closely related Iberian ibex (Capra pyrenaica) and domestic goat (Capra hircus). The source population of all Alpine ibex maintained genetic diversity comparable to Iberian ibex, which experienced less severe bottlenecks. All reintroduced Alpine ibex populations had individually and combined lower levels of genetic diversity than the source population. The reintroduction strategy consisted of primary reintroductions from captive breeding and secondary reintroductions from established populations. This stepwise reintroduction strategy left a strong genomic footprint of population differentiation, which increased with subsequent rounds of reintroductions. Furthermore, analyses of genomewide runs of homozygosity showed recent inbreeding primarily in individuals of reintroduced populations. We showed that despite the rapid census recovery, Alpine ibex carry a persistent genomic signature of their reintroduction history. We discuss how genomic monitoring can serve as an early indicator of inbreeding.

(e.g., habitat, diseases) and genetics seem to influence the outcome (Armstrong & Seddon, 2008). Genetic aspects influence not only the short-term reintroduction success (due to the genetic architecture of the founder individuals), but also the long-term reintroduction success (due to adaptive evolutionary potential and inbreeding). Several studies identified genetic consequences imposed by reintroductions based on microsatellite loci (e.g., Biebach & Keller, 2009;Diefenbach, Hansen, Bohling, & Miller Butterworth, 2015;Jamieson, 2011;Strzała, Kowalczyk, & Łukaszewicz, 2015), but only very few studies surveyed genomewide effects or monitored long-term consequences.
Reintroductions are generally composed of two key stages. The first stage is the transfer of individuals and successful reproduction in a new area; such events were causing founder effects in reintroduced populations with a likely loss of genetic diversity. The maintenance of genetic diversity during the first step is a major challenge during species reintroductions (Jamieson & Lacy, 2012;Tracy, Wallis, Efford, & Jamieson, 2011).
The second stage of a reintroduction is the minimization of longterm demographic effects on populations after the reintroduction event. Major concerns are low effective population size due to slow population growth rates (mainly caused by postrelease mortality and low reproduction rates), lack of gene flow with other populations, and inbreeding. The individuals' ability to move in space and contribute to gene flow among populations is key for the long-term survival of reintroduced populations. Hence, reintroductions in a continuous, connected landscape may be less likely to fail than reintroductions in a strongly subdivided landscape due to anthropogenic or topographical factors.
The different stages of a species reintroduction are predicted to show different signatures of genetic subdivisions. At the early stage of species reintroductions, population differentiation is expected to be weak if reintroduced individuals were drawn at random from a source population. Subsequently, restricted gene flow in a fragmented habitat will lead to genetic drift and reduced genetic diversity. The effect of habitat fragmentation is an important factor for the long-term viability of reintroduced populations. The most immediate consequence of reduced gene flow for a small population is an increase in inbreeding (i.e., mating among related individuals). Inbreeding can have harmful consequences on individual survival and reproductive success (Keller & Waller, 2002). Negative effects of inbreeding on fitness are called inbreeding depression and were widely reported in populations experiencing elevated levels of inbreeding (Keller & Waller, 2002). For instance, in a small island population of song sparrows, several different traits including lifetime reproductive success, survival, and immune response were susceptible to inbreeding (Keller, 1998;Keller, Arcese, Smith, Hochachka, & Stearns, 1994;Reid, Arcese, & Keller, 2003). An increase in inbreeding of 10% led to a decreased male lifetime reproductive success of 25% (Keller, Marr, & Reid, 2006).
Inbreeding depression is an important concern for species conservation (Allendorf, Hohenlohe, & Luikart, 2010;Allendorf, Luikart, & Aitken, 2012;Hoffman et al., 2014;Xue et al., 2015). However, inbreeding depression may only be detectable many generations after an increase in inbreeding. Conservation management can prevent inbreeding depression by translocating individuals from genetically distant populations. Translocations also harbor risks by increasing the likelihood of disease transmission or causing outbreeding depression by breaking up locally adapted allele combinations. Therefore, a careful assessment of the inbreeding risk is an important aspect of successful conservation management actions. Inbreeding is generally estimated using pedigrees (Wright, 1922), which are established from detailed population surveys over multiple generations. In a species with long generation times, pedigree estimations and fitness estimates are difficult to obtain and time-consuming. Hence, such data are only available for few large mammals (e.g., Overall, Byrne, Pilkington, & Pemberton, 2005;Walling et al., 2011). For conservation management practices, pedigree and fitness estimates are impractical to obtain.
Genetic markers can serve to estimate inbreeding directly without the need for pedigrees. Levels of multilocus heterozygosity can be estimated based on neutral genetic markers (e.g., microsatellites).
However, traditional genetic markers such as microsatellites suffer from a number of biases and shortcomings. Microsatellite markers were generally developed to maximize polymorphism in a focal population, which may lead to ascertainment bias in other populations. Recently developed genomic tools based on high-throughput sequencing such as restriction-associated DNA sequencing (RAD-seq) alleviate many of the shortcomings of microsatellite markers. The unbiased genomewide estimates of multilocus heterozygosity generally reflect more accurately levels of inbreeding (Kardos, Taylor, Ellegren, Luikart, & Allendorf, 2016). Recent inbreeding can be estimated by genomewide analyses of runs of homozygosity (ROH, Broman & Weber, 1999;Lencz et al., 2007;Keller, Visscher, & Goddard, 2011;Pemberton et al., 2012).
ROH are estimated from physical marker locations along the genome.
Long ROH are a strong indicator for recent inbreeding as long ROH indicate that the parental individuals were related. While ROH have been extensively studied in humans and domestic animals (e.g., Lencz et al., 2007;McQuillan et al., 2008;Pemberton et al., 2012;Purfield, Berry, McParland, & Bradley, 2012), this measure has rarely been used to estimate inbreeding as a consequence of reintroductions and species conservation. The recent advances in genome sequencing make it feasible to estimate ROH in many nonmodel species (e.g., Xue et al., 2015). Estimating inbreeding based on ROH rather than pedigrees has several advantages (Kardos, Luikart, & Allendorf, 2015;Keller et al., 2011). ROH estimate the fraction of the genome that is identical by descent (i.e., autozygous) and serve to estimate the age of inbreeding events. The power to detect early signs of inbreeding through ROH may help to prevent inbreeding depression before reaching detectable fitness consequences in populations.
Alpine ibex (Capra ibex) are an ideal system to test for the longterm population genetic consequences of population reintroductions.
Alpine ibex are a species of wild goat living at high elevation across the European Alps. The species was overhunted over centuries and went nearly extinct in the early 19th century. An extensive conservation program started in 1906 that reintroduced Alpine ibex to their original distribution range re-establishing the census size from ca. 100 individuals to the current census of ca. 50,000 (Brambilla pers. comm.).
The reintroduction of Alpine ibex to the European Alps was likely the most successful reintroduction of a large mammal. After recovery of the sole remaining population in the Gran Paradiso area of northern Italy to well over 1,000 individuals, 55 male and 45 female ibex were transferred from the Gran Paradiso area to two Swiss zoos to initiate a captive breeding program (Giacometti, 2006). At least 11 individuals died before reproducing, and thus, the maximum captive founder stock consisted of 88 individuals (Stuwe & Nievergelt, 1991). The individuals were first bred in captivity in two Swiss zoos. As all individuals used for the captive breeding originated from the same source population (Gran Paradiso), the zoo-bred population likely started with no genetic subdivisions. Offspring of the zoo-bred population were introduced to different locations in the Swiss Alps. Three reintroduced populations reproduced well (here called "primary" populations) and served as source populations for multiple, secondary reintroductions (here called "secondary" populations). Most of these secondary reintroductions were each initiated with individuals from a single primary population. Secondary populations originating from more than one source population are here called mixed populations (see Table S1 for more details on the reintroductions). Microsatellite-based analyses showed that the stepwise reintroduction strategy left a strong genetic footprint in extant populations (Biebach & Keller, 2009). The repeated establishment of new Alpine ibex populations from previously established populations caused up to four serial bottlenecks (captive breeding, primary, one or more secondary reintroductions). Bottlenecks increase the likelihood of inbreeding due to the reduction in effective population size and increase the risk of inbreeding depression. Indeed, populations with higher levels of inbreeding grew more slowly in harsh environments (Bozzuto pers. comm.). In the Gran Paradiso population, heterozygosity and fitness-related traits were negatively correlated, which suggests that inbreeding can indeed lower fitness (Brambilla, Biebach, Bassano, Bogliani, & von Hardenberg, 2015). Despite effective management and hunting restrictions across the European Alps, some Alpine ibex populations recently declined in size.
The precisely documented reintroduction strategy and the large number of available samples provide a framework to test multiple hypotheses on the genetic trajectory of reintroduced populations.
The Gran Paradiso source population can still be sampled, and hence, genetic diversity in reintroduced populations can be compared to a population that was maintained over centuries. We will use the term autochthonous for such extant, nonreintroduced populations. The primary and secondary population reintroduction history is well documented. Hence, changes in genetic diversity and inbreeding can be assessed as a function of the reintroduction history. Furthermore, the population genetic diversity and structure of a related European ibex species can be used as a comparison. The Iberian ibex (Capra pyrenaica) is the closest relative to Alpine ibex (Pidancier, Jordan, Luikart, & Taberlet, 2006) and was historically found throughout the Iberian Peninsula, from Portugal to southwest France (Grubb, 2005). Similar to Alpine ibex, Iberian ibex suffered from population bottlenecks due to habitat fragmentation and overhunting over the last two centuries (Perez et al., 2002). About ten small populations survived and were legally protected from 1900 to 1960. From 1960s onwards, several Iberian ibex populations were founded by reintroductions (Perez et al., 2002;Refoyo, Olmedo, & Muñoz, 2016). Present-day Iberian ibex populations are mainly found in the mountain ranges of the Sierra Nevada and the Iberian System. Two consecutive and severe bottlenecks (ca. 30 individuals) were reported for the Maestrazgo population between 1940and 1960(Couturier, 1962. The Maestrazgo population recovered to 7,000 individuals (Perez et al., 2002). The Sierra Nevada population suffered from a population bottleneck of ca. 600 individuals in the 1960s and is currently the largest Iberian ibex population with about 16,000 individuals (Couturier, 1962;Perez et al., 2002). Sierra Nevada is also thought to be the genetically most diverse population (Manceau, Crampe, Boursot, & Taberlet, 1999). Despite severe demographic changes over the past century, Iberian ibex were not affected by similar range and census reductions as Alpine ibex. The less severe bottlenecks experienced by Iberian ibex are reflected in higher genetic diversity at the major histocompatibility complex (MHC) (Alasaad et al., 2012;Amills et al., 2004).
We aimed to identify the population genomic consequences of the Alpine ibex reintroduction strategy using Iberian ibex and domestic goat (Capra hircus) as a reference. We addressed the following questions: (i) What was the loss of genomewide diversity due to the reintroduction of Alpine ibex from a single source population? How do the levels of genetic diversity compare with the major Iberian ibex populations? (ii) Did the stepwise reintroduction strategy of primary and secondary populations lead to a stepwise increase in genetic subdivision with subsequent rounds of reintroductions? (iii) Did the reintroductions lead to genomewide signatures of recent inbreeding? Finally, we discuss how genomewide data can be used to assess the long-term viability of reintroduced species.

| Sampling
We obtained 190 samples from 10 Alpine ibex populations, two Iberian ibex (C. pyrenaica hispanica) populations, and four Swiss domestic goat breeds. Nine Alpine ibex populations were reintroduced to the Swiss Alps over the past century, and one was the population from the Gran Paradiso National Park in Italy. The sample size per population was 15 for the reintroduced Alpine ibex and Iberian ibex populations, 16 for the source population Gran Paradiso, and two to three per domestic goat breed (Tables 1 and S2).

| RAD library preparation and SNP calling
Genomic DNA was extracted from tissue or blood samples using a BioSprint 96 kit (QIAGEN). The RAD libraries were prepared as described in Grossen, Keller, Biebach, and Croll (2014). In short, 1.35-1.5 μg of DNA was digested using the restriction enzyme SbfI (New England Biolabs). Adapters containing unique barcodes were ligated to the digested DNA. Ligated DNA was sheared using a COVARIS ultrasonicator and size-selected at 300-700 bp on a CALIPER LabChip XT. After ligation of P2 adapters, a PCR was performed to amplify the library. Quality-checked libraries were sequenced on an Illumina HiSeq 2000 platform to generate 100-bp paired-end reads.
We obtained a total of 468 million high-quality reads and a median of 2.3, 2.1, and 1.6 mio read pairs per Alpine ibex, Iberian ibex, and domestic goat individual, respectively (Fig. S1).
Raw sequencing reads were demultiplexed using the FASTXtoolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html) according to the P1 adapter barcode. Trimmomatic (Bolger, Lohse, & Usadel, 2014) was used for quality and Illumina adapter trimming, as well as the removal of the RAD-seq barcodes. Reads were aligned to the reference genome of the domestic goat using Bowtie2 version 2.2.4 with "end-to-end" and "very-sensitive" settings (Langmead & Salzberg, 2012 McKenna et al., 2010). SNPs were hardfiltered using VariantFiltration and SelectVariants. SNPs were removed if any of the following criteria were met: QD < 2.0, MQ < 30.0, −12.5 > MQRankSum > 12.5, FS > 60.0, ReadPosRankSum < −8.000, QUAL < 30.0, AN < 40. Vcftools (Danecek et al., 2011) was used for all further filtering steps. Genotypes with genotype quality (GQ) <20 were set to missing. A total of 16 individuals were removed from the dataset due to low read number (less than 900,000) and/or high missing rates (higher than 80% for unfiltered SNPs). One additional individual was removed due to very high levels of heterozygosity (0.24), which could neither be explained by introgression nor gene flow among divergent populations (see Table S2). Autosomal SNPs were retained if the genotyping rate was at least 45%. Furthermore, we filtered out potential paralogous regions by requiring that the normalized coverage was less than twofold the mean coverage of other SNPs (normalized by individual). The resulting SNP panel was called "general." Further filters were applied to produce different SNP panels aimed at maximizing the information content for different phylogenetic and population genetic analyses. For phylogenetic and comparative population genomic analyses among species, the SNP panel "multispecies" was produced.
For this, only SNPs with a genotyping rate of at least 50% within each species and 60% across species were retained. The SNP panel "intraspecies" was produced for all analysis performed for each species individually. For this, SNPs with a genotyping rate of >80% within a given species were retained. For all panels, monomorphic SNPs within a given group were removed (i.e., among all species together for "multispecies" or per species for "intraspecies"), as well as if the minor allele was found in only one individual (singletons and private doubletons).

| Phylogenetic and genetic diversity analyses among species
Two maximum-likelihood phylogenetic trees were built using RAXML (7.3.0, Stamatakis, 2006) based on a supermatrix of the "multispecies" SNP panel. The first tree explored the phylogenetic relationships among the three species and included one Alpine ibex, one Iberian T A B L E 1 Genetic diversity of domestic goat breeds, Alpine ibex, and Iberian Ibex populations. The sequencing coverage is shown as the median base pairs sequenced at restriction sites. The SNP density is shown as SNPs per kilobase sequenced. The estimated SNP density based on resampling corresponds to n = 10 individuals per population ibex of each of the two populations, and one domestic goat. Two individuals representing Iberian ibex were used for this first tree because a relatively high differentiation was expected between the two populations from a previous study (Angelone-Alasaad pers. comm.). A second tree explored the evolutionary relationship among all Alpine and Iberian ibex in our sample. We applied the GTR+gamma model for sequence evolution. A rapid bootstrap analysis including the search of the best maximum-likelihood tree was run with 100 replicates (Stamatakis, Hoover, & Rougemont, 2008). The trees were visualized using FigTree v1.4.2 (http://tree.bio.ed.ac.uk).
The R package {SNPRelate} was used to perform a principal component analysis (PCA) of all three species based on the genetic covariance matrix calculated from the genotypes (Zheng et al., 2012). For this, we used the SNP panel "multispecies" but retained only individuals genotyped at least at 80% of all SNPs (94 individuals). SNPs were then further filtered for a genotyping rate >80% and a minor allele frequency (MAF) of >5%. A total of 10,881 SNPs were retained. These additional filters were applied to avoid potential confounding effects of variable genotyping rates (see also Fig. S2C-F).

| Population structure within species
Principal component analysis as described above were also performed for each species individually. For this, the SNP panels "intraspecies" were used. As above, only individuals genotyped at least at 80% of all SNPs were retained (100 individuals) and SNPs were further filtered for a genotyping rate >80% and a minor allele frequency (MAF) of >5%. The function stat_ellipse() in the R package {ggplot2} was used to draw a 95% confidence ellipse for each population assuming a multivariate t-distribution.
For a complementary analysis of population structure, we used the Bayesian clustering method STRUCTURE (Pritchard, Stephens, & Donnelly, 2000). The same marker and individual set was used as for the PCA except for requiring a minimum of 250 kb pairwise distance between adjacent SNPs (3,167 SNPs left) to reduce potential nonindependence among loci. We varied K from K = 2 to K = 4 and ran for each K a burn-in of 10,000 and 50,000 Monte Carlo replicates using an admixture model with correlated allele frequencies.
For a quantitative assessment of the population structure of Alpine ibex, we performed a hierarchical AMOVA as implemented in the R package {poppr} with the method ade4. With the exception of two, all secondary reintroduced populations in this study were established from one primary reintroduced population. For the lower hierarchical level, we grouped each secondary reintroduced population with its primary population if there was such a simple 1-to-1 relationship. For the higher hierarchical level, we grouped all 1-to-1 pairs of reintroduced populations mentioned above and compared these with the source population (Gran Paradiso).

| Genetic diversity
The SNP density per population and per species, respectively, was calculated as the number of SNPs per kilobase sequenced with the SNP panel "general." As variation in sequencing coverage among populations could lead to variation in the number of SNPs passing the filters described above, we decided to use SNP density instead of SNP number to control for variation in sequence coverage. RAD-seq most consistently covers the region around the defined restriction enzyme cut sites. Hence, we restricted our genetic diversity analyses to these regions. Cut sites were identified in the goat genome by an in silico digest using restrict from EMBOSS (Rice, Longden, & Bleasby, 2000).

| Estimates of recent individual inbreeding
Runs of homozygosity (ROH) are defined as tracts of homozygosity along physical positions of a chromosome. ROH arise if both parents share identical chromosomal segments (i.e., segments of identity by descent, Lencz et al., 2007;Pemberton et al., 2012). Long ROH can be used to estimate the inbreeding history of individuals in a population (e.g., McQuillan et al., 2008). To identify ROH, we used the SNP panel "multispecies" but retained only individuals genotyped at more than 80% of the SNPs. Furthermore, only SNPs genotyped in ≥90% of the individuals were retained. The remaining dataset contained 41,907 SNPs and 94 individuals. This filtering reduced confounding effects of variable rates of missingness among populations (Fig. S4).
We identified chromosomal tracts composed of ROH using the option -homozyg in Plink (Purcell et al., 2007) as follows: We used a sliding window of 25 SNPs (homozyg-window-snp 25). The window was defined as homozygous if there was not more than one heterozygous site (homozyg-window-het 1) and not more than five missing sites (homozyg-window-missing 5). If at least 5% of all windows that included a given SNP were defined as homozygous, the SNP was defined as being in a homozygous segment of a chromosome (homozygwindow-threshold .05). This threshold was chosen to ensure that the edges of a ROH are properly delimited. A homozygous segment was defined as a ROH if all of the following conditions were met: the segment included ≥25 SNPs (homozyg-snp 25), covered ≥1,000 kb (homoz-kb 1,000), and contained a maximum of one heterozygous site (homoz-het 1). Furthermore, the minimum SNP density was one SNP per 50 kb (homozyg-density 50) and the maximum distance between two neighboring SNPs was ≤1,000 kb (homozyg-gap 1,000).

| Effective population size in recent past
We estimated effective population size using the LD-based method (Waples & Do, 2008) implemented in NeEstimator V2.01 (Do et al., 2014). SNP panel "general" was used with an additional missingness filter of ≤10% for a given population. The intermarker distance had to be at least 250 kb to avoid strong nonindependence of markers.

| RESULTS
We used RAD-sequencing to generate high-density, genomewide markers for Alpine ibex from the Gran Paradiso source population, the three primary Alpine ibex populations established during the reintroductions (Albris, Brienzer Rothorn, Pleureur) as well as subsequently established secondary populations (Figure 1a, Tables S1 and S2). We included Iberian ibex individuals of the subspecies C. pyrenaica hispanica from two populations: Maestrazgo Natural Park (north-eastern Spain) and Sierra Nevada Natural Space (southern Spain, Figure 1a).
Both populations were affected by bottlenecks but no extinction events ("autochthonous" populations). We also analyzed domestic goat from four Swiss breeds.
Using RAD-seq, we obtained overall read alignment rates varying from 95% to 96% among species. The more stringent alignment rates The median base pairs covered at restriction sites was 8.9, 8.3, and 7.6 Mb for Alpine ibex, Iberian ibex, and domestic goat, respectively.
The coverage amounted to 0.28%-0.33% of the total genome. The median SNP density was 2.3, 2.9, and 6.0 SNPs per kilobase sequenced for Alpine ibex, Iberian ibex, and domestic goat, respectively.

| Phylogenetic relationships and genetic differentiation among European ibex
As expected, maximum-likelihood phylogenetic analyses showed that Alpine and Iberian ibex were closely related (Figure 1b (35.1%) clearly separated the two populations and therefore explained considerably more variation than PC2 (4.7%). Individuals from the Sierra Nevada population were mostly differentiated by PC1, while individuals from the Maestrazgo population were differentiated by both PC1 and PC2. A PC analysis for domestic goats alone showed a clustering according to breeds (Fig. S2B).

| Genetic structure of reintroduced Alpine ibex populations
The  A Bayesian clustering analysis using STRUCTURE (Pritchard et al., 2000) confirmed the grouping identified by the PCA and identified the mixed ancestry of individuals from the ab and sm populations (Fig. S7).
The optimal K was K = 2 according to the delta K method (Evanno, Regnaut, & Goudet, 2005). The source population of all Alpine ibex (gp) was composed of multiple clusters. This is consistent with the fact that individuals of the source population were used to reintroduce all primary and secondary populations. Hence, clusters identified among reintroduced populations should also be represented in the source population.
The hierarchical AMOVA mirrored the results of PCA and STRUCTURE analysis and found significant differentiation among the groups (comprising the respective primary and secondary reintroduced populations) and among all individual populations (Table S3). However, no significant differentiation (p = .26) was found between source and reintroduced populations.

| Genetic diversity in populations
The Gran Paradiso Alpine ibex population and the Iberian ibex populations give insights into levels of genetic diversity in populations of ibex that have experienced bottlenecks but not reintroductions.
Comparing the genetic diversity of the Gran Paradiso population with the reintroduced populations gives insights into the impact of reintroduction bottlenecks. SNP density (SNPs per kilobase sequenced) was similar in the Gran Paradiso population and the two Iberian ibex populations. However, SNP density was lower in all reintroduced Alpine ibex populations ( Figure 3). As estimates of genetic diversity can be affected by sampling effort, we also used a resampling procedure to compare genetic diversity among groups. We found higher SNP densities in the autochthonous Gran Paradiso population and Iberian ibex than in the pool of the reintroduced populations. In comparison, domestic goat had the highest SNP density (Fig. S8A). Based on resampling 12 individuals, we found a median SNP density of 1.3 SNPs per kilobase of sequence across all reintroduced Alpine ibex populations, of 1.6 and 1.7 SNPs/kb in the two autochthonous populations of Iberian ibex, and of 1.8 SNPs/kb in Alpine ibex from Gran Paradiso ( Figure 4). Resampling curves of SNP densities estimated for individual populations confirmed higher SNP densities in primary than in secondary reintroduced populations (Fig. S8B).
Expected heterozygosity was highest in the two Iberian ibex populations (0.039 and 0.040 for Sierra Nevada and Maestrazgo, respectively), intermediate in Gran Paradiso (0.032) and lowest in the reintroduced populations (0.023-0.026; filled circles in Figure 3b and Table 1). Domestic goats had the highest expected heterozygosity.

| Comparisons of effective population sizes
Estimates of effective population sizes revealed that all secondary reintroduced Alpine ibex populations had a N e of 40-120. The Iberian ibex population Maestrazgo had a comparable N e ( Figure 5 and Table   S4)  Table S4 for N e estimates based on different MAF thresholds). Note that Waples and Do (2008) showed that N e estimations performed better for MAF > 0.1 than for lower thresholds (i.e., higher percentage of confidence intervals included the true Ne).
Although the upper bounds of the estimates were often infinite, comparisons of the N e estimates of the secondary reintroduced populations with observed levels of genetic differentiation among populations suggest that the point estimates of N e are reasonable.
Assuming no differentiation at the time of reintroduction (i.e., assuming F 0 = 0), pairwise F ST is predicted to increase over time according (Wright, 1969)

| Individual-based levels of recent inbreeding
We estimated inbreeding among ibex individuals using runs of homozygosity (ROH) along the genome (Figure 6). Using stringent filters, we predominantly identified ROH in regions of high SNP densities (Fig. S9). Long ROH indicate recent inbreeding (see Figure 7 for examples). Total lengths of all ROH including short ROH, which are not indicative of recent consanguinity, were largely identical among ibex populations (Fig. S10A). However, the total length of long ROH (>5 Mb), which are indicative of recent inbreeding, varied among populations (Fig. S10B). In Alpine ibex, we found that the population median of total length of long ROH (>5 Mb) was at least 216 Mb (Gran Paradiso) and covered at least 9% of the autosomal genome (Fig. S10B). Finally, we correlated alternative estimates of inbreeding with ROH. The molecular inbreeding coefficient FROH 10,000 is defined as the total length of ROH longer than 10 Mb divided by the autosomal genome length (Purfield et al., 2012). Multilocus homozygosity and FROH 10,000 were positively correlated across individuals (R 2 = .29, p < .0001, n = 92) (Fig. S11). Biebach and Keller (2010) used population-specific F ST to estimate population inbreeding. The F ST used for these estimates reflected the differentiation of the focal population compared to all other populations. We found a weak but significant positive correlation between FROH 10,000 and population inbreeding estimated based on F ST (R 2 = .047, p = .04, Fig. S12).

| DISCUSSION
Sustainable species conservation depends on understanding the genetic consequences of management strategies. Here, we described

| The application of RAD-seq to recapitulate demographic processes in a reintroduced species
Our analyses were based on high-throughput RAD-seq genotyping.
High-throughput genotyping such as RAD-seq and other barcoded next-generation sequencing techniques can introduce potential biases in the genotypic dataset. We performed a series of analyses to identify issues arising specifically from barcoded sequencing libraries that combine a large number of individuals of potentially different DNA quality. We found that uneven coverage among individuals (and populations) can introduce biases if these are not properly controlled for. For example, sequencing coverage was inversely correlated with genotype missingness (data not shown). Hence, we required a minimum coverage for any individual to be included in the analyses. We chose SNP density as a measure of genetic diversity rather than raw SNP counts to control for remaining variation in per-individual coverage. Variation in the proportion of genotyped markers (i.e., genotyping rate) can influence measures of heterozygosity and population differentiation (e.g., Fig. S3). Therefore, we applied strict missingness filters to include only the best-sequenced SNPs in such analyses. We also applied strict filters for marker density to produce conservative estimates of runs of homozygosity (ROH). Ascertainment bias may further be introduced by the fact that reads were aligned to a single reference genome from the domestic goat. However, we found that reads obtained from different species aligned nearly equally well to the reference genome. The ascertainment bias from aligning to the goat genome should therefore be small. Interestingly, the alignment rates of sequences obtained from domestic goat were lower on average than sequences obtained from ibex species. This is despite the fact that the reference genome was built from a domestic goat. The lower alignment rates of domestic goats are most likely explained by the substantial genetic diversity among domestic goat (samples representing Swiss breeds mapped to reference built from Chinese breed).
High-throughput sequencing protocols for reduced representation libraries such as RAD-seq include PCR cycles to amplify the small pool of primary DNA fragments. PCR produces duplicate fragments, which can be retained in the final pool of sequenced fragments. We found indeed PCR duplicates among the reads obtained from each individual. The proportion of duplicates was in the range of typical RAD-seq studies (for a review see Andrews, Good, Miller, Luikart, & Hohenlohe, 2016). PCR duplicates should not lead to systematic (i.e., genomewide) genotyping errors. However, PCR duplicates can be an issue if high genotyping accuracy is required at individual loci (e.g., for genomewide association studies). As our analyses were based on genomewide estimates of population genetic parameters, we retained PCR duplicates in our analyses. We found that the presence of PCR duplicates had no discernable impact on estimates of population differentiation (see Fig. S5).

| The impact of demography on reintroduced populations
The primary aim of our study was to quantify the loss of genomewide diversity due to the reintroduction of Alpine ibex from a single source population (Gran Paradiso). We found strong evidence for withinspecies substructure that reflects the reintroduction strategy.

| Comparative population genomic analyses of European ibex and domestic goat
Alpine ibex suffered a unique population bottleneck compared to other ibex species due to the near extinction in the 19th century. However, the Iberian ibex occupy an alpine habitat that is similarly fragmented by topography. Both autochthonous populations of Iberian ibex and the source population of Alpine ibex were confronted with strong barriers to dispersal. At high population densities, some barriers might be overcome due to high dispersal pressure. Such dispersal may lead to the colonization of new territories or to population expansion (Refoyo, Olmedo, Polo, Fandos, & Muñoz, 2015;Refoyo et al., 2016). However, dispersal barriers between the two Iberian ibex populations seem to be effective and demographic bottlenecks are evident leading to a high differentiation from each other. The genetic differentiation is higher than the most differentiated Alpine ibex populations (pairwise F ST = 0.39 vs. This may seem surprising because most domesticated plant and animal species suffered from considerable losses in genetic diversity (e.g., Blackman et al., 2011;Renaut & Rieseberg, 2015;Ross-Ibarra, Morrell, & Gaut, 2007;Zeder, Emshwiller, Smith, & Bradley, 2006), deleterious mutations accumulation (e.g., sunflower: Renaut & Rieseberg, 2015;dogs: Marsden et al., 2016), and inbreeding (e.g., Marsden et al., 2016). However, the domestication of goats was likely more complex than in most other domesticated animals as the extant domestic goat may have multiple origins of domestication (Luikart et al., 2001). In addition, neither of the two ibex species under study are direct sister species of the domestic goat. Genetic diversity data of the wild sister species of domestic goat, bezoar (Capra aegagrus), may provide further insights into the origins of the domestic goat diversity.

| Genetic diversity in autochthonous ibex populations
All populations in this study, including Iberian ibex populations, suffered from bottlenecks. Nevertheless, we observed that the levels of genetic diversity estimated by SNP density and multilocus heterozygo- Although SNP density was comparable between Gran Paradiso and the Iberian ibex, expected heterozygosity was considerably lower in Gran Paradiso. This discrepancy may be explained by a higher proportion of low-frequency polymorphism in the Gran Paradiso (Fig. S13).
Population genetics theory predicts that severe bottlenecks over few generations can have a similar impact on genetic diversity to mild bottlenecks over longer periods. This may explain why we observed levels of genetic diversity among the Iberian ibex and Gran Paradiso populations that did not match predictions made from demographic information about bottlenecks. Furthermore, census data reported for specific bottlenecks can be a poor predictor of effective levels of genetic diversity. This may be due to the fact that the reported data were inaccurate or that additional demographic information, such as the number of reproducing individuals and the occurrence of immigration, is needed.
Autochthonous populations (Gran Paradiso and Iberian ibex) had consistently higher genetic diversity than reintroduced populations.
One Gran Paradiso individual showed extraordinarily high levels of heterozygosity (>0.08, Figure 3b), which increased the SNP density of the Gran Paradiso population. However, even after excluding the individual, SNP density was still higher than the pool of reintroduced populations (Fig. S8C). This suggests that reintroductions to the Swiss Alps led to consistent losses in genetic diversity. A previous study based on microsatellites was unable to resolve these differences between the source and the pool of reintroduced populations (Biebach & Keller, 2009).

| Runs of homozygosity to accurately assess levels of recent inbreeding
A major concern for the management of reintroduced species is that the reduced levels of genetic diversity cause harmful inbreeding.
However, detecting early inbreeding in natural populations can be challenging. Traditionally, inbreeding was estimated using pedigree data (which are difficult to obtain from wild populations) or multilocus heterozygosity. The affordability of RAD-seq or genotyping by sequencing led to the implementation of genomewide inference in a number of recent conservation studies despite considerable challenges (Shafer et al., 2015). Genomewide markers allow powerful estimates of inbreeding based on runs of homozygosity (ROH), which have so far mostly been performed in studies of human and domestic animals (e.g., Lencz et al., 2007;McQuillan et al., 2008;Pemberton et al., 2012;Purfield et al., 2012). A key requirement for ROH estimates is knowledge of physical marker positions. Such information is often only available if a reference genome of the same or a closely related species exists. The rapid increase in the number of available genome sequences from wild species should facilitate ROH analyses in many species of conservation concern. Although pedigree, ROHbased inbreeding estimates, and multilocus homozygosity are often correlated (e.g., Purfield et al., 2012) (Purfield et al., 2012). Furthermore, we only considered regions for ROH with above-average SNP densities. Therefore, we may have underestimated the extent of ROH in certain genomic regions (see Fig. S9).
All individuals in our study were found to have ROH up to 5 Mb.
The total length of ROH was similar among all Alpine and Iberian ibex populations and comparable to estimates based on whole genome sequencing in other mammals. In inbred Gorilla populations, ROHs of 2.5-10 Mb covered between 8% and 20% of the total genome (Xue et al., 2015). ROH tracts of the same length category covered 14%-17% of the total autosomal genome length in the two ibex species.
ROH tracts longer than 20 Mb were only observed in the reintroduced Alpine ibex populations. Hence, recent inbreeding was more frequent in the reintroduced populations than the Gran Paradiso source population. We observed individual ROH-based inbreeding estimates (FROH 10,000 ) up to 0.1. Values in this range were shown to correspond to pedigree inbreeding estimates of matings among half sibs (pedigree F > 0.125; Purfield et al., 2012). The high level of inbreeding in reintroduced Alpine ibex populations suggests that inbreeding depression may be expressed. Indeed, a study of Alpine ibex in the Gran Paradiso revealed inbreeding depression in several fitness-related traits (Brambilla et al., 2015).

| Applying genomic surveys to assess the genetic status of reintroduced populations
Many endangered species suffered from population bottlenecks in their recent history. Species with a significantly diminished distribution range can be conserved by population reintroductions. However, population reintroductions repeatedly failed to meaningfully extend lost species ranges (Armstrong & Seddon, 2008). Despite the large number of population reintroductions, few studies investigated the long-term genetic consequences of reintroductions. Such consequences critically include the potential loss of genetic diversity and increased inbreeding levels due to genetic bottlenecks. Following bottlenecks, populations may have a reduced potential to adapt to a changing environment or resistance to new diseases (O'Brien & Evermann, 1988). The population genetic consequences of bottlenecks were convincingly demonstrated using microsatellites and mitochondrial markers in multiple species (e.g., Biebach & Keller, 2009Diefenbach et al., 2015;Jamieson, 2011;Strzała et al., 2015).
However, only few studies included genetic data from the source population and monitoring data of the entire population reintroduction history. Using genomewide inference, our study showed that population reintroductions of the Alpine ibex lowered genetic diversity and increased recent inbreeding. Both Iberian ibex populations and the Alpine ibex source population maintained higher levels of genetic diversity and lower levels of recent inbreeding. However, our estimates of genetic diversity in these autochthonous populations did not correspond to census information on historic population bottlenecks. Census information can be a poor predictor of effective levels of genetic diversity, because reported data may be inaccurate and not fully informative about demography (e.g., migration and sex ratios).
Hence, genomic surveys serve as an important monitoring tool even if census data are available. Consideration of the population genomic consequences of conservation efforts should become an integral component of managing reintroductions.

| Management recommendations
Our study identified low genetic diversity and high inbreeding levels in reintroduced populations of Alpine ibex. Effects were even more pronounced in secondary vs. primary reintroduced populations. Low genetic diversity might reduce fitness in the long-term due to reduced evolutionary potential (Allendorf et al., 2012;Biebach, Leigh, Sluzek, & Keller, 2016). Inbreeding levels were in a range where inbreeding is expected to impact fitness (e.g., Huisman, Kruuk, Ellis, Clutton-Brock, & Pemberton, 2016;Keller & Waller, 2002). Thus, existing Alpine ibex populations would benefit from introducing individuals from genetically differentiated populations to increase genetic variation and reduce inbreeding. Such measures may be necessary to preserve the long-term survival of Alpine ibex. Results of this study show that secondary reintroductions from a sole primary reintroduced population should be avoided. Instead, new populations should be founded with individuals from the source population (GP) or from multiple primary reintroduced populations. This is in line with a previous study with microsatellites where the degree of admixture of the founder group had a higher impact on genetic variation than the number of founders .
The re-establishment of Alpine ibex across the Alps was one of the most successful species reintroductions. Alpine ibex recovered from about 100 individuals to ca. 50,000 individuals across the Alps.
Our study showed that although demographic risks are currently minimal, genetic (i.e., inbreeding) risks due to the reintroduction history persist. Therefore, future reintroductions should be genetically monitored even after the successful establishment of a population, with a particular focus on early signs of inbreeding. Such monitoring would ensure the early detection of genetic problems of endangered species and would allow initiating the necessary management actions.