Ancient mitochondrial and modern whole genomes unravel massive genetic diversity loss during near extinction of Alpine ibex.

Population bottlenecks can have dramatic consequences for the health and long-term survival of a species. Understanding of historic population size and standing genetic variation prior to a contraction allows estimating the impact of a bottleneck on the species' genetic diversity. Although historic population sizes can be modelled based on extant genomics, uncertainty is high for the last 10-20 millenia. Hence, integrating ancient genomes provides a powerful complement to retrace the evolution of genetic diversity through population fluctuations. Here, we recover 15 high-quality mitogenomes of the once nearly extinct Alpine ibex spanning 8601 BP to 1919 CE and combine these with 60 published modern whole genomes. Coalescent demography simulations based on modern whole genomes indicate population fluctuations coinciding with the last major glaciation period. Using our ancient and historic mitogenomes, we investigate the more recent demographic history of the species and show that mitochondrial haplotype diversity was reduced to a fifth of the prebottleneck diversity with several highly differentiated mitochondrial lineages having coexisted historically. The main collapse of mitochondrial diversity coincides with elevated human population growth during the last 1-2 kya. After recovery, one lineage was spread and nearly fixed across the Alps due to recolonization efforts. Our study highlights that a combined approach integrating genomic data of ancient, historic and extant populations unravels major long-term population fluctuations from the emergence of a species through its near extinction up to the recent past.

INTRODUCTION

The ongoing crisis of biodiversity loss is often referred to as the human‐caused sixth mass extinction (Ceballos et al., 2015, 2020). However, not only are species disappearing at an increasingly fast rate, anthropogenic pressures on population size and connectivity affect many more species, which are currently still at fairly good numbers. Hence, it becomes crucial to understand the scale and the genetic consequences of small population size as well as population fragmentation in the wild. Theory predicts progressive loss of genetic diversity, increased inbreeding and accumulation of deleterious mutations in small and isolated populations (Frankham et al., 2002; Hedrick & Garcia‐Dorado, 2016). Long standing theoretical and empirical studies indicate the risk for reduced long‐term survival and adaptive evolvability caused by low genetic diversity (Frankham, 2010; Hasselgren & Norén, 2019; Hedrick & Garcia‐Dorado, 2016; Keller & Waller, 2002; Wright, 1921). In this context, it is interesting that currently the protection status of a species does not necessarily predict its level of genetic diversity (Díez‐Del‐Molino et al., 2018; Grossen et al., 2020, but see also Canteri et al., 2021). Furthermore, current genetic patterns may have different explanations. Low genetic diversity is usually explained by human induced population fragmentation and a sudden strong bottleneck, but also a long history of small population size and restricted gene flow may be its cause. Disentangling these two scenarios is important for proper species conservation measurements. Hence, this underlines the importance of taking into account the past demographic history of a species to predict its long‐term viability. However, past bottlenecks (times of small population size) are not trivial to estimate. Recent methods based on coalescence theory such as pairwise sequentially Markovian coalescent (PSMC) (Li & Durbin, 2011) and more recently multiple sequentially Markovian coalescent (MSMC) (Schiffels & Wang, 2020) have been widely used to estimate the demographic history of species using contemporary DNA and investigate the impact of environmental changes on trajectories of past population size (Kozma et al., 2016, 2018; Palkopoulou et al., 2015; Pečnerová et al., 2021). However, the study of recent DNA is expected to be limited in power due to uncertainty of past diversity, most of all after strong bottlenecks leading to significant loss of signal. Furthermore, estimates for the recent past, which was mostly affected by human presence on Earth, are not possible (Li & Durbin, 2011; Nadachowska‐Brzyska et al., 2016).

Ancient genomics is a promising state of the art method to solve this issue (Díez‐Del‐Molino et al., 2018). The study of ancient DNA (aDNA) allows quantifying population size over time and retracing changes of diversity through near extinction events thereby shedding light on the most important factors shaping current genetic diversity patterns. As expected, several aDNA studies focusing on prebottlenecked diversity suggested a major role for human impact on current patterns of genetic diversity (e.g., Casas‐Marce et al., 2017). The Iberian lynx (Lynx pardinus) for instance, displayed a clear reduction in population size due to overhunting and divergence into two distinct subpopulations with drops to alarmingly low genetic diversity within a century (Casas‐Marce et al., 2017). However, humans may not always be to blame. Studies on Musk ox (Ovibos moschatus) and Kea (Nestor notabilis) interestingly identified environmental changes as the most important drivers of diversity loss suggesting prolonged histories of small population size (Campos et al., 2010; Dussex et al., 2015), while Larsson et al. (2019) revealed a complex interplay of climatic and anthropogenic factors in Arctic fox (Vulpes lagopus). There is also contradicting evidence and hence an ongoing debate on the role of humans in the Late Quaternary megafauna mass extinction (Campos et al., 2010; Koch & Barnosky, 2006; Lord et al., 2020; Sandom et al., 2014; Stewart et al., 2021). Due to their imminent exposure to changing glaciation levels, demographic histories of species living in arctic or alpine habitats are expected to have been considerably affected by past climatic changes (Sommer, 2020). While recent aDNA studies have mainly shed light on extinct arctic or alpine species (Barlow et al., 2020; Gretzinger et al., 2019; Lord et al., 2020; Palkopoulou et al., 2015; Ramos‐Madrigal et al., 2021), fewer studies have investigated alpine species which survived the megafauna extinction event of the Pleistocene–Holocene boundary (e.g., Ureña et al., 2018).

An alpine species with currently low genetic diversity, high mutation load, high levels of inbreeding and signs of inbreeding depression, is the Alpine ibex (Capra ibex) (Biebach & Keller, 2009; Bozzuto et al., 2019; Brambilla et al., 2018; Grossen et al., 2018; Grossen et al., 2020). Historic records suggest that Alpine ibex were intensely hunted, presumably since the 15th century and encountered their most severe bottleneck in the 19th century. The species survived near extinction in a single small population in a region which is today known as the Gran Paradiso National Park in Italy (Stüwe & Nievergelt, 1991) before intense conservation efforts led to a fast recovery of the species during the 20th century. As predicted by theory, the severe bottleneck of close to 100 individuals left clear genetic footprints in contemporary populations (Biebach & Keller, 2009; Grossen et al., 2018). Yet, the eradication of the species from almost its entire species range also eradicated signals of the past demography and it remains currently unclear what level of diversity was present before the near extinction. Because of the specific ecological needs of Alpine ibex (Grignolio et al., 2003, 2007), past populations may also have been small and isolated. Given its current distribution in previously glacier‐covered habitats, it is plausible that environmental changes in the late Pleistocene and Holocene also played a substantial role in shaping genetic patterns (Seersholm et al., 2020).

Here, we combine diversity estimates covering the last eight millennia to obtain insight into the demographic and genetic history of the nearly extinct Alpine ibex. Taking advantage of published whole genome sequences, we quantify the effective population size of Alpine ibex over time and compare our observations with related species. The analysis of ancient and historic mitogenomes allows retracing changes of diversity and the demographic trajectory through the species near extinction. We quantify the current haplotype diversity in the light of past diversity, identify lost haplotypes and most affected gene regions. Finally, we investigate to what degree demographic modelling based on current genetic data can recapitulate insights gained from prebottleneck aDNA sampling.

MATERIALS AND METHODS

PSMC

To reconstruct the demographic trajectories of the Alpine ibex during the late Pleistocene, we incorporated a PSMC approach (Li & Durbin, 2011; Nadachowska‐Brzyska et al., 2016). The PSMC infers historical recombination events within a diploid genome and facilitates differences in heterozygosity within one individual. It has improved statistical strength to infer deep‐in‐past events of coalescence by inferring the most recent common ancestor (TMRCA) and thereby the ancestral effective population size, depending on generation time, over the last 2*103 to 3*106 years (Li & Durbin, 2011; Mather et al., 2020; Nadachowska‐Brzyska et al., 2016). We used whole genome data, which has previously been published (Alberto et al., 2018; Grossen et al., 2020), representing domestic goats (Capra hircus,N = 3), Bezoar (Capra aegagrus,N = 2), Nubian ibex (Capra nubiana,N = 2), Sibiran ibex (Capra sibirica,N = 2), Iberian ibex (Capra pyrenaica,N = 2), and six Alpine ibex from the source population in Gran Paradiso, Italy. The whole genomic data of the domestic goat, Bezoar and domestic sheep are available through the NextSeq Consortium (see also Alberto et al., 2018; Grossen et al., 2020). The data of the ibex species were obtained by Grossen et al. (2020). The reads were trimmed using trimmomatic v.0.36 (Bolger et al., 2014) and subsequently mapped with bwa‐mem v.0.7.17 (Li, 2013) to the domestic goat reference genome (ARS1, Bickhart et al., 2017). Duplicated reads were marked with MarkDuplicates from picard v.1.1301 and hence ignored for all further analysis. The mean genome‐wide coverage was >99.38% for all samples and the average depth ranged from 6.9× to 46× (Table S1). To produce the input data set for the PSMC analysis, we followed the general pipeline suggested in Palkopoulou et al. (2015). In detail, we used samtools mpileup (Li et al., 2009) in combination with the bcftools call command, keeping reads with a minimum mapping quality (‐q) and minimum base quality (‐Q) of 30 to produce an alignment. Next, we called consensus sequences using bcftools ‐c and performed a final filtering step using ‘vcf2fq’ from vcfutils.pl with options ‐d 5 and ‐D 34 (minimum and maximum coverage) and ‐Q 30 (minimal mean squared mapping quality). This pipeline has the advantage that the aligner does not assume Hardy–Weinberg equilibrium and does not rely on population frequencies for variant calling (Nadachowska‐Brzyska et al., 2016). We then used fq2psmcfa to produce input data for psmc v.0.6.5‐r67, which was run with standard parameters as suggested by Li (2016). Specifically, the limit of TMRCA and the maximum number of iterations were left at the default values ‐t 15 and ‐N 25. N e was inferred across 64 free atomic time intervals. The initial population size parameter (option ‐p) was set to “4 + 25*2 + 4 + 6”. This corresponds to four atomic time intervals followed by 25 parameters spanning two intervals followed by two parameters spanning four and six intervals, respectively and allowed for 28 (1 + 25 + 1 + 1) free interval parameters. The psmcfa output was visualized in r, using a modified version of the plotPsmc.r script supplied by Liu and Hansen (2017) with mutation rates of 2.23e‐09 sites/year (1.784e‐08 sites/generation) inferred for the Siberian ibex (Chen et al., 2019). To explore the effect of different mutation rates on the demographic trajectories, we also doubled and halved the mutation rate for the PSMC analysis (Figure S1A and B).

Dating of split between Alpine and Iberian ibex

We took advantage of the whole genome data to date the split time between the Alpine (N = 29) and Iberian ibex (N = 4) and get another estimate of the long‐term N e of both species by facilitating a fast sequential Markov coalescent simulation using the tool fastsimcoal2 v.2.7.05 (Excofffier et al., 2021). To this end, we used the samtools model (‐GL 1) to estimate the genotype likelihood at variable sites for all autosomes and calculated the full site frequency spectra (‐dosaf 1) for the derived alleles using the goat reference genome (ARS1, Bickhart et al., 2017) as reference and ancestral (‐ref and ‐anc). We filtered out reads with base quality lower 20 (‐minQ 20), mapping quality lower than 30 (‐minq 30) and reads with two perfect matches (‐uniqueOnly) and also removed bad reads (‐remove_bads 1). We processed the resulting saf.idx with angsd's realSFS (Korneliussen et al., 2014) to create the site frequency spectra (SFS). The realSFS output was modified using a custom‐made script to the SFS format required by fastsimcoal2. The model was set up to estimate the effective population size of the Iberian ibex (predefined range: 4000 to 400,000), the Alpine ibex (predefined range: 20 to 20,000) as well as the time of divergence between these two groups (predefined range: 2500 to 60,000 generations). Furthermore, we selected the FREQ data type. We then ran 100 independent simulations with the setting: ‐n 200,000, ‐L 50, ‐d, ‐M. We chose the best run based on the smallest AIC and highest maximum likelihood value. Next, we followed the fastsimcoal2 manual to create non‐parametric bootstraps. With fastsimcoal2 we generated 100 SFS each based on 200,000 non‐recombining DNA segments of length 1000 bp using the options: ‐n 1, ‐b99, ‐j, ‐d, ‐s0, ‐l, ‐q. The generated SFS were then used for the parameter estimation in fastsimcoal2 with the settings: ‐L 30, ‐initvalues, ‐n 100,000, ‐d, ‐M.

Sampling

To analyse the ancient genetic diversity of Alpine ibex, we identified and collected samples which originated from before and during the last strong species bottleneck during the 19th century. We obtained a total of 22 (sampled in 2017) and four (sampled in 2020) specimens from Swiss museums, archaeological institutions and cave excavations (Table S2). After screening for endogenous DNA content, a total of 15 samples were chosen for subsequent analysis. Samples originating from the last millennium are from now on referred to as historic samples. Their age (here given as CE, Common era) was inferred by consulting registry entries of the museums of origin (except for one C14‐date, Table S2) and ranged from 1000 CE to 1919 CE. These specimens originated from Italy, France and Switzerland. The ancient samples were AMS‐C14 dated and hence their age is given as BP (before present defined as 1950). They ranged from 8601 ± 33 BP to 3302 ± 26 BP and were found in caves, except one specimen (Sample ID: Gro1) found in a glacier field in the Austrian Alps (7212 ± 27 BP) (Table S2).

Since whole genome sequence data allows assembling mitogenomes, we took advantage of previously published whole genome sequencing data to compare pre‐ and postbottlenecked genetic diversity of the species. The downloaded data represented recent populations of Alpine ibex (N = 29), Iberian ibex (N = 4), Nubian ibex (N = 2), Markhor (C. falconeri,N = 1), Siberian ibex (N = 2), Bezoar (N = 6) and domestic goat (N = 16) (Alberto et al., 2018; Grossen et al., 2020). Additionally, previously published whole genome sequencing data of five sheep individuals representing four species (Ovis sp.,Table S2, Alberto et al., 2018) were used as an outgroup for the phylogenetic analysis.

We produced a detailed map displaying the origin of recent, historic and ancient Alpine ibex specimens in qgis v.3.0.2. For specimens where the exact location was not known (most of the historic samples), assigned coordinates are based on the region of origin (Table S2).

Sample collection and DNA extraction

To ensure the accuracy of the aDNA results, established quality standards were incorporated (Cooper & Poinar, 2000). To account for the usually minute quantity and degraded state of DNA in ancient samples, special care in sample handling was taken. (1) DNA extraction and library preparation were performed in a dedicated aDNA laboratory, (2) a one‐way work flow from pre‐PCR to post‐PCR laboratory was integrated, (3) blank controls were used, (4) and used equipment was decontaminated with 7% sodium hypochlorite/NaClO and/or UV‐irradiation.

For all sampling runs, DNA was extracted in a specialized clean‐laboratory at the University of Zurich. In detail, we extracted DNA from teeth, skulls, long bones, petrous bones and horn material. During the sampling in 2017, surface sterilization of the samples was performed by washing the fragments in 1% sodium hypochlorite/NaClO for 1 min, followed by three washing steps with ddH2O. We extracted DNA from teeth, by detaching the cervix from the crown of the tooth with a precision drill (Dremel 8200) at a low rotation rate of c. 7000 rpm (speed level 10). To minimize friction, diamond cutting blades (Dremel SC545) were used. Approximately 0.75 cm3 (horn) or 1 cm3 (bone or tooth material) were extracted with the Dremel. All 2017 samples were then ground to bone powder with an analytical mill (IKA A11 basic) and stored in micro test tubes at room temperature and under absence of light. For the 2020 sampling run, we sterilized the bone and horn surface with UV for 30 min per side, cleaned the surface with 7% sodium hypochlorite/NaClO and removed the most outer part of the bone or horn with a dental drill and tungsten steel drill bits (Alpine Orthodontics, H1‐014‐HP). After removing an initial surface layer of bone, we extracted bone powder at low rotation rate.

To extract and purify mitochondrial and nuclear DNA from the bone and horn powder, we used a QIAquick PCR purification kit and applied an established protocol for DNA extraction as described in Dabney et al. (2013). The following modifications were made: we lysed 100 mg bone powder by adding 50 μl proteinase K, 800 units/ml (New England Biolabs Inc.) to 950 μl EDTA and incubated at 37°C overnight. The overnight digestion was centrifuged at 14 krpm for 5 min in a table top centrifuge, mixed with 10 ml binding buffer (without tween) and 400 μL 3 M sodium acetate. We centrifuged the resulting suspension at 1500 rpm for 15 min. We stored the eluted samples at 4°C. To get a first impression of sample quality of the 2017 samples, we performed a PCR‐based screening on two microsatellites and two mitochondrial regions and only kept samples for which at least one PCR product was visible for further analysis (Appendix S1, Table S2).

Library preparation and sequencing

Double‐stranded DNA libraries were built for all samples under strict precautions to avoid contamination. All preamplification steps for constructing the libraries were performed in a dedicated aDNA clean‐laboratory and nontemplate controls were included. For the 2017 samples, initial libraries were constructed for screening purposes at the Institute of Evolutionary Medicine at the University of Zurich following Kircher et al. (2012) and sequenced on an Illumina HiSeq 2500 system (1 lane, 200 cycles, paired‐end) at the Functional Genomics Center Zurich (FGCZ). The resulting data was mapped to the domestic goat reference genome CHIR1 (Dong et al., 2012). Samples with very low endogenous DNA content (mapping rate of <4%) were excluded for further analysis (Table S2). After screening, another set of libraries were built as individually double‐indexed libraries optimized for aDNA (Kircher et al., 2012; Meyer & Kircher, 2010), but here notably we corrected for post‐mortem‐damage by adding the uracil‐specific excision reagent (USER, New England Biolabs Inc.) during the blunt end repair step (for further detail see Appendix S1). The libraries were constructed in facilities of the Swedish Museum of Natural History in Stockholm, Sweden and deep sequenced on an Illumina HiSeq X system (one lane per individual, 300 cycles, paired‐end) in a facility of the National Genomics Infrastructure (NGI), Sweden. Double‐stranded DNA libraries for the 2020 samples (N = 4) were built in a specialized aDNA clean‐laboratory at the Institute for Evolutionary Medicine, Zurich (Appendix S1) following the protocol described in (Kircher et al., 2012; Meyer & Kircher, 2010). The sequencing was performed on one run of an Illumina NexSeq‐500 in mid‐output mode with 150 cycles (2*75 + 8 + 8) and sequenced at the Functional Genomics Center Zurich (FGCZ).

Raw data analysis

Trimming of adapter sequences and read quality filtering was performed using trimmomatic v.0.36 (Bolger et al., 2014) with the following settings: ILLUMINACLIP:2:30:10:1:TRUE, LEADING:3, SLIDINGWINDOW:4:15, MINLEN:25. Trimmed reads were mapped using the Burrows‐Wheeler Aligner (Li & Durbin, 2009) bwa‐mem v.0.7.15‐r1142 to the Capra ibex mitochondrial reference genome NC_020623.1 (NCBI). Reads were sorted using samtools v.1.10 (Li et al., 2009). Duplicates were identified using MarkDuplicates from Picard (http://broadinstitute.github.io/picard, v.2.8.3). Summary statistics were produced using samtools v.1.10. Post‐mortem‐damage of the 2017 samples was already confirmed at the screening step (Figure S2) and the samples were USER‐enzyme treated for the deep sequencing. Hence, no further post‐mortem assessment was carried out for these libraries. Post‐mortem‐damage of the 2020 samples (Figure S2) was assessed with mapdamage v.2.0 (Jónsson et al., 2013). Base quality scores were corrected for post‐mortem‐damage by running mapDamage ‐‐rescale‐only for the paired and unpaired reads separately. The published data representing all recent samples was also quality trimmed in trimmomatic v.0.36 (Bolger et al., 2014) with settings ILLUMINACLIP:2:30:10 LEADING:5 TRAILING:5 SLIDINGWINDOW:4:15 MINLEN:50 and bam files were generated as described above.

We used HaplotypeCaller from gatk v.4.1.7 (Van der Auwera and O'Connor) to discover variant sites and produce g.vcf files for all recent, historic and ancient samples, before combining them and perform joint genotyping using GenomicsDBimport and GenotypeGVCFs. gatk Variantfiltration was used for hard filtering applying the following criteria to retain a site: Quality by depth (QD) >2.0, mapping quality (MQ) >20.0, mapping quality rank sum (MQRankSum) >−3 or <3, Fisher Strand (FS) <40.0, strands odds ratio (SOR) <5.0, ReadPosRankSum >−3.0 and <3.0. We furthermore removed low quality genotypes with a genotype quality (GQ) below 20 (vcftools, Danecek et al., 2011).

We produced two mitogenome data sets, which differed in the number of individuals/species (Capra_all and Capra_ibex). For both data sets, we required a minimal genotyping rate per site of 90% (vcftools, Danecek et al., 2011). Hence, also the number of sites included differed between the two. The data set Capra_all contained 80 individuals representing Alpine ibex (N = 44: 29 recent, 8 historic, 7 ancient), Iberian ibex (N = 4), Bezoar (N = 6), Siberian ibex (N = 2), Markhor (N = 1), Nubian ibex (N = 2), domestic goat (N = 16) andsheep (Ovis sp., N = 5). This data set contained a total of 15,986 of a possible 16,157 sites. We furthermore produced a second data set Capra_ibex, which only contained Alpine ibex specimens (N = 44: 29 recent, 8 historic, 7 ancient). It was composed of 16,135 out of a total of 16,157 sites.

Phylogeny

To explore the phylogenetic relationship among pre‐ and postbottlenecked Alpine ibex and related species, we built a maximum likelihood tree. We also used it to confirm the species of our specimens, because the unambiguous morphological identification of small remains of Alpine ibex in relation to other Capra species can be challenging. The tree was built using the data set Capra_all. We used vcf2phylip to convert VCF to PHYLIP and defined all five sheep as the outgroup. We used the program raxmlv.8.2.10 with the ‐GTRGAMMA model, which is a GTR model of nucleotide substitution under the Γ‐model of rate heterogeneity (Stamatakis, 2014). We performed a rapid bootstrap analysis with 100 repetitions and a search for the best‐scoring maximum likelihood tree.

For a more in‐depth analysis of the relationships among the ancient and historic samples and to date the approximate split between Alpine and Iberian ibex we took advantage of our heterochronous data set and conducted a serial birth‐death skyline analysis on 15,975 known sites. We ran the analysis on all Alpine ibex and Iberian ibex mitogenomes with beast v.2.5.2 (Bouckaert et al., 2014) and the beast package bdsky v.1.4.5 (Stadler et al., 2012 ). We used the default parameters except for the original node which was set to 10,000 years as the starting point for the birth‐death‐skyline‐serial prior (the original node must be older than the oldest sample, here 8600 years old). We used model averaging by applying the BEAST Model Test as site model (Bouckaert & Drummond, 2017), estimation of mutation rate under strict clock assumption and executed a MCMC chain of 1.5*107 samples with a burnin of 10% on. We sampled a total of 15,001 trees and used treeannotator v.1.10.4 to summarize the posterior estimates and HPD limits of the node ages of the target tree. The results were visualized with figtreev.1.4.3 (Rambaut, 2012) and edited in Adobe InDesign to improve readability.

Haplotype networks

To infer and visualize the haplotype diversity among recent and past Alpine ibex, we built a haplotype network including all Alpine ibex specimens (ancient, historic and recent, data set Capra_ibex). The tool vcf‐to‐tab (Chen, 2014) and the Perl script vcf_tab_to_fasta_alignment.pl (Chen, 2014) were used to convert the VCF format to FASTA. The FASTA sequences were then aligned using clustal‐omega v.1.2.4 (Sievers et al., 2011) which incorporates a clustalW algorithm. The resulting FASTA alignment was converted into PHYLIP using the Perl script Fasta2Phylip.pl by Hughes (2007). The r package TempNet (Prost & Anderson, 2011) was then used to build the haplotype network using the incorporated TCS algorithm. TCS calculates an absolute pairwise distance matrix of all haplotypes and connects the haplotypes according to the parsimony criterion to minimize mutation steps between haplotypes. The resulting haplotype network was edited for improved readability in Adobe InDesign.

Description of mitogenome

We calculated mitochondrial diversity statistics with the r package popgenome v.2.7.5 applied to data set Capra_all (Pfeifer et al., 2014). In PopGenome, by default, only sites genotyped in all individuals are retained. The gff3 file corresponding to the Alpine ibex mitogenome sequence NC_020623.1 was adapted for compatibility with the package PopGenome and was used with the command set.synonym to infer neutrality statistics. The neutrality statistics Tajima's D and the number of segregating sites were determined with neutrality.stats (GENOME.class,detail = TRUE). We calculated the nucleotide diversity for each species (in Alpine ibex per sample age group) with the function GENOME.class@Pi/GENOME.class@n.sites and the haplotype diversity within each population with GENOME.class@hap.diversity.within. Finally, we calculated nucleotide diversity within coding regions with the split_data_into_GFF_features and command GENOME.class@nucleotide‐diversity.within. Using the data set Capra_ibex, we furthermore constructed circos plots for the historic, ancient and recent Alpine ibex specimens to visualize the genetic variation along the mitogenome of the species with the r package circlize v.0.4.12 and a custom script programmed in r, published in Gu et al. (2014). Allele frequencies for all mitochondrial biallelic SNPs polymorphic among Alpine ibex (data set Capra_ibex) were calculated using the option ‐freq in vcftools (Danecek et al., 2011). SNP density among Alpine ibex (data set Capra_ibex) along the mitogenome was computed using the option ‐SNPdensity in vcftools (window size of 500 bp) and visualized with the r package ggplot2.

Skyline plot

We inferred the demographic trajectories for the female N e of Alpine ibex during the Holocene with a Bayesian Skyline approach applied to the mitogenomes (Drummond et al., 2005). This method facilitates Markov chain Monte Carlo (MCMC) sampling to infer a posterior distribution of effective population size through time by sampling directly from gene sequences. We used the data set Capra_ibex, including all Alpine ibex samples. We performed the analysis with beast2 v.2.5.2 (Bouckaert et al., 2014) and used model averaging by applying the BEAST Model Test as site model (Bouckaert & Drummond, 2017), enabled estimation of mutation rate under strict clock assumption and ran a MCMC chain of 1.5*107 samples and a burnin of 10%. Change in initial mutation rates yielded similar results (not shown). Tracer v.1.7.1 was used for visualization and Adobe InDesign to improve readability.

RESULTS

Demographic history of alpine ibex and related species

To explore past demographic fluctuations, we first took advantage of a published whole genome data set (Figure 1a) including Alpine ibex (N = 29, recent), Iberian ibex (N = 4), Bezoar (N = 6), Siberian ibex (N = 2), Nubian ibex (N = 2) and domestic goat (N = 16) (Alberto et al., 2018; Grossen et al., 2020). Average depth of coverage ranged from 6.9× to 46× (Table S1). The population dynamics inferred using the PSMC method suggested comparable trajectories for most species under study with two rises in N e estimates each followed by a N e contraction (Figure 1b, Figure S1). The first general contraction led to a local N e minimum between 300 and 400 kya for all species except for the Iberian ibex (at ~120 kya) and one Bezoar individual (~220 kya). One of the two Bezoar individuals showed a slightly later increase after the drop and both individuals suggested large N e estimates towards the more recent past. Both these individuals had relatively low whole genome sequencing coverage (6.9× and 10.8×, Table S1), which can lead to noise when estimating PSMC trajectories. Furthermore, PSMC analyses get less reliable towards the more recent past (Li & Durbin, 2011; Mather et al., 2020; Nadachowska‐Brzyska et al., 2016). While the observed trajectories were similar among its relatives, all Alpine ibex individuals stuck out with a nearly flat trajectories after their N e trajectory separated from Iberian ibex (approximately 200 kya, Figure 1b) with one shallow and one steeper contraction, suggesting that their demographic history considerably differed from the other species with relatively low N e estimates (Figure 1b). The steeper contraction (after ~25 kya, Figure 1b) coincides with the second contraction of Iberian ibex and the last glacial cycle 14–29 kya (Seguinot et al., 2018). During the peak of the last glacial cycle, the so‐called last glacial maximum (LGM, 19 to ~25 kya) nearly the entire Alpine ridge was covered by glaciers (Seguinot et al., 2018).

FIGURE 1

(a) Capra species distribution (except domestic goat) as stated by the IUCN and their phylogenetic relationship (maximum likelihood tree, RAxML) based on 80 whole mitogenomes. Alpine ibex (C. ibex,N = 44), Iberian ibex (C. pyrenaica,N = 4), Bezoar (C. aegagrus,N = 6), Siberian ibex (C. sibirica,N = 2), Markhor (C. falconeri,N = 1), Nubian ibex (C. nubiana,N = 2) and the domestic goat (C. hircus,N = 16) were included in the phylogeny. Sheep (Ovis sp., N = 5) were used as the outgroup. Nodes with bootstrap support lower than 100 are explicitly stated and all species branches were collapsed. N indicates the number of mitogenomes used, dotted lines indicate a cytonuclear discordance. (b) Pairwise sequentially Markovian coalescent approach (PSMC) analysis based on Alpine ibex specimens from the Gran Paradiso National Park (N = 6) and five related Capra species (N = 2 per each species). The PSMC was constructed over 29 autosomes, a generation time of 8 years and a mutation rate of 1.784e‐08 sites per generation were assumed. The last glacial cycle is depicted in blue and indicates the most recent cold glacial period. LGM, last glacial maximum (19–25 kya). (c) Demographic model used to estimate the split date between Alpine ibex and its sister species Iberian ibex using coalescent simulations in fastsimcoal2 (Excoffier et al. 2021 ) applied to the recent whole genome data

The long‐term effective population size of diploid individuals based on the coalescent simulations with fastsimcoal2 was estimated to 3874 (95% percentile CI: 2816–5086) for Alpine ibex and 112,454 (95% percentile CI: 110,941–113,759) and the split time between the two species to 28.2 kya (95% percentile CI: 20.3–37.3 kya) (Figure 1c).

Analysis of recent whole genome data allows exploring estimates for past effective population size (Li & Durbin, 2011). However, changes of N e estimates over time do not necessarily reflect actual species estimates of N e, but can be confounded by population structure (overestimation of N e) or restricted geographic distribution of direct ancestors (Mazet, Rodriguez & Chickhi, 2015; Wakeley & Aliacar, 2001). Furthermore, from recent data alone, it is difficult to estimate historic diversity, most of all if a species went through a strong bottleneck and lost large parts of its past diversity and hence signal (Mazet, Rodriguez & Grusea et al., 2015). To get a better understanding of the past demography of Alpine ibex, we therefore sequenced and analysed 15 ancient and historic mitogenomes.

Retracing the mitogenome diversity through a near extinction

To approximate the ancient genetic diversity of the Alpine ibex, we identified and collected samples which originated from several 1000 years before to within the last strong bottleneck during the 19th century (Figure 2a). We analysed whole mitogenomes using shotgun sequencing for seven ancient specimens (8601 ± 33 BP to 3302 ± 26 BP, six from caves, one from a glacier field) and eight historic specimens (1000 ± 50 CE to 1919 CE) originating from museums and archaeological excavations (Table S2, Figure 2a). We found that ancient samples performed better in terms of DNA preservation than historic samples. This is probably due to treatments used for preserving the historic specimens, which can negatively influence aDNA yield (McDonough et al., 2018). Post‐mortemdamage patterns were as expected for the respective sample age (Figure S2). The average endogenous DNA content was 62% (ranging from 28.9% to 92.6%, Table S2). For a broader view on past and present diversity, we complemented our data set with 65 additional published mitogenomes representing recent populations of six Capra species and domestic sheep (Figure 1a). The mean depth per individual ranged from 4.4× to 2449× for the historic and ancient samples and from 91× to 7093× for the recent samples (Table S2). Coverage along the mitogenome was >99.9% in all samples and 2220 sites were found to be biallelic. 45 segregating sites were found among all Alpine ibex (including ancient, historic and recent specimens) after filtering for genotype quality and missingness.

FIGURE 2

(a) Sample locations of ancient, historic and recent Alpine ibex specimens used for the study. Specimens which originated from the northern side of the main Alpine divide are coloured in dark tones, whereas specimens originating from the southern part are coloured in light tones. Specimens which could originate from both sides of the Alpine ridge are shown in both colour tones. The diameter of symbols indicates the sample size. Furthermore, triangles indicate ancient, diamond shapes historic and circles recent Alpine ibex specimens. (b). Serial birth‐death skyline tree based on 15,975 known sites for seven ancient (8601 ± 33 BP to 3302 ± 26 BP), eight historic (1000 ± 50 CE to 1919 CE), 29 recent Alpine ibex and four recent Iberian ibex mitogenomes. Branch labels indicate the posterior probabilities, node labels indicate the inferred split time and blocks around nodes the 95% HPD (highest posterior density). The last glacial cycle is indicated in blue. LGM: Last glacial maximum. Branch tips are colour coded as in 2a. *) indicates one individual with unknown origin. The analysis was performed in BEAST under strict clock assumption, an MCMC chain of 1.5*107 samples and a burnin of 10%. (c) Haplotype network including seven ancient, eight historic and 29 recent Alpine ibex samples. Little dots indicate mutational steps, size of pie chart indicates number of specimens representing the respective haplotype, colour coding as described above (see Figure S4 for further detail). (d) Mitogenome nucleotide diversity π per each Capra species with N ≥ 2 (Alpine ibex per sample age group) based on 15,986 known sites and shown as box plots. Siberian ibex are not included due to a cytonuclear discordance (see main text). The solid line indicates the median, the box spans from 25% to 75% of the interquartile ranges and upper and lower whisker spanning 1.5 * interquartile range

We first investigated the phylogenetic relationship between recent, historic and ancient Alpine ibex specimens and other Capra species by performing a maximum likelihood phylogenetic analysis (Figure 1a). All Alpine ibex samples built a well‐supported, monophyletic branch (lowest bootstrap value of 98), which was clearly distinct from the sister species Iberian ibex and all other Capra species in the tree (Figure 1a). The observed phylogenetic relationships among the other species (for instance Bezoar and domestic goat were more closely related to Alpine ibex than Nubian ibex) confirmed previous findings based on mitochondrial DNA by Pidancier et al. (2006), except for one of the two Siberian ibex clustering with Markhor (Figure 1a, Figure S3).

Next, to infer the mitochondrial split times between Iberian and Alpine ibex as well as between the major clades within species, we conducted a birth‐death serial skyline analysis based on the mitogenomes of Alpine and Iberian ibex only. We estimated the split time between Iberian and Alpine ibex to be approximately at 47 kya (95% HPD of 35–55 kya, Figure 2b). All recent Alpine ibex samples formed a monophyletic branch (posterior probability of 0.97) together with the historic southern samples and hence the known source population (light purple) and the oldest historic northern sample (Figure 2b). The most basal split separated two ancient and a historic northern sample from the other branches. This split was dated to shortly after the end of the last glacial cycle 14 kya ago (Figure 2b). Similarly, the most basal split among Iberian ibex was dated to 11 kya (Figure 2b).

For a more comprehensive description of the intra‐species diversity in Alpine ibex and to identify relationships in respect to their age and origin, we constructed a haplotype network of all Alpine ibex mitogenomes (Figure 2c, Figure S4, Table S2). We identified 12 distinct haplogroups (grouping haplotypes differing by up to two mutational steps) and 21 distinct haplotypes (Figure 2c). Most noticeable was the low haplotype diversity among the recent Alpine ibex populations with only two haplogroups (Figure 2c). The less frequent haplogroup was represented by six recent and two historic Alpine ibex individuals. Only one of the individuals originated from the northern part of the Alps, one from the main ridge (bicolour, Figure 2c) and the remaining samples from south of the Alps (recent and historic). The one from the main ridge is assumed to be the last Swiss Alpine ibex individual shot in 1846 (Urs Zimmermann, personal communication). The second major haplogroup was represented by a total of 26 individuals, both from south and north of the Alps and including three historic samples, all from the southern part of the Alps (Figure 2c). Within this haplogroup, 19 individuals carried the same haplotype (Figure 2C). Also two historic individuals (both dated to 1806 CE) from the Gran Paradiso National Park, source of all recent Alpine ibex populations, were assigned to the most abundant haplogroup. This may suggest that this haplotype was already common in this region during the time of the near extinction. All remaining historic and ancient specimens carried private haplotypes with up to 20 steps separating them from neighbouring haplotypes (Figure 2C).

When compared to related wild species, mitochondrial nucleotide diversities suggested lowest diversity in Alpine ibex, interestingly also for historic and ancient specimens (Figure 2d). Compared to ancient and historic samples as well as all other Capra species, we found in recent Alpine ibex the lowest number of segregating sites, the lowest haplotype diversity a nd the lowest nucleotide diversity (all measured across the entire mitogenome, Table 1, Figure 2d). Furthermore, recent Alpine ibex were monomorphic in eight out of 13 mitogenomic coding regions, which is in stark contrast to the ancient and historic specimens of the Alpine ibex and all other Capra species showing a maximum of two monomorphic coding regions (Table 1). A peak of diversity was found in the D‐loop control region with the highest SNP density among the ancient samples (Figure S5). The recent samples only showed 10 segregating sites, of which six within coding regions (Table 1). When comparing ancient and historic with recent diversity, a clear pattern of drift was observed, which led to the loss of most variants (Figure 3a). Variants with intermediate frequency (between 0.25 and 0.75) among historic samples were also more likely to still be segregating among recent samples (Figure 3b). Frequency among ancient samples did not necessarily predict frequency among recent samples (Figure 3a and Figure S6). This is not unexpected, because the ancient samples were mostly sampled north of the Alps, while a large proportion of the historic samples originate from the source of all recent Alpine ibex populations (Gran Paradiso National Park).

TABLE 1

Summary statistics calculated based on 15,986 mitochondrial sites per each Capra species with N > 2 (alpine ibex per sample age group). Siberian ibex are not included due to a cytonuclear discordance (see main text). Age group specifies the age of the samples: Ancient (8601 ± 33 BP to 3302 ± 26 BP), historic (1000 ± 50 CE to 1919 CE), and recent individuals

Age group.	Species	[N] Individuals	[N] Segregating sites	[N] Segregating sites in CDS	[N] Monomorphic CDS	[N] Haplotypes	Haplotype diversity	Nucleotide diversity	Nucleotide diversity in CDS
Ancient	Alpine ibex	7	34	22	2	6	0.86	8.09E‐04	9.11E‐04
Historic	Alpine ibex	8	37	26	2	6	0.75	6.61E‐04	7.23E‐04
Recent	Alpine ibex	29	10	6	8	5	0.17	1.71E‐04	1.54E‐04
Recent	Domestic goat	16	91	70	0	16	1	7.78E‐04	8.93E‐04
Recent	Iberian ibex	4	60	47	1	4	1	1.89E‐03	2.13E‐03
Recent	Nubian ibex	2	49	36	2	2	1	2.97E‐03	3.19E‐03
Recent	Bezoar	6	331	275	0	5	0.83	9.85E‐03	1.27E‐02

[N], number of individuals; [N] segregating sites, number of segregating sites in the mitogenome; [N] segregating sites in CDS, number of segregating sites in coding regions; [N] monomorphic CDS, number of monomorphic coding regions; [N] haplotypes, number of haplotypes; haplotype diversity, number of haplotypes/number of individuals; nucleotide diversity, nucleotide diversity across whole mitogenome; nucleotide diversity in CDS, nucleotide diversity in coding regions.

FIGURE 3

(a) Allele frequency of polymorphic sites along the mitogenome in ancient, historic and recent Alpine ibex. Coding regions are indicated in green. Vertical lines join the same site among sample age groups. Colours indicate frequency differences compared to recent Alpine ibex, with large differences indicated in red. Circle size represents absolute allele frequencies. The D‐loop is zoomed in for better readability (b) comparison of allele frequencies between recent and historic Alpine ibex. Circle size indicates the number of observations of a certain frequency combination. All sites segregating among Alpine ibex (ancient, historic or recent) were used for this graph. Hence, sites with a frequency of zero in both groups were segregating among ancient samples. A few alleles were only observed among recent (N = 29) but not among historic samples (N = 7). Rather than newly appeared mutations, these were more likely low‐frequency mutations in the past

We furthermore compared ancient, historic and recent Alpine ibex specimens by visualizing segregating sites along the mitogenome as a circos plot, further illustrating the low genetic diversity among recent specimens (Figure S7).

Stable demography before collapse during the last millennium

Historic records of Alpine ibex suggest a strong species bottleneck approximately 200 years ago (Brambilla et al., 2020; Grodinsky & Stüwe, 1987; Stüwe & Nievergelt, 1991). To determine the prebottlenecked demography of the species after the last glacial cycle (which would have been too close to the present for the PSMC to be reliable), we used the whole mitogenome data of all Alpine ibex individuals (seven ancient, eight historic and 29 recent mitogenomes) applying the Bayesian skyline plot approach (BSP) implemented in BEAST2 (Bouckaert et al., 2014). The estimates of past effective population size suggested no fluctuations between 5 and 13 kya with an average N e of approximately 35,000. A shallow decrease after ~6 kya to about 20,000 individuals at 2 kya was then followed by a steep drop to ~2500 individuals (Figure 4). Hence, the demography seems relatively stable after the end of the last glacial cycle around 14 kya (Seguinot et al., 2018; Spötl et al., 2021), in a time of significant environmental change and deglaciation (Figure 4). The first occurrence of presumed pastoralism above timberline in the Swiss Alps was dated to about ~6.9 kya (Schwörer et al., 2013) and there is evidence for hunting of Alpine ibex about 5 kya (Alpine ibex DNA found in the stomach of the iceman, Maixner et al., 2018). Evidence of significantly raising human activity along the European Alps can be found since 2.4 kya and intensified during the last 1000 years (Boxleitner et al., 2017, Figure 4). The large‐scale disappearance of the species is historically reported since the middle of the 16th century, which coincides with an elevated human population growth (Ziegler, 1963; Stüwe & Nievergelt, 1991; Head‐König, 2011) and the development of firearms (Westwood, 2005).

FIGURE 4

Bayesian skyline‐plot based on seven ancient (8601 ± 33 BP to 3302 ± 26 BP) and eight historic (1000 CE to 1919 CE) as well as 29 recent Alpine ibex mitogenomes with a total of 16,135 input sites. We used model averaging, incorporated a strict clock assumption, MCMC chain of 1.5*107 samples and used a burnin of 10%. Sample ages are indicated by ibex silhouettes along the time axis. Colour tones indicate geographic origin (as in Figure 2), lower x‐axis indicate the century. The invention of firearms and increased population growth are indicated schematically (see also main text)

DISCUSSION

In the light of a rapid global change and an ongoing mass extinction, it is crucial to understand the demographic and genetic consequences of bottlenecks in wild populations (Ceballos et al., 2020). However, in species which suffered a recent bottleneck, it becomes difficult to disentangle underlying causes of the low genetic variability as this requires the reconstruction of the prebottlenecked population (van der Valk et al., 2019; Casas‐Marce et al., 2017; Kirchman et al., 2020).

In the present study, we retraced patterns of genetic diversity through the near extinction of Alpine ibex and compared our findings with current population samples. We analysed 44 high‐quality mitogenomes spanning 8600 years and the current species range to get insight into the demographic and genetic history of the species. As expected from a near extinction event, we found a massive loss of haplotype diversity when comparing recent with prebottlenecked populations and identified 13 previously unknown haplotypes. However, already a few thousand years ago, estimates of mitogenomic diversity of Alpine ibex were lower than current estimates from related species. The analysis of published whole genome sequences revealed a relatively recent split from the sister species Iberian ibex and low long‐term population size in comparison to related species. This suggests that both long‐ and short‐term factors, environmental and human‐induced, have caused the genetic depletion of Alpine ibex. Our study shows how combining prebottlenecked with contemporary sampling provides a more comprehensive understanding of current patterns of diversity.

Phylogenetic analysis confirms previously discovered cytonuclear discordance

Our phylogenetic analysis across species based on whole mitogenomes confirmed a cytonuclear discordance previously discovered based on cytochrome b and Y‐chromosome sequences (Pidancier et al., 2006). Nuclear data places Bezoar and domestic goat basal to Siberian and Nubian ibex and the latter two species closer to the two sister species Iberian and Alpine ibex (Grossen et al., 2020 ; Pidancier et al., 2006). However, our phylogeny based on mitogenomes placed Bezoar and domestic goats next to Alpine and Iberian ibex (Figure 1a, Figure S3). Wild goat species can interbreed and hence Pidancier et al. (2006) suggested mitochondrial introgression among ancestral taxa to explain this cytonuclear discordance. Interestingly, our data set includes a Siberian ibex, which, based on whole genome data, was clearly grouped with a second Siberian ibex (Grossen et al., 2020), but here we found that its mitochondrial haplotype clustered with the Markhor (Figure 1a, Figure S3, Table S2). Hence, as a case example, it may have received mtDNA through introgression from Markhor, which is also plausible given the close spatial proximity of the two species distributions (Figure 1a). Cytonuclear discordance has been reported for a large number of organisms, including mammals and birds (Toews & Brelsford, 2012).

Alpine ibex show unique signals of demographic history compared to related species

The analysis of the recent Alpine ibex specimens using the PSMC method suggested a unique pattern of demographic history for Alpine ibex (Figure 1b, Figure S1). The estimates of effective population size of the related ibex species showed a relatively steep increase starting about 500,000 years ago to a local maximum between 170,000 (Siberian ibex) and 40,000 (Iberian ibex) before a further decrease. The absolute time estimates changed proportionally to the mutation rates used for the simulations (Figure S1). For instance, the local maximum of Iberian ibex was estimated to ~80,000 years ago when the mutation rate was halved (Figure S1). In contrast to all related species in this study, estimates for Alpine ibex only barely increased again after the drop of N e around ~250,000 years ago. Changes of coalescent‐based N e estimates over time can have different interpretations (Mather et al., 2020). An increase in N e estimates can for instance demonstrate increased population structure (Heller et al., 2013; Mather et al., 2020; Mazet, Rodriguez & Chikhiet al., 2015). A decrease in N e estimates, as observed in most species for the more recent past, is often simply the result of recent ancestors having lived in closer proximity and hence in more closely related populations (Wakeley & Aliacar, 2001). Thus, the rather flat N e trajectories of Alpine ibex over the last 200,000 years may be explained by long‐term local fidelity of the ancestor population. Alternatively, it could suggest relatively small species N e over a long timescale, in accordance with the small long‐term N e estimate of about 3874 since the split from the Iberian ibex 28.2 kya (the long‐term N e estimate for Iberian ibex was with 112,454 about 26 times higher). Similar declines of N e estimates after a lineage split were for instance found in Northern lions (de Manuel et al., 2020) and Amur leopards (Pečnerová et al., 2021) and were suggested to be the result of founding bottlenecks. The low long‐term N e estimates of Alpine ibex may also be explained by environmental factors. The Late Pleistocene was determined by large ice shields covering large parts of Northern Europe and the Alps (Seguinot et al., 2018). The last large glaciation period occurred 13.7 to 115 kya and the steepest drop in N e estimates of Alpine ibex (after ~25 kya) coincides with the last glacial cycle (14–29 kya), which had its peak (the last glacial maximum, LGM) between 19 and 25 kya (Seguinot et al., 2018). Considering the geographic distribution of Alpine ibex across the Alps (Figures 1a and Figure 2a) and their specific ecological needs (rocks and alpine meadows; Grignolio et al., 2003, 2007), it is plausible that the species was more strongly affected by these the last glacial cycle than related ibex species. Past climatic changes, as for instance the last glacial maximum, have been suggested to have influenced past demographies of a number of species including an increase in population size in Caribou (Rangifer tarandus) (Taylor et al., 2021), a marginal effect in the Pygmy hog (Porcula salvania) (Liu et al., 2020) and a decrease in population size in Musk ox (Ovibos moschatus) (Hansen et al., 2018; Campos et al., 2010) as well as Woolly mammoths (Mammuthus primigenius) (Palkopoulou et al., 2015). Notably, the N e trajectories of the endangered Iberian lynx show a comparable decline coinciding with the last glacial cycle (Abascal et al., 2016). Given these differing responses, there is an ongoing debate on the role of climatic changes in the Late Quaternary megafauna mass extinction, which identifies either environmental factors, very early human influence or the succession of both as the main drivers of the extinction events (Koch & Barnosky, 2006; Lord et al., 2020; Sandom et al., 2014; Seersholm et al., 2020; Stewart et al., 2021).

In accordance with prolonged times of relatively small population size of Alpine ibex is the observation of relatively low genetic diversity in the several 1000 years old ancient Alpine ibex in comparison to the diversity found in related ibex species (Figure 2d). Interestingly, the domestic goats showed similar mitogenomic diversity to the ancient Alpine ibex and higher diversity than recent Alpine ibex confirming previous results based on recent whole genome data (Grossen et al., 2020). In contrast to most domesticated species, domestic goats are assumed to have had genetic contributions of multiple wild populations over time (Alberto et al., 2018; Daly et al., 2018; Zheng et al., 2020) probably explaining overall higher than expected genetic diversity. However, while genome‐wide diversity in domestic goats was found to be highest among related ibex species (Grossen et al., 2020), mitogenomic diversity in domestic goats was (except for Alpine ibex) lower compared to all other ibex species in this study (Figure 2d). This discrepancy could be explained by presumed admixture events between goats from different sources (Daly et al., 2018) which counteracted the reduction of genome‐wide diversity usually induced by selective breeding (Bosse et al., 2018), while only one mitochondrial type became dominant worldwide (Daly et al., 2018).

Massive loss of mitochondrial diversity and increasing human activity

The PSMC analyses allowed investigating the demographic trajectories of Alpine ibex since the split from its sister species Iberian ibex to about 10–20 kya, when PSMC trajectories become less reliable (Nadachowska‐Brzyska et al., 2016). The analysis of the temporal mitogenomic data using the skyline approach then allowed more recent population size estimates from about 13 kya onwards (Figure 4), hence after the end of the last glacial cycle (Seguinot et al., 2018). The skyline plot indicated a stable effective population size until 5–6 kya, when a slow contraction starts followed by a drastic decline 1–2 kya (Figure 4). It is assumed that the Alpine ibex was bound to rocky areas and mountain meadows (Grignolio et al., 2007). Hence, the species was probably existing in a relatively narrow niche, which was defined by the upper timberline and the retreating ice shields of the glaciers. This could explain the rather shallow demographic trajectories of the species during the first few millennia after the last glacial cycle despite a fast changing environment. As violations of the scenario of a single, isolated and panmictic population in coalescent‐based demographic inferences can lead to spurious demographic signals, these results have to be interpreted cautiously (Heller et al., 2013). Nevertheless, the first contraction (5–6 kya) coincides with the first occurence of presumed pastoralism in the Swiss Alps above the timberline (c. 6.9 kya, Schwörer et al., 2013) and evidence of hunted Alpine ibex (Maixner et al., 2018). The human population grew since the early Bronze ages (about 3.5 kya) and started to work on alpine meadows and forests for agriculture up to high altitudes (Dietre et al., 2016; Röpke et al., 2011). During that time, Alpine ibex populations seem to have stayed relatively stable. The drastic drop starting 1–2 kya coincides with an increasing human activity in the Alps (Boxleitner et al., 2017; Stüwe & Nievergelt, 1991). The large‐scale disappearance of the species is historically reported since the middle of the 16th century, which coincides with elevated human population growth (Ziegler, 1963; Stüwe & Nievergelt, 1991; Head‐König, 2011) and the more frequent use of firearms (Westwood, 2005). Our N e estimates of the more recent past were based on historic samples from north and south of the Alps and therefore are probably a little overestimated. This, together with possible underestimations of census size in the past, may explain why our minimal N e estimates are higher than historical records reporting a reduction of the Alpine ibex census size down to about 100 individuals at the beginning of the 19th century (Brambilla et al., 2020; Grodinsky & Stüwe, 1987).

The diversity comparisons with related species and the PSMC analyses suggest that population size of Alpine ibex may have been relatively low already before the onset of intensive human activity. However, the main cause for the near extinction of Alpine ibex was most likely overhunting and possibly also an increased competition with domestic ungulates (Acevedo & Cassinello, 2009; Chririchella, Apollonio & Putman, ; Stüwe & Nievergelt, 1991). Alpine ibex are among a large number of species, in particular large mammals, which went through human‐induced bottlenecks up to extinctions (Ripple et al., 2016). Species assumed to have gone extinct at least partly due to early anthropogenic pressure include the cave bear (Gretzinger et al., 2019) and the Woolly mammoth (Palkopoulou et al., 2015; Fordham et al., 2021). Similar erosion of mitochondrial diversity when comparing contemporary to prebottlenecked diversity was for instance also found in the Iberian lynx attributed to intensive persecution followed by drift (Casas‐Marce et al., 2017). However, some of the erosion of Iberian lynx diversity already happened much earlier on (Abascal et al., 2016) similar to Alpine ibex and the Arctic fox (Dalén et al., 2005; Larsson et al., 2019).

The serial birth‐death skyline approach estimated the split dates between the mitochondrial lineages of Alpine and Iberian ibex to about 46 kya (Figure 2b), while the coalescent‐based analysis of the recent whole genome data suggested a more recent split around 28.2 kya (Figure 1c). This discrepancy is not necessarily surprising. In contrast to nuclear DNA, mitochondrial DNA is maternally inherited, does generally not recombine and is more strongly affected by drift than nuclear DNA (Fay & Wu, 1999). Introgression, which is known to have occurred repeatedly among wild goats (e.g., Pidancier et al., 2006; Grossen et al., 2014), would also lead to a younger split of nuclear compared to mitochondrial data. Both estimated dates are younger than a previously estimated split time of approximately 50–90 kya based on 780 bp of the cytochrome b gene (Ureña et al., 2018).

The phylogenetic analysis revealed several divergent lineages among Alpine ibex and the serial birth‐death skyline approach identified a well supported split within the Alpine ibex around 12 kya, after the end of the last glacial cycle. The basal split is in terms of timing and support comparable with the deepest split found among the Iberian ibex individuals (11 kya, Figure 2a). Interestingly, our study confirms previous findings of higher differentiation between two recent populations of the subspecies Capra pyrenaica hispanica than the differentiation to the other subspecies Capra pyrenaica victoriae (Groves & Grubb, 2011; Angelone‐Alasaad et al., 2017; Grossen et al., 2018). Among recent Alpine ibex, on the other hand, we find little differentiation. This is not surprising given the distinct demographic histories of the two species. Alpine ibex went through a very severe bottleneck (approximately 100 individuals) and are assumed to have only survived in one single population in Northern Italy (today the Gran Paradiso National Park; Stüwe & Nievergelt, 1991). The Iberian ibex went through a less severe bottleneck (approximately 1000 individuals) and survived in several isolated populations (Angelone‐Alasaad et al., 2017). As a consequence, the mitochondrial diversity remained larger in Iberian ibex than Alpine ibex confirming previous results based on genome‐wide diversity measures (Grossen et al., 2018, 2020a).

Interestingly, some lineages with deeper splits were still represented by historic Alpine ibex suggesting that some of the ancient diversity was still present very close to the near extinction. This is also in accordance with nucleotide and haplotype diversity of historic samples being intermediate between ancient and recent Alpine ibex (Figure 2d, Table 1). Such findings give hope for other bottlenecked species with long generation times, because a prompt recovery from a bottleneck may save considerable diversity. However, while the recent samples shared haplogroups with historic specimens, none of the ancient haplogroups were observed among the recent samples. Accordingly, historic allele frequencies were correlated with recent allele frequencies, but ancient allele frequencies were only very marginally correlated with the recent ones. This is probably explained both by space and time. First, it is not surprising that the recent samples are genetically more similar to the historic samples just because there was less time in between for the mitogenomes to evolve. Second, a large proportion of the historic samples originates from the Gran Paradiso population, the source of all recent populations. More data will be needed to investigate how strong past population structure was, if the Alps formed a barrier to gene flow and if Alpine ibex populations on both sides of the Alps were substantially differentiated from each other.

All the mitogenomic diversity observed among the recent samples was represented by only two haplogroups with nine mutational steps in between. Although 29 specimens is not a very large sample size, the sampling was explicitly chosen to represent the current species diversity (Biebach & Keller, 2009; Grossen et al., 2018, 2020; Kessler et al., 2020). Hence, it is unlikely that several common recent haplogroups remain undetected. As expected from the species history, individuals from Gran Paradiso (recent and historic) were found in either haplogroup. All recent individuals sampled North of the Alps (except for one) belonged to the most abundant haplogroup. The second haplogroup was represented by several recent and historic individuals, in particular from the population Alpi Marittime in Italy, a reintroduced Alpine ibex population occurring at the southern edge of the species distribution (Figures 1a and Figure 2a). This population was reintroduced in the early 20th century, based on only about six founders. Accordingly, it shows high inbreeding and is highly divergent from all other existent Alpine ibex populations (Grossen et al., 2020; Kessler et al., 2020).

The most striking loss of diversity was observed in the D‐loop (Figure 3). The D‐loop is the regulatory region of the mitochondrial DNA and responsible for its replication and transcription (Nicholls & Minczuk, 2014). It contains two hypervariable regions (HV‐I, HV‐II), which in humans have a 100 to 200 times higher mutation rate than the nucleus (Sharawat et al., 2010). Due to its high substitution rate, the D‐loop can help to resolve differences between closely related individuals (Kundu & Ghosh, 2015). The substantial differences in diversity in this region also underlines the severity of the most recent bottleneck which erased much of the rapidly evolving genetic diversity of the D‐loop.

Conclusions

The complementary analysis of ancient, historic and recent data allowed us to investigate the demographic history of a species from its emergence through its near extinction to the very recent past. The demographic modelling using whole genome data from recent populations and diversity comparisons with related species, suggest that Alpine ibex population size was reduced over prolonged times soon after the species' emergence and previous to a large‐scale anthropogenic land‐use in the alpine region. Hence, Alpine ibex demography was probably affected by long‐term environmental processes such as glaciation. While the PSMC analysis of the whole genome data of recent individuals give insight into the early demographic history, sufficient resolution in the more recent demography is lacking. Here, access to ancient and historic mitogenome data enabled the investigation of the more recent demographic history of the species. We found massive loss of mitogenome diversity and identified an increase in anthropogenic pressure during the last 1–2 millennia as the main cause of the low genetic diversity of contemporary Alpine ibex populations. Our study underlines the value of a combined approach of ancient and historic mitogenomes, demographic modelling based on contemporary and related species data to understand past population fluctuations and their consequences on contemporary patterns of genetic diversity. Such complementary studies of severely bottlenecked and successfully recovered species are crucial to understand the scale and the genetic consequences of small population size and population fragmentation in the wild giving important insights for species conservation.

AUTHOR CONTRIBUTIONS

Giada Ferrari and Christine Grossen conceived the project. Mathieu Robin acquired all samples. Mathieu Robin, Giada Ferrari, Johanna von Seth and Gülfirde Akgül carried out the sample and library preparation. Love Dalén and Verena J. Schuenemann supported the project by enabling access to laboratory, consumables and sequencing facilities and giving advice. Mathieu Robin, Xenia Münger and Christine Grossen conducted all bioinformatic analysis. Mathieu Robin and Christine Grossen wrote the manuscript with input from all coauthors.

CONFLICT OF INTEREST

The authors declare no competing interests.

Appendix S1

Click here for additional data file.