An Assessment of Applicability of SNP Chip Developed for Domestic Goats in Genetic Studies of Caucasian Tur ( caucasica

: Caucasian tur ( Capra caucasica ) is native to Greater Caucasus Mountain Chain from Azerbaijan and Georgia in the East to Krasnodar region of Russia in the West. This species is divided into two subspecies (by some authors into species)—East-Caucasian tur and West-Caucasian tur and a subpopulation referred to as Mid-Caucasian tur. Up to date most of the genetic studies of Caucasian tur are based on mitochondrial DNA sequences and comprehensive investigation based on nuclear DNA is required for clariﬁcation of its genetic diversity and population structure. In our work, we assessed the applicability of Illumina Goat SNP50 BeadChip for genetic studies of Caucasian tur. Total of 15 specimens of Capra caucasica including East-Caucasian tur from Dagestan (E_TUR, n = 5), West-Caucasian tur from Karachay-Cherkessia (W_TUR, n = 5), and Mid-Caucasian tur from Kabardino-Balkaria (M_TUR, n = 5) were genotyped. After quality control, 5544 polymorphic loci, which were distributed all over 29 autosomes, were detected. The lowest number of SNPs was found on the 25th chromosome—68, and the highest on the 1st chromosome—348. It was shown that all the three groups of Caucasian tur clustered separately. A total of 2061 SNPs were common for all the populations, 594 were found only in W_TUR, 689 in E_TUR, and 530 in M_TUR. Individual heterozygosity ranged from 0.273 to 0.282 in W_TUR, from 0.217 to 0.253 in E_TUR, and from 0.255 to 0.283 in M_TUR. A clinal pattern of genetic variation was revealed. It was suggested to consider Caucasian tur a single species with several ecotypes. Thus, in our study we demonstrated that the Illumina Goat SNP50 BeadChip developed for domestic goats can be used as a useful tool for genetic studies of Caucasian tur.


Introduction
Caucasian tur (Capra caucasica) is an ungulate species from the subfamily Caprinae ( Figure 1). Its habitat is limited to Greater Caucasus Mountain Chain from Azerbaijan and Georgia in the East to Krasnodar region of Russia in the West (Figure 2). Its range is one of the smallest habitats of all ungulates-around 770 km in length and 80 km in width [1]. Taxonomically, the species is divided into two subspecies: West-Caucasian tur (C. c. severtzovi) and East-Caucasian tur (C. c. cylindricornis). Some authors even considered these populations as different species [2][3][4][5]. The areas of these subspecies are united, and  in the territory between them a hybrid form, which is referred to as a Mid-Caucasian tur, is distinguished [6].  The International Union for Conservation of Nature (IUCN) has listed East-Caucasian tur with the total numbers of about 31,000-32,000 animals as near threatened [7] and West-Caucasian tur with census size only around 5000 animals as endangered [8]. Of particular concern is an increasing anthropogenic pressure (poaching, habitat destruction, grazing competition with livestock, etc.,), resulting in migrations to higher mountains and adaptation to nocturnal behavior in some populations.
Little is known about the genetic diversity and evolutionary history of this animal. Only several studies based on Y-chromosome SRY gene and mitochondrial cytochrome b and control region were conducted [9][10][11][12]. To clarify the taxonomy of Caucasian tur and investigate the diversity of its populations molecular genetic studies based on nuclear DNA are needed.
Currently, one of the most effective technologies for molecular genetic research is the genome-wide analysis of single-nucleotide polymorphisms (SNPs). Recent advances in the development of high-throughput genotyping platforms (SNP chips) have turned SNPs into a powerful tool for genetic studies of domestic animals [13,14]. Although such chips are not available for most of the wild animals, the use of the SNP chips developed for their domestic relatives were found to be suitable.
The first SNP chip for domestic goats-Goat SNP50 BeadChip, containing more than 50,000 markers was developed by the International Goat Genome Consortium in 2011 [22]. In terms of wild species, Goat SNP50 BeadChip, so far was only used for genotyping of Bezoar goat (Capra aegagrus) which is considered the ancestor of domestic goats [23].
It should be noted that for non-model species the genotyping call-rate decreases approximately 1.5% per each million years of divergence time between species and the number of polymorphic sites decline exponentially leveling off after about 5 million years of divergence [20]. Therefore, this approach is not suitable for some analyses, i.e., based on linkage disequilibrium (LD), but meanwhile could be successfully used for investigation of phylogenetic relationships, genetic diversity, admixture, and introgression.
According to the web resource TimeTree (http://www.timetree.org, accessed on 15 May 2021) [24] the estimated median time of Caucasian tur and domestic goat divergence is 1.36 million years ago, which is much less than for the majority of other non-model animals. It gives a great opportunity to investigate phylogeny and genetic characteristics of Capra caucasica.
In this regard, the aim of our study was to assess the applicability of Illumina Goat SNP50 Bead-Chip, developed for domestic goats, for genetic studies of Caucasian tur (Capra caucasica).

Materials and Methods
A total of 15 specimens of Capra caucasica, including East-Caucasian tur from Dagestan (E_TUR, n = 5), West-Caucasian tur from Karachay-Cherkessia (W_TUR, n = 5), and Mid-Caucasian tur from Kabardino-Balkaria (M_TUR, n = 5) were selected for this study. The sampling sites are presented in Figure 2. All the samples were taken from trophy hunters of the Mountain Hunters Club (www.kgo-club.ru), who were licensed to hunt Caucasian tur. DNA extraction was carried out using Nexttec columns (Nexttec Biotechnology GmbH, Leverkusen, Germany) in accordance with the manufacturer's recommendations. SNP genotyping was performed using the Illumina Goat SNP50 BeadChip containing 53,347 SNPs. SNP quality filtering was performed in PLINK 1.9 [25]. SNPs with unknown positions, located on sex chromosomes as well as those that were genotyped in less than 90% of individuals (− −geno 0.1), with a minor allele frequency (MAF) < 1% (− −maf 0.01) and with deviations from Hardy-Weinberg equilibrium (− −hwe 1e−6) were removed.
For the population genetic analyses (PCA, Admixture, f3 statistics, etc.,) SNPs in linkage disequilibrium (− −indep-pairwise 50 5 0.5) were pruned. The observed (H O ) and unbiased expected (H E ) heterozygosity [26], inbreeding coefficient (F IS ), and rarified allelic richness (A R ) were calculated in the R package "diveRsity" [27]. Pairwise F ST genetic distances [28] were calculated in the R package "StAMPP" [29]. Principle component analysis was performed with PLINK 1.9 (− −pca 4) and visualized in the R package "ggplot2" [30]. An individual tree based on the pairwise identity-by-state (IBS) distance matrix (− −distance 1-ibs) was constructed using the Neighbor-Net algorithm implemented in SplitsTree 4.14.6 [31]. To estimate and visualize the distribution of heterozygosity at the individual level, multilocus heterozygosity (MLH) was calculated in the R package "inbreedR" [32]. Venn diagram was constructed in the R package "VennDiagram" [33]. Cluster analysis was performed in Admixture 1.3 software [34], and visualized in the R package "pophelper" [35]. f3 statistics [36] was computed using ADMIXTOOLS [37], as implemented in the R package "admixr" [38]. The map with sampling sites was created using R packages "maps" [39] and "ggplot2". As long as the Illumina Goat SNP50 BeadChip was developed for domestic goats (Capra hircus), the positions on chromosomes in this paper correspond to the reference genome of a domestic goat.

Results
From the initial set of 53,347 SNPs available in the Illumina Goat SNP50 BeadChip after quality control procedures (− −geno 0.1, − −maf 0.01) 5544 variants (10.4%) passed all the filters while 7002 variants were removed due to missing genotype data and 37,407 variants were removed due to minor allele threshold. No SNPs were removed due to Hardy-Weinberg exact test (Table S2) Mean minor allele frequency (MAF) computed for all the populations was 0.214 ± 0.002. About 2474 SNPs were highly informative with MAF > 0.3 and 1895 SNPs had minor allele frequency less than 0.1 ( Figure S2).
While 2061 polymorphic SNPs were common for all the Caucasian tur groups, the number of unique SNPs for W_TUR, E_TUR and M_TUR was 594, 689, and 530, Mean minor allele frequency (MAF) computed for all the populations was 0.214 ± 0.002. About 2474 SNPs were highly informative with MAF > 0.3 and 1895 SNPs had minor allele frequency less than 0.1 ( Figure S2).
While 2061 polymorphic SNPs were common for all the Caucasian tur groups, the number of unique SNPs for W_TUR, E_TUR and M_TUR was 594, 689, and 530, respectively. The most shared SNPs between two groups was found in W_TUR-M_TUR pair-869, and the least number of shared loci was observed in W_TUR-E_TUR pair-420 ( Figure 4). Mean minor allele frequency (MAF) computed for all the populations was 0.214 ± 0.002. About 2474 SNPs were highly informative with MAF > 0.3 and 1895 SNPs had minor allele frequency less than 0.1 ( Figure S2).
While 2061 polymorphic SNPs were common for all the Caucasian tur groups, the number of unique SNPs for W_TUR, E_TUR and M_TUR was 594, 689, and 530, respectively. The most shared SNPs between two groups was found in W_TUR-M_TUR pair-869, and the least number of shared loci was observed in W_TUR-E_TUR pair-420 ( Figure 4). For all the next analyses LD pruning was performed, after which 3885 SNPs were selected.
The results of the principal component analysis (PCA) revealed that all the studied groups of Caucasian tur clustered separately ( Figure 5).  For all the next analyses LD pruning was performed, after which 3885 SNPs were selected.
The results of the principal component analysis (PCA) revealed that all the studied groups of Caucasian tur clustered separately ( Figure 5). The first component (PC1), which explained 2.48% of genetic variability, divided Caucasian tur into two groups. The first group consisted of E_TUR (PC1 < 0) and the second one included W_TUR and M_TUR (PC1 > 0).
The level of genetic differentiation was estimated with pairwise FST distances The first component (PC1), which explained 2.48% of genetic variability, divided Caucasian tur into two groups. The first group consisted of E_TUR (PC1 < 0) and the second one included W_TUR and M_TUR (PC1 > 0).
The level of genetic differentiation was estimated with pairwise F ST distances (Table 1). We observed the highest F ST value between W_TUR and E_TUR-0.161. M_TUR was genetically closer to W_TUR than to E_TUR, but in both cases the differentiation was moderate-0.099 and 0.135, respectively.
The above-mentioned results were also supported with Neighbor-Net dendrogram ( Figure 6). All the studied individuals were placed in clades according to their geographical origin. This analysis also revealed that there were no close relatives among the 15 genotyped animals.
Expected heterozygosity (HE) and allelic richness (AR) in W_TUR were significantly higher than in E_TUR (for both indicators t-test p-values were less than 0.0001) ( Table 2). Slightly higher genetic diversity parameters were observed in M_TUR. The within population inbreeding coefficient (FIS) was close to zero in W_TUR and M_TUR which showed that these groups were in a fairly stable state, and slightly higher than zero in E_TUR, which indicated weak heterozygote deficit.  All the studied individuals were placed in clades according to their geographical origin. This analysis also revealed that there were no close relatives among the 15 genotyped animals.
Expected heterozygosity (H E ) and allelic richness (A R ) in W_TUR were significantly higher than in E_TUR (for both indicators t-test p-values were less than 0.0001) ( Table 2). Slightly higher genetic diversity parameters were observed in M_TUR. The within population inbreeding coefficient (F IS ) was close to zero in W_TUR and M_TUR which showed that these groups were in a fairly stable state, and slightly higher than zero in E_TUR, which indicated weak heterozygote deficit. Individual heterozygosity values, which were represented by multilocus heterozygosity (MLH) showed that all the E_TUR individuals had lower heterozygosity (from 0.217 to 0.253) than both W_TUR (from 0.273 to 0.282) and M_TUR (from 0.255 to 0.283) (Figure 7). Although mean values of multilocus heterozygosity were almost equal in W_TUR-0.277 ± 0.002 and in M_TUR-0.269 ± 0.005, the range of values was wider in M_TUR than in W_TUR. In E_TUR mean MLH was significantly lower-0.233 ± 0.006.

Discussion
In our study we assessed the applicability of Illumina Goat SNP50 BeadChip, developed for domestic goats (Capra hircus), for genetic studies of Caucasian tur (Capra caucasica). About 15 individuals that belong to three groups of Caucasian tur (West-Caucasian tur, East-Caucasian tur, and Mid-Caucasian tur) were genotyped with 5544 polymorphic loci (3,885 SNPs after LD pruning) distributed all over 29 autosomes. This number of polymorphic sites is sufficient for carrying out phylogenetic studies and assessing the genetic diversity. For instance, in the work on European bison only 960 SNPs were used [15], 1121 polymorphic loci were detected for snow sheep (Ovis nivicola) [21], and based on 7303 SNPs, population structure and genetic diversity of reindeer were investigated [18].
Regarding intraspecific taxonomy of Caucasian tur, there are two main questions to which there is still no answer: 1) Whether the West-Caucasian tur and the East-Caucasian tur are one species or should they be considered different species; 2) is the Mid-Caucasian tur a separate subspecies or a hybrid form. The works of 19th and 20th centuries, based on the morphological characteristics could not resolve these problems. The majority of authors such as Lydekker [2], Tsalkin [3], Geptner [4], Damm [5] suggested division into two species-C. cylindricornis and C. caucasica. Unlike the above-mentioned authors, Sokolov [42] and Tembotov [43] proposed to recognize one species-Capra caucasica with several subspecies C. c. severtzovi, C. c. cylindricornis, C. c. caucasica, and C. c. dinniki. Couturier [44] in 1962 suggested the presence of East West morphological cline in Caucasian tur and indicated that there is only one taxon in the Caucasus (C. aegagrus caucasica). Thirty years later Ayunts [45] showed a smooth transition in the variability of horns in tur, the main distinguishing feature of these populations and noted that Caucasian tur was represented by a number of slightly different groups and that the Although mean values of multilocus heterozygosity were almost equal in W_TUR-0.277 ± 0.002 and in M_TUR-0.269 ± 0.005, the range of values was wider in M_TUR than in W_TUR. In E_TUR mean MLH was significantly lower-0.233 ± 0.006.

Discussion
In our study we assessed the applicability of Illumina Goat SNP50 BeadChip, developed for domestic goats (Capra hircus), for genetic studies of Caucasian tur (Capra caucasica). About 15 individuals that belong to three groups of Caucasian tur (West-Caucasian tur, East-Caucasian tur, and Mid-Caucasian tur) were genotyped with 5544 polymorphic loci (3885 SNPs after LD pruning) distributed all over 29 autosomes. This number of polymorphic sites is sufficient for carrying out phylogenetic studies and assessing the genetic diversity. For instance, in the work on European bison only 960 SNPs were used [15], 1121 polymorphic loci were detected for snow sheep (Ovis nivicola) [21], and based on 7303 SNPs, population structure and genetic diversity of reindeer were investigated [18].
Regarding intraspecific taxonomy of Caucasian tur, there are two main questions to which there is still no answer: (1) Whether the West-Caucasian tur and the East-Caucasian tur are one species or should they be considered different species; (2) is the Mid-Caucasian tur a separate subspecies or a hybrid form. The works of 19th and 20th centuries, based on the morphological characteristics could not resolve these problems. The majority of authors such as Lydekker [2], Tsalkin [3], Geptner [4], Damm [5] suggested division into two species-C. cylindricornis and C. caucasica. Unlike the above-mentioned authors, Sokolov [42] and Tembotov [43] proposed to recognize one species-Capra caucasica with several subspecies C. c. severtzovi, C. c. cylindricornis, C. c. caucasica, and C. c. dinniki. Couturier [44] in 1962 suggested the presence of East West morphological cline in Caucasian tur and indicated that there is only one taxon in the Caucasus (C. aegagrus caucasica). Thirty years later Ayunts [45] showed a smooth transition in the variability of horns in tur, the main distinguishing feature of these populations and noted that Caucasian tur was represented by a number of slightly different groups and that the division into subspecies was not justified.
Currently in the era of molecular genetic studies, many taxonomic riddles were resolved. Clarification of the phylogeny in the subfamily Caprinae made it possible to determine the status of three subspecies of tahr, which are now considered different genera-Hemitragus, Nilgiritragus, and Arabitragus [46]. New subspecies of blue sheep (Pseudois nayar) [47] and snow sheep (Ovis nivicola) [21,48] were identified.
To date molecular genetic studies of Caucasian tur were based on the investigation of mitochondrial DNA and Y-chromosome. Manceau et al. [9]  In all the abovementioned studies strong East-West differentiation in the Caucasian tur was demonstrated and therefore all the authors supported evidence for the existence of two species-C. caucasica and C. cylindricornis. In the work of Zvychaynaya [11], some mtDNA and Ychromosome haplotype mixture was observed in the Central Caucasus that was explained as the evidence of C. caucasica and C. cylindricornis hybridization. In all the mtDNA studies the authors emphasized that further research on nuclear markers should be performed.
Our preliminary analysis of nuclear DNA markers confirmed clear genetic differentiation between all the three studied populations. PCA analysis showed a clear division into clusters, however the percent of genetic variation was very low-2.48. Pairwise F ST genetic distance between the most distant populations was only 0.161, which is too low for dividing them into species and even subspecies (e.g., in Yakut snow sheep (Ovis nivicola lydekkeri) pairwise F ST values ranged from 0.044 to 0.144 within the subspecies, and from 0.112 to 0.205 with the proposed new subspecies [21]). Moreover, this was also confirmed by cross-validation (CV) error calculated for the cluster analysis ( Figure S3), which suggested that K (optimum number of populations in our dataset) was 1 (Figure 8). In all the above-mentioned studies strong East-West differentiation in the Caucasian tur was demonstrated and therefore all the authors supported evidence for the existence of two species-C. caucasica and C. cylindricornis. In the work of Zvychaynaya [11], some mtDNA and Y-chromosome haplotype mixture was observed in the Central Caucasus that was explained as the evidence of C. caucasica and C. cylindricornis hybridization. In all the mtDNA studies the authors emphasized that further research on nuclear markers should be performed.
Our preliminary analysis of nuclear DNA markers confirmed clear genetic differentiation between all the three studied populations. PCA analysis showed a clear division into clusters, however the percent of genetic variation was very low-2.48. Pairwise FST genetic distance between the most distant populations was only 0.161, which is too low for dividing them into species and even subspecies (e.g., in Yakut snow sheep (Ovis nivicola lydekkeri) pairwise FST values ranged from 0.044 to 0.144 within the subspecies, and from 0.112 to 0.205 with the proposed new subspecies [21]). Moreover, this was also confirmed by cross-validation (CV) error calculated for the cluster analysis ( Figure S3), which suggested that K (optimum number of populations in our dataset) was 1 ( Figure 8). Mid-Caucasian tur was genetically a little closer to West-Caucasian tur (FST 0.99) than to East-Caucasian tur (FST 0.135) but we did not find any evidence that this population arose as a result of hybridization. Heterozygosity in Mid-Caucasian tur was a little lower than in West-Caucasian tur. The number of unique SNPs was even lower than in both West-Caucasian and East-Caucasian tur. The values of f3 statistics, that is a test revealing if the population of interest is the result of admixture between two other populations,  Mid-Caucasian tur was genetically a little closer to West-Caucasian tur (F ST 0.99) than to East-Caucasian tur (F ST 0.135) but we did not find any evidence that this population arose as a result of hybridization. Heterozygosity in Mid-Caucasian tur was a little lower than in West-Caucasian tur. The number of unique SNPs was even lower than in both West-Caucasian and East-Caucasian tur. The values of f3 statistics, that is a test revealing if the population of interest is the result of admixture between two other populations, were not significant and therefore also rejected the hypothesis of hybrid origin of Mid-Caucasian tur (Table S1). Based on the above we suggest that Caucasian tur genetic variation follows a clinal pattern [49] and therefore it should be considered a single species-Capra caucasica, which includes several ecotypes. The division into subspecies does not seem to be justified. Thus, the results of our investigation based on nuclear molecular genetic markers agree with Ayunts's study [45] on variability of horns and confirming the hypothesis proposed by Couturier [44]. Further research with the collection of more samples will help to give a final answer regarding the intraspecific taxonomy of the Caucasian tur, reconstruct the natural history of the species formation, and assess the genetic diversity.

Conclusions
Here we demonstrated that the Illumina Goat SNP50 BeadChip developed for domestic goats could be used as a useful tool for genetic studies of Caucasian tur. It was shown that West-Caucasian tur, East-Caucasian tur, and Mid-Caucasian tur clustered separately from each other. Genetic diversity parameters in East-Caucasian tur were lower than in the other groups. A clinal pattern of genetic variation was revealed. It was suggested to consider Caucasian tur a single species with several ecotypes. To obtain more accurate information about population structure and genetic diversity of Caucasian tur, further studies based on Illumina Goat SNP50 Bead-Chip with larger number of samples are needed.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/d13070312/s1. Figure S1: Proportion of polymorphic sites per chromosome in Caucasian tur (Capra caucasica). Figure S2: Minor allele frequency (MAF) in Caucasian tur (Capra caucasica). Figure S3: Admixture clustering results at K = 2 and K = 3 for the three populations of Caucasian tur. Table S1: The results of f3 statistics computed for the three populations of Caucasian tur. Table S2: Hardy-Weinberg equilibrium p-values and minor allele frequency in SNPs selected for the study of Caucasian tur (Capra caucasica).