Population genetics of the main population of brown bears in southwest Asia

Samples

We collected 154 samples from the Grand Kaçkar Mountains (GKM) region in two neighboring Turkish provinces, Artvin and Bayburt (Fig. 1; Table 1), respectively located on the eastern and the southwestern edges of the GKM. Both are small montane provinces in the Lesser Caucasus Region characterized by steep mountains and large fast-flowing rivers. We chose to sample in these provinces because current dam construction activities (including new roads and forest logging) may reduce connectivity between eastern and western GKM (sub)populations even further in the near future. Obtaining samples from these locations enabled us to assess population connectivity between areas separated by up to 300 km, a distance near the limit of the homing ability of brown bears (Conover, 2001). The bear population in Artvin is connected to bear populations in Georgia, and the bear population in Bayburt is connected to bear populations in Eastern and Central Anatolia (Ambarlı, Ertürk & Soyumert, 2016).

Noninvasive genetic sampling was carried out from 2008 to 2014 (from May to October). In the Artvin province, we collected hair samples from rubbing trees and power poles impregnated with creosotes (N = 105), corral hair-traps (N = 7), and from barbed wires placed on rubbing trees (N = 12), at altitudes between 700–2,500 m above sea level (Table 1). Additionally, we used three invasive hair samples from captured bears (N = 3). Samples from the Bayburt province were collected from rubbing trees and power poles (N = 20). Hair samples were placed in paper envelopes without contacting human skin and were then stored at room temperature in zip lock bags with silica gel (Roon, Waits & Kendall, 2003). In addition, fresh tissue samples (N = 7) were obtained from bears captured in the Artvin province between 2010 and 2011 (Permission No:B.23.0.DMP.0.13.02–445.05-36125, Ministry of Forestry and Water Affairs) (Ambarlı, 2012), and were stored in 98% ethanol. In total, we collected 147 hair samples (Noninvasive N = 144, invasive N = 3) (Artvin N = 127, Bayburt N = 20) and seven tissue samples (Artvin N = 7) (Table 1).

DNA extraction and genotyping

For both hair and tissue samples, DNA was extracted using a commercially available kit (GEN-IAL GmbH, Troisdorf, Germany) following the manufacturer’s instructions without any modifications.

We genotyped all samples at 10 microsatellite loci: G10C, G10X, UarD1585, UarT739, UarD3139, UarD3684, G1D, Mu05, Mu23, Mu50, and at the sexing locus SRY (primer details in Table S1). One primer of each primer pair was 5′  labeled with a fluorescent dye (6-FAM, HEX, or NED), so that differently labeled primer pairs could be used in a multiplex approach for genotyping. Loci were amplified using the Qiagen Multiplex PCR Kit (Qiagen, Hilden, Germany) in three multiplexes of 10 µL final reaction volume, following the manufacturer’s recommended conditions (multiplexes detailed in Table S1).

Following Miller, Joyce & Waits (2002), we used a Maximum Likelihood approach for genotyping noninvasively collected samples: each sample was genotyped in parallel (N = 2 replicates per sample) and if a mismatch occurred between replicates, results were discarded and genotyping was repeated on a new DNA extraction of that sample, resulting in 2–4 replicates per sample. The low number of hairs collected per sample limited us to two DNA extractions per sample. Thus, if a mismatch occurred in the second round of parallel genotyping, the sample was discarded.

As bears are highly mobile, we had to consider multiple sampling from a single individual. To remove such bias, we looked for matching genotypes using the R package (R Core Team, 2015) allelematch (v2.5; Galpern et al., 2012) applying the option ‘alleleMismatch=2’. Matching genotypes were always assigned to the same individual, mismatching genotypes were only assigned to the same individual if either one of the following criteria was met: (i) only one allele difference was observed between genotypes and could be attributed to large allele drop-out or allele-size shift due to stuttering (N = 3), or (ii) two allele differences between genotypes were observed but could be attributed to large allele drop-out, allele-size shift or missing data (N = 6). Genotypes from tissue samples were included in this analysis.

Genotypic analyses

The probability of null alleles being present in the data set was assessed using Micro-Checker (v.2.2; Van Oosterhout et al., 2004). To assess if the number of loci and alleles was sufficient to discriminate individuals, we estimated the cumulative values of the unbiased probability of identity (P_IDunb) and of the probability of identity given siblings (P_IDsib) using the software package Gimlet (v.1.3.3; Valière, 2002).

We examined the genotypes of the GKM bear population using the following software packages: Popgene (v.1.32; Yeh et al., 1997) to test for deviation from Hardy–Weinberg equilibrium, and to estimate observed heterozygosity (H_O) and expected heterozygosity (H_E), and fstat (v.2.9.3.2; Goudet, 2002) to estimate the inbreeding coefficient F_IS.

Genetic structure among sampling sites was examined using multiple approaches because these can differ in their assumptions about markers or populations, such as random mating, absence of selection or absence of mutation. We thus included individual-based analyses that do not make such assumptions, as well as analyses that infer population genetic structure using allele frequencies that do make such assumptions.

The following analyses are individual-based and are not be impacted by, for example, uneven sampling of populations or deviation from HWE. The R package memgene (v.1.0; Galpern et al., 2014) was used to calculate the proportion of shared alleles among samples (following Bowcock et al., 1994), and the package ape (v.5.1; Paradis, Claude & Strimmer, 2004) to construct unrooted neighbor-joining trees from the distance matrix. We used adegenet (v.2.0.1; Jombart, 2008) to conduct principal component analyses (PCA). PopGenReport (v.2.1; Adamack & Gruber, 2014) was used to calculate pairwise genetic dissimilarity between samples (following Kosman & Leonard, 2005) and to plot these values against the geographic distance between samples. This allowed us to retain information regarding the resampling of individuals (i.e., unique genotypes) in a geographic context, as it was conducted on the dataset consisting of 71 samples that had been scored at a sufficient number of microsatellite loci (Table 1). Thus, for example, the genotypes that had been detected at more than one location could be included, together with spatial data regarding sampling location. Native functions of R were also used to plot graphs.

The following analyses infer population genetic structure using allele frequencies that do make such assumptions, and can be impacted by deviations from these assumptions (e.g., if they are not in HWE, by Wahlund effect). Arlequin (v.3.5.2.1; Excoffier & Lischer, 2010) was used to calculate population pairwise F_ST; significance was tested by permuting individuals between populations (10,000 permutations).

We used Bayesian inference implemented in the software Structure (v.2.3.4; Pritchard, Stephens & Donnelly, 2000; Falush, Stephens & Pritchard, 2003) to estimate the number of genotypic clusters in our study area. We used the admixture model as prior and ran the software for 600,000 steps following a burn-in of 200,000 steps. The most likely number of genotypic clusters (K) was evaluated in a range from K = 1 to K = 6 and was assessed in ten independent runs per K. We also used Bayesian inference implemented in the software tess (v.2.3; Chen et al., 2007) to infer population structure, and included spatial data regarding sampling location. We used the admixture model and ran the software for 100,000 sweeps, discarding the first 50,000 as burn-in. We ran tess 100 times for each K in a range from K = 2 to K = 6, and averaged results (per K) over the 20% runs with the lowest Deviance Information Criterion (DIC). Results were visualized using clumpp (v.1.12; Jakobsson & Rosenberg, 2007) and distruct (v.1.1; Rosenberg, 2004).

To examine the impact of sampling date and sample source on genotyping success, we compared results for hair samples collected from different noninvasive sources (N = 144) and years (2008–2014) using Chi-square and contingency tables. Understanding how age and source of samples impact genotyping success can inform improving sampling strategy for further studies.

Microsatellite variation

The microsatellite loci used were all highly polymorphic, with a range of seven to ten alleles (Table 2). In total, 64 out of the 147 hair samples (43.5%) yielded data for at least eight out of the ten microsatellite loci (Table 1). Applying the criteria outlined above (Materials & Methods), we scored 48 unique genotypes among these 64 hair and the seven tissue samples (71 samples total). Based on the SRY data, these consisted of 40 males (83.3%) and eight females (16.7%). However, we also found that SRY worked for a tissue sample from a captured female bear and for one hair sample (H58) that had four exact matches to a female (Samples H20, H51, H110, H114, see dataset in the Supplemental Files). We detected 14 unique genotypes (12 males and two female genotypes) during a single sampling period (2010) in a small part of the Artvin sampling area (<100 km²).

In order to reduce the impact of missing data on downstream analyses, we discarded six unique genotypes because there were missing data at more than one locus (Table 1). Out of the remaining 42 unique genotypes, 36 were detected in the Artvin province, five were detected in the Bayburt province, and one genotype was observed in both provinces (Table 1).

Evaluation of microsatellite allele distribution data using MicroChecker revealed that five loci had a significant probability for the presence of null alleles (G10C, UarD1585, G1D, Mu05, and Mu23; Table 2). These five loci also displayed significant levels of inbreeding (measured as F_IS), and three of them (G10C, UarD1585, and G1D) showed significant deviation from Hardy–Weinberg equilibrium (Table 2).

While inclusion of these loci did not alter the results of analyses on population structure, we restricted analyses to the five remaining loci (Mu50, UarT739, G10X, UarD3139, UarD3684) as these were still capable of distinguishing 41 out of the 42 unique genotypes identified using the full set of markers. Reducing the number of loci resulted in the cumulative estimates of ‘probability of identity’ dropping from P_IDunb = 2.7 × 10⁻¹³ and P_IDsib = 5.2 × 10⁻⁵ (10 loci) to P_IDunb = 8.7 × 10⁻⁷ and P_IDsib = 7.7 × 10⁻³ (5 loci).

Over all genotypes (N = 42), expected and observed heterozygosity was high (H_E = 0.78, H_O = 0.74; Table 2). When considering genotypes detected in only one province (i.e., excluding the genotype found in both provinces), we observed virtually identical values for genotypes from the Artvin province (N = 36, H_E = 0.79, H_O = 0.74), and slightly lower values for genotypes from the Bayburt province (N = 5, H_E = 0.68, H_O = 0.72).

Population structure

Despite only obtaining a moderate number of genotypes from the Bayburt province, we nevertheless attempted to resolve brown bear population structure in our study area. However, excluding the genotype found in both provinces, we found no genetic differentiation between the brown bear populations from the Artvin and Bayburt provinces (F_ST = 0.004, P = 0.47; N = 41). When examining the proportion of shared alleles between individuals (following Bowcock et al., 1994; N = 42), no clustering of genotypes based on sample origin was apparent (Fig. 2A). Similarly, the principal component analysis (PCA; N = 42) did not show a clustering of genotypes based on sample origin (Fig. 2B). Bayesian inference using Structure (N = 42) suggested the presence of a single population (Fig. 2C; for Structure plots of K = 1 to K = 6 see Fig. S1).

To take spatial information into account, we firstly examined genetic dissimilarity versus geographic distance (Fig. 2D). Within the geographic range of the study, we found both samples with no genetic dissimilarity at short distances (up to ∼15 km apart, resampled genotypes within provinces), and at long distances (160 km, the genotype detected in both provinces, indicated by an arrow in Fig. 2D). Pairwise values among samples from the two different provinces (geographic distance of ∼150 km to 300 km) were mostly on the same order as pairwise values for samples from the same province (geographic distance of 0 km to ∼100 km). Thus, genetic dissimilarity did not increase with increasing geographic distance at the spatial scale of our study (horizontal line in Fig. 2D).

We also assessed population structure by Bayesian inference implemented in tess, incorporating spatial data. As this software does not calculate the likelihood of all individuals belonging to a single population (K = 1), we conducted the analysis for K = 2 to K = 6. Neither the Deviance Information Criterion (DIC) nor the plots of individual cluster membership provided unambiguous evidence for a population subdivision. Only up to five genotypes had a cluster membership probability above 0.8 (for K = 2) in any analysis (Fig. S2).

Source of hair samples

Hair samples collected from rubbing trees (with or without barbed wire) yielded DNA extracts that performed significantly better for genotyping than did samples collected from power poles treated with creosotes (Chi-square = 15.81, Df = 1, P < 0.001) regardless of the collection year (Tables 3 and 4). Overall, we observed a significant difference in genotyping success among different noninvasive sources of hair samples (Chi-square = 19.09; Df = 3; P = 0.0026) (Table 4).

Genetic diversity

We found high levels of allelic variation among bears in our study area (H_E = 0.78, H_O = 0.74; Table 2), which were also apparent when the two areas were examined separately (Artvin province: N = 36, H_E = 0.79, H_O = 0.74; Bayburt province: N = 5, H_E = 0.68, H_O = 0.72). While our estimates are not directly comparable with those obtained from other studies employing different microsatellite loci (Table 5), they do provide the first evidence that the brown bears in the GKM in northeastern Anatolia are not genetically impoverished. Measures of genetic diversity for threatened populations, for example for the brown bear population in Spain (Table 5), are far lower than those we report here. As the main source population for brown bears in northeastern Anatolia and surrounding regions, the GKM population represents an important genetic reservoir that needs to be preserved.

Unfortunately, we were not able to satisfactorily resolve population structure in our study area because we were only able to obtain a moderate number of genotypes from the Bayburt province. None of the analyses could reject the hypothesis of a single population in our study area (e.g., Fig. 2), and we did detect one bear in both provinces, indicating that (current) landscape features do not hinder the movement of brown bears between these two areas. However, due to the uneven sampling between areas, which can be problematic for some analyses of population structure (e.g., Puechmaille, 2016), we prefer to refrain from definitively stating that there is no population subdivision. Consequently, we cannot claim that the bears in the GKM can be managed as a single conservation unit. Further work is clearly needed to address this, particularly with respect to additional sampling; ideally, also increasing the geographic scope of the current study.

Genetic monitoring of this population is also important to gauge the impact of habitat loss and habitat fragmentation in the region (Şekercioğlu et al., 2011; Balkenhol & Waits, 2009). These result from the ongoing development of infrastructure (big dams, HEPPs, roads and high-voltage transmission pylons) in areas where still continuous pristine habitat exists (e.g., natural old forests) for large carnivores. Extensive loss of habitat, and concomitant fragmentation of the remaining habitat, has already been reported for Eurasian lynx, golden jackals and wild goats in the same region (Özdemirel, Turak & Bilgin, 2016). Due to the importance of the brown bear population in the GKM for on-going and future conservation actions, the impacts of habitat alterations need to be understood, and when possible minimized, in order to preserve both the current genetic variation in these brown bears and the connectivity of their populations.

Sources of noninvasively collected samples

It was recently shown that hair samples obtained from rubbing trees perform better for genotyping than hair samples collected by other noninvasive methods, such as hair from corral traps (Berezowska-Cnota et al., 2017). We similarly observed that hair samples from rubbing trees performed best among the hair sample sources used in our study (Table 4). Our hair sampling strategy relied mostly on rubbing objects (such as trees and power poles) and barbed wire on rubbing trees. The significant reduction in genotyping success of samples collected from power poles that we observed may be due to a higher proportion of shed guard hair without follicles that got stuck to the creosote. Moreover, hair on power poles were generally exposed to direct sunlight for long periods of time, which most probably caused higher levels of DNA degradation (Stetz, Seitz & Sawaya, 2014), when compared to hair obtained from trees where direct sunlight is almost absent. While power poles with creosote should not be discounted as a source for hair samples in future studies, it is advised to supplement these and rubbing trees with other sources from which samples could be collected noninvasively. For example, adding some barbed wire to power poles may increase the retention of underfur hair.

As reported for North American bears (Lamb, Walsh & Mowat, 2016), rubbing trees in the GKM also appear to be mostly used by male bears as only 12 of 71 genotypes belongs to female bears (Table S1). Thus, population estimates relying on noninvasive sampling from rubbing trees will be male-biased. As we detected discrepancies in the sexes of individuals when using only SRY, using a more sophisticated sexing procedure (e.g., a multilocus approach with internal controls/validation, Sastre et al., 2009) would be safer for noninvasive samples to determine sex of individuals. Considering this, and the fact that female brown bears in Turkey have very small home ranges (Ambarlı, 2012), greater effort needs to be allocated to sample female bears with cubs, as well as sub-adult bears. This may be accomplished by including fecal samples in the analyses or by establishing hair traps in more remote areas.