Population structure and genetic diversity of non-native aoudad populations

The aoudad (Ammotragus lervia Pallas 1777) is an ungulate species, native to the mountain ranges of North Africa. In the second half of the twentieth century, it was successfully introduced in some European countries, mainly for hunting purposes, i.e. in Croatia, the Czech Republic, Italy, and Spain. We used neutral genetic markers, the mitochondrial DNA control region sequence and microsatellite loci, to characterize and compare genetic diversity and spatial pattern of genetic structure on different timeframes among all European aoudad populations. Four distinct control region haplotypes found in European aoudad populations indicate that the aoudad has been introduced in Europe from multiple genetic sources, with the population in the Sierra Espuña as the only population in which more than one haplotype was detected. The number of detected microsatellite alleles within all populations (< 3.61) and mean proportion of shared alleles within all analysed populations (< 0.55) indicates relatively low genetic variability, as expected for new populations funded by a small number of individuals. In STRUCTURE results with K = 2–4, Croatian and Czech populations cluster in the same genetic cluster, indicating joined origin. Among three populations from Spain, Almeria population shows as genetically distinct from others in results, while other Spanish populations diverge at K = 4. Maintenance of genetic diversity should be included in the management of populations to sustain their viability, specially for small Czech population with high proportion of shared alleles (0.85) and Croatian population that had the smallest estimated effective population size (Ne = 5.4).

www.nature.com/scientificreports/ of aoudad in its native range are declining 11,12 , classifying this species as Vulnerable in the IUCN Red List 13 . The aoudad is a generalist herbivore and extremely plastic in its utilization of available food resources 14 . It is a polygynous species with high reproductive rate. Large and strong individuals have high fitness and reproductive success 15 , increasing their potential to colonise different localities whenever conditions are appropriate 11 . In Croatia, during 2002 a number of aoudads of unknown origin were illegally released in the southern Dinaric region (Mosor Mountain). Current data obtained by GPS tracking of aoudad movements showed that the population is limited to the Mosor Mountain range and it is highly unlikely that it will expand its range beyond that area 16 . The population in the Czech Republic was established following the escape of several individuals from the Plzeň Zoo in 1976, and it consisted of a few dozen animals 17 , while today this population is believed to be extinct (Cupic S., pers. comm.). In Italy, one small population of about 20 individuals is known to be present at the Beigua Regional Park in the province of Savona (Liguria) since 2007 15 . Previously, another population was present in the province of Varese (Lombardy), where a group of six individuals (one male and five females) escaped from a private enclosure in 1993 and established a small breeding population that reached the maximum size of about 20 individuals. This population has been eradicated by the personnel of the Province of Varese accordingly to a plan approved by the Italian National Institute for Environmental Protection and Research (ISPRA) since 2005 (Martinoli A., pers. comm.). In Sierra Espuña, Spain, 36 aoudads, from zoos in Germany and Morocco, were intentionally introduced by the regional administrations in 1970. This founding population reproduced with great success and naturally dispersed very rapidly from their release area to nearby montane enclaves. Its current population, estimated at around 2000 individuals, is still expanding 11,18,19 . In 1972 sixteen animals from Sierra Espuña were released into La Palma island, Canary Islands, and established successfully since then 10,18 .
So far, there is no evidence of the negative impact of the aoudad on host ecosystems in mainland Spain, except for La Palma island, where they have critically affected the survival and diversity of native, endemic flora and caused high levels of soil erosion 11,20 . There is no information on environmental impacts for the other nonnative aoudad populations 15 .
When the introduction of species into new areas is human-mediated, the new populations are often founded by only a few individuals that are completely isolated from the source populations. In the case of the European aoudad populations, sudden and substantial reduction in effective population size during first introductions and lack of gene flow into the established populations are assumed to have led to the loss of genetic variation through genetic drift. Severe reductions in genetic diversity may limit viability and adaptive potential of introduced populations under environmental change 5 , as adaptation in such species occurs mainly through selection on pre-existing genetic variation 21 . Since scientific literature about the genetic structure of European aoudad populations is scarce, it is thus important to study genetic diversity and structure of those populations to assess their sustainability. Unfortunately, no comprehensive genetic analyzes have yet been carried out on original aoudad populations from their native range, except for a few studies of specific populations (e.g., Derouiche et al. 22 ). In addition, insight into the genetic structure and origin of European aoudad populations might contribute to the ex-situ conservation of the species that is threatened in its native range.
Here we report the first attempt to characterize the genetic diversity and population genetic structure of all non-native European aoudad populations. We used neutral genetic markers, mitochondrial DNA (mtDNA) control region sequence and microsatellite loci, to characterize and compare genetic diversity and spatial pattern of genetic structure on different timeframes among study populations. While variability of maternally inherited and more conserved mtDNA control region reflect maternal lineages present in the founder individuals, microsatellites, due to their codominant nature and high variability can reflect more recent events that shaped current genetic structure 23 . We aimed to understand how current levels of genetic diversity and structuring vary among European aoudad populations that differ with regard to the time and source of introduction.
Based on samples collected from all known locations in Europe, the objectives of this study were to characterize the patterns of neutral genetic structure of recently established non native populations of aoudad and gain insights into the number of maternal lineages of these populations.

Results
Mitochondrial sequence analysis. Sequencing of the mtDNA control region was successfully performed on 69 samples from five different populations (Croatia n = 32, Czech Republic n = 1, Sierra Espuña n = 14, Almeria captivity n = 7, La Palma n = 15). None of the samples collected in Italy were successfully sequenced, due to the low quality of DNA. Low quality of the DNA isolated from faecal samples is known to dramatically reduce the genotyping success of such samples 24 . In the full sample, four mtDNA haplotypes (GenBank accession numbers MW349820-MW349823, Fig. 1) were detected.
Haplotype Amle01 was the most frequent (70%) and the only one found in Croatian, Czech and La Palma samples. Samples from Almeria were monomorphic for Amle02 haplotype, therefore haplotype and nucleotide diversity was zero in these populations. Samples from Sierra Espuña comprised of three haplotypes (Amle01, 03 and 04) with haplotype diversity of 0.69 and nucleotide diversity of 0.06. Observed four haplotypes were defined by 59 variable sites of which 11 were parsimony informative. The least number of mutation (12) was present between haplotypes Amle02 and Amle03, while the haplotype Amle04 stood out as remarkably divergent from all others, with the maximum of 55 variable sites separating Amle04 and Amle02 (Fig. 2). Presence of four distinct mtDNA haplotypes in our samples indicates that at least four maternal mtDNA lineages are present in current aoudad populations in Europe.
Microsatellite genotyping. We succesfully genotyped 85 individuals. Eighty three of them were genotyped at all 15 microsatellite markers, while two samples were genotyped at fourteen and thirteen markers, www.nature.com/scientificreports/ respectively. In the Croatian population, two multilocus genotypes were shared between four individuals. The unbiased probability of identity, or the probability that two individuals randomly drawn from the population had the same multilocus genotype, was 1.47 × 10 −7 (Supplementary Table S1). Therefore, shared genotypes were included only once. In addition, one sample containing more than three null alleles was excluded from the dataset. Italian samples had 25% of successfully genotyped markers, probably due to low amount of DNA obtained from fecal samples and were thus excluded from the following analyses. Accordingly, a final data set consisted of 82 unique multilocus genotypes (Croatia n = 31, Czech Republic n = 4, Sierra Espuña n = 17, Almeria captivity n = 15, La Palma n = 15).

Hardy-Weinberg equilibrium and within-population genetic diversity. The fifteen microsatellite
loci yielded a total of 76 alleles, ranging from 2 (BM302, MAF70 and TGLA073) to 9 (INRA040), with a mean number of alleles per locus of 5.07 (Supplementary Table S2). The PIC values ranged from 0.32 to 0.82, with an average of 0.57 (Supplementary Table S2). No evidence of scoring errors due to stuttering or large allele dropout was found in the whole data set. According to MICRO-CHECKER, the presence of null alleles was suggested in one locus/population combination (BM143 in Croatia) with the frequency of null alleles equal to 0.096 (Supplementary Table S3). However, Croatian population showed no departure from Hardy-Weinberg equilibrium (HWE) at this locus. The allelic richness varied from 1.93 (Almeria, Spain) to 3.67 (Sierra Espuña,  www.nature.com/scientificreports/ Spain) ( Table 1). Private alleles were observed in each population, with a total of 24 private alleles detected in four populations. The highest number of private alleles (9) was observed in captive population from Almeria.
Observed heterozygosity values were between 0.358 in Almeria and 0.541 in Sierra Espuña, while unbiased expected heterozygosity ranged from 0.337 in Almeria to 0.564 in Sierra Espuña (Table 1). Linkage disequilibrium was observed for: 14 pairs of loci in Croatia, four pairs in Sierra Espuña, and for one pair of loci in La Palma and Almeria. After applying FDR corrections, linkage disequilibrium was significant only for MM12/SR-CSRP12 and TGLA073/SR-CSRP12 in Croatia. The multilocus value of the f estimator ranged between − 0.029 (Almeria) to 0.070 (Sierra Espuña) with a mean positive value of 0.302 ± 0.057.
Estimated effective population size was the smallest in Croatian population reflecting the small number of founding individuals, and the largest in Sierra Espuña.
Mean proportion of shared alleles between all individuals within each population was: 0.56 in Sierra Espuña, 0.57 in La Palma, 0.66 in Croatia, 0.75 in Almeria and 0.85 in the Czech population.
Genetic differentiation and structure. Global F ST value was 0.296 (95% CI = 0.249-0.347). The lowest pairwise F ST value was observed between populations from Czech and Sierra Espuña (F ST = 0.109), while the highest value was found between Croatian and Almeria populations (F ST = 0.459) ( Table 2). Global and all pairwise F ST values were significantly different from zero (P < 0.01).
Among STRU CTU RE runs, the highest ΔK value was observed for K = 3, followed by that for K = 4 (Fig. 3). Captive population from Almeria was unequivocally discriminated from the individuals belonging to freeranging populations at each presented K value. After accounting for this major split, our results suggest the presence of two clusters in the four free-ranging populations surveyed. The first cluster is composed of Croatian and Czech populations, indicating that they belong to the same ancestral population. The second cluster includes Sierra Espuña and La Palma populations, while at K = 4, they appear as separate.

Discussion
While this research aimed to present the genetic structure of all aoudad populations in Europe, we were not able to obtain any results for the Italian population due to the very low quality of DNA sampled from fecal samples.
Four distinct haplotypes found in European aoudad populations indicate that aoudads have been introduced in Europe from at least four maternal lineages. Sierra Espuña is the only population in which more than one haplotype was detected. The presence of multiple maternal lineages in this population, along with largest estimate of effective population size of all sampled populations, indicates that it was established by females with at least three different mtDNA haplotypes. Higher genetic diversity of Sierra Espuña population than of all other European populations is further supported by its lowest mean proportion of shared alleles (0.56) between individuals from this population and the largest estimated effective population size of 36.
Sierra Espuña population is known to have originated from individuals from the Frankfurt Zoo in Germany and the Ain Sebad Zoo in Casablanca, Morocco 10,25 . There is no available data on subspecies/lineages of aoudads used in introductions to Sierra Espuña (or in preceding introductions to Frankfurt and Ain Sebad Zoos), nor published mtDNA sequences from original populations that would allow for identification of subspecies/lineages Table 1. Genetic diversity assessed from fifteen microsatellite markers in four European aoudad populations. n = number of successfully genotyped individuals per sampling population, N av = mean number of alleles per locus, N ar and N par = rarefied allelic and private allelic richness (smallest sample size = 28 alleles, as explained in materials and methods section), N pr = total number of private alleles, H o = observed and H e = expected heterozygosity, f = estimator of the inbreeding coefficient (none of the values were statistically different from 0, with P < 0.05).  www.nature.com/scientificreports/ present in this population, but variability of detected haplotypes might reflect introductions from multiple distinct lineages. This is especially valid for haplotype Amle04, which differs by at least 46 mutations from other haplotypes detected in this population, indicating its origin from evolutionary distinct units. Reported numbers of mutations between haplotypes from same subspecies of ungulates are usually lower. For example, in the study of chamois (Rupicapra sp.) 26 , native European mountain ungulate, D-loop haplotypes from the same subspecies were up to 15 mutations apart. In its native range in North Africa, six aoudad subspecies have been described, based on their distribution and morphological differences in coat color and horns 27 . Derouiche et al. 22 aimed to further evaluate the validity of the subspecies and determine the geographical limits of the valid ones using mitochondrial cytochrome b gene. Based on the haplotypes obtained, their results indicated a deep Mediterranean-Saharan genetic break in the species, suggesting the presence of two highly distinct evolutionary lineages. All analyzed aoudad populations in Europe were only recently introduced, with relatively small number of founding individuals, so it is expected that their genetic diversity is low in comparison to the viable free mating populations of native ungulates. This is confirmed by the number of detected alleles within all populations (N ar was less than 3.61 in all analysed populations) and mean proportions of shared alleles within all populations (> 0.55) which confirm their poor genetic diversity. The theoretical expectations for the loss of allelic variation after bottlenecks are described by Nei et al. 28 who showed that critical factors determining this loss are the size www.nature.com/scientificreports/ of founding population and the growth rate of a newly established population. They indicated that the average number of alleles per locus will be more sensitive to founding population size than heterozygosity. This is illustrated in most of sampled aoudad populations, in which allelic richness is substantially low (Table 1), while deficit of heterozygotes is not significant. Similarly low levels of genetic variability were previously reported for native European mountain ungulates (e.g. Šprem and Buzan 29 , for chamois), where it was attributed to isolation and poor management practices. Largest number of detected alleles and the highest observed heterozygosity in Sierra Espuña population is another indication of the highest genetic variability in this population, supporting the theory of multiple origins. In all presented STRU CTU RE results, Croatian and Czech populations cluster in the same genetic cluster, indicating joined origin. Among three populations from Spain, Almeria population shows as genetically distinct from others in results from K = 2 to K = 4, while Sierra Espuña and La Palma diverge only at K = 4. It is worth pointing out that the population from Almeria originated from a founding couple captured in the Atlantic Sahara (south of Morocco) back in 1975; so that it presumably belongs to the Saharan aoudad subspecies, while the other ones analyzed in our study might be of admixed origin, based on the fact that it is a population formed by individuals coming from the Casablanca Zoo in Morocco (probably belonging to the Atlas subspecies, the most abundant subspecies in that country), and from the Frankfurt Zoo, of unknown origin. In addition, Almeria population had the second highest mean proportion of shared alleles (0.75), lower only than the small Czech population.
Management of European aoudad populations varies between countries and the views change with new insights of their coexistence with native species and ecosystems 30 . Implications of detected low genetic diversity on management of these populations can be discussed only after management goals are clearly defined. Commonly accepted practice is to increase genetic variability of introduced species by introducing more individuals with different genetic background. The effects of such practices were reviewed by Dlugosch and Parker 31 who quantitatively summarized the genetic diversity data available at the time for 80 introduced or invasive species of animals, plants and fungi. Their review shows that multiple introductions had only small positive effects on diversity on average, disapproving the argument that multiple introductions have been critical in providing genetic rescue from severe and deleterious founder effects in most cases. According to this review, increases in genetic diversity caused by multiple introductions and/or gene flow occur after many decades, during which time the range expansion of the founder population is successful. Later introductions of genetically distinct individuals can increase genetic diversity through admixture, and the strength of this effect will depend on the reproductive rate of the individuals from different "waves" of introductions.
In Croatia, 34 tissue samples were collected during regular hunting. In the Czech Republic, four blood samples were collected from live aoudad because they are protected and not hunted 17 . Since the aoudad population in Italy is small and we were not able to collect any tissue samples, we collected seven fresh faeces samples. In Spain 47 tissue samples were collected from three locations: Sierra Espuña n = 17, La Palma n = 15, during regular hunting, and from captivity in Almeria n = 15. Geographic locations of the sampled populations are presented in Fig. 1.
Sampling Mitochondrial control region was amplified using tRNAPro and tRNAPhe primers 34 , following the PCR protocol given in Mereu et al. 34 . The PCR products were sequenced in forward direction using an ABI 3730 automated DNA sequencer (Applied Biosystems). Electropherograms were visually inspected using Applied Biosystems SEQSCAPE software and the sequences were trimmed to 767 bp which could be unequivocally called.  45 . Allelic richness [N ar (g)] and private allelic richness [N par (g)] were estimated using a rarefaction procedure implemented in HP-RARE 46 , where g represents the minimum number of alleles observed at a locus in one of the populations (i.e., twice the number of genotypes). The minimum number of genes in analysed populations was 28 (since locus locus INRA005 failed to aplify for one individual in the Almeria population), so this was used as a basis for rarefraction.
We estimated effective population size (Ne) for each population using the linkage disequilibrium method 47 , as implemented in NEESTIMATOR 2.01 48 . As suggested by Waples and Do 49 , we excluded alleles with frequencies below 0.02 in order to avoid bias caused by rare alleles.
Finally, to estimate the genetic similarity of individuals within populations, we estimated the mean proportion of shared alleles between all individuals within each population using custom script in Microsoft Excel 2016 (Microsoft Corporation, Redmond, WA, USA, https:// www. micro soft. com/ en-ca/ micro soft-365/ excel).
Genetic differentiation and structure. We estimated overall and pairwise genetic differentiation with the θ estimator of F ST 44 and determined their statistical significance by 1000 permutations in FSTAT. Finally, to infer the number of ancestral populations among non-native aoudad populations, we identified genetic structure using a model-based clustering method implemented in STRU CTU RE v. 2.3.4 50 . Ten runs per cluster (K), with K ranging from 1 to 7, were carried out with 10 6 iterations after a burn-in period of 10 5 iterations. We considered no admixture model and uncorrelated allele frequencies. The most likely number of clusters (K) was estimated by estimation of the ΔK statistic that reflects the rate of change of log-likelihood values between sets of runs with successive K values, using STRU CTU RE HARVESTER 51 . To account for label switching between results of different runs with same K, results from 10 runs with selected K were combined using CLUMPP v.1.1.2 52 . Results were displayed graphically using DISTRUCT v.1.1 53

Data availability
Data generated and analysed during this study are included in the Supplementary Information file and are available in the GenBank repository (Accession Nos. MW349820-MW349823).