**3. 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). Total genotyping rate for the 15 samples of Caucasian tur was 94.8%. The selected polymorphic loci were distributed all over 29 autosomes. The lowest number of SNPs was found on the 25th chromosome—68, and the highest on the 1st chromosome—348 (Figure 3A). A similar pattern was observed for the three populations examined separately (Figure 3B–D), (Figure S1).

**Figure 3.** Number of polymorphic sites per chromosome in Caucasian tur (*Capra caucasica*). (**A**). Number of SNPs in all the studied populations. (**B**). Number of SNPs in West-Caucasian tur. (**C**). Number of SNPs in East-Caucasian tur. (**D**). Number of SNPs in Mid-Caucasian tur.

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).

**Figure 4.** Venn diagram representing the number of unique and shared polymorphic SNPs in the three groups of Caucasian tur (*Capra caucasica*). W\_TUR = West-Caucasian tur, E\_TUR = East-Caucasian tur, M\_TUR = Mid-Caucasian tur.

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).

**Figure 5.** Principal component analysis (PCA) of the groups of Caucasian tur (*Capra caucasica*): (**A**) the first two components (PC1 and PC2), (**B**) the first and third components (PC1 and PC3). W\_TUR = West-Caucasian tur, E\_TUR = East-Caucasian tur, M\_TUR = Mid-Caucasian tur.

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).

**Table 1.** Pairwise *F*ST genetic distances between the groups of Caucasian tur (*Capra caucasica*).


Notes: W\_TUR = West-Caucasian tur, E\_TUR = East-Caucasian tur, M\_TUR = Mid-Caucasian tur. *F*ST values in the range 0–0.05 indicate low genetic differentiation; value between 0.05 and 0.15, moderate differentiation; values between 0.15 and 0.25, high differentiation; and values above 0.25, very high genetic differentiation [40,41].

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).

**Figure 6.** Neighbor-Net dendrogram of 15 Caucasian tur (*Capra caucasica*) individuals, based on IBS distances.

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.


**Table 2.** Genetic diversity in the groups of Caucasian tur (*Capra caucasica*).

Notes: *n*—number of samples, *H*O—observed heterozygosity, *H*E—expected heterozygosity, *F*IS—inbreeding coefficient, *A*R—allelic richness, SE—standard error, CI—confidence interval, W\_TUR = West-Caucasian tur, E\_TUR = East-Caucasian tur, M\_TUR = Mid-Caucasian tur.

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).

**Figure 7.** Individual multilocus heterozygosity (MLH) in the groups of Caucasian tur (*Capra caucasica*).

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.

#### **4. 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] examined concatenated fragments of cytochrome b and control region (500 bp in length) of three samples of East-Caucasian tur and nine samples of West-Caucasian tur. Kazanskaya et al. [10] examined fragments of cytochrome b and control region of eight samples of East-Caucasian tur and nine samples of West-Caucasian tur. In the research of Zvychaynaya [11] the cytochrome b fragment (1128 bp) and Y-chromosome intron of gene SRY (1832 bp) of West-Caucasian tur (*n* = 6) and East-Caucasian (*n* =6) tur were analyzed. Kashinina and Kholodova [12] examined a 715 bp fragment of cytochrome b gene in six samples of East-Caucasian tur, 11 samples of West-Caucasian tur and eight samples of Mid-Caucasian tur. 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).

**Figure 8.** Admixture cross-validation (CV) error indicating the optimum number of clusters in the dataset.

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.
