GPS tracking and movement modeling confirm the existence of two life history strategies.
To determine the movement behavior of brown bears in Eastern Türkiye we used GPS-collar data from 31 brown bears (10 female, 21 male) captured around the Sarıkamış, Kars province of Türkiye (Fig. 1) and tracked successfully for an average of 370 days (range: 126–726) (Table S1). From this data, we calculated the net squared displacement (NSD) of each individual by measuring the square of the distance from the animal's starting point to successive daily locations (Bunnefeld et al., 2011). The NSD data were then integrated into a nonlinear hierarchical modeling framework using the R package migrateR (Spitz et al., 2017) to classify individuals into five a priori non-linear models, representing (1) year-round residency, (2) nomad, (3) disperser, (4) migrant, and (5) mix-migrant movement behavior (Bunnefeld et al., 2011; Börger & Fryxell, 2012).
Out of the 31 bears, 15 were classified as resident, four as migrant, nine as mix-migrant, and two as nomadic. None of the individuals were classified as dispersers (Table S1, Movement Figures S1-S31). Following Cozzi et al. (2016), and to simplify genetic analyses we grouped nomadic and resident bears into a single category called “sedentary”, and migrant and mix-migrants into a second category called “migratory” (Table S1). Out of the 16 individuals modeled for movement patterns by Cozzi et al. (2016), we converged on the same results for 12 of them (Table S1). Two individuals (BTR009 and BTR018) classified as migratory in Cozzi et al. (2016) were re-classified as sedentary as these bears showed no evidence of leaving the Sarıkamış forest (Table S1). For the remaining two individuals (BTR016 and BTR019), GPS data were not complete in Cozzi et al. (2016) and thus were not classified previously. Our analysis of NDS patterns confirms the presence of two different life-history strategies, sedentary vs. migratory, in Eastern Türkiye brown bears.
RADseq provides high quality genome-wide information for comparing Eastern Türkiye brown bears to their conspecifics.
To produce a high-quality genomic resource for the 31 GPS-tracked brown bears, we sampled the genome of each individual by paired end RAD sequencing (RADseq; Miller et al., 2007; Baird et al., 2008; Etter et al., 2012) using the Sfb1 restriction enzyme and the Illumina sequencing platform. Sequenced reads were aligned to the chromosome length assembly of the American black bear genome, which had the most complete assembly among available bear genomes (Dudchenko et al., 2017; Dudchenko et al., 2018; Srivastava et al., 2019).
Mapping success to the reference genome was high with over 99% of the reads mapped per individual (Table S1). Additionally, whole-genome sequences (WGS) of brown bears from Georgia, Russia, Greece, Slovakia, Slovenia, Austria, and Alaska (Cahill et al., 2013; Benazzo et al., 2017; Barlow et al., 2018) were downloaded from the Sequence Read Archive (SRA) at NCBI (Table S2) and aligned to the same reference genome.
In total we interrogated 21,745,475 sites that overlapped between the RADseq and WGS data, had an average per individual read depth of 6X and were present in at least 50% of the eastern Türkiye brown bears. Polymorphic sites (i.e., SNPs) were identified using the probabilistic framework in the program ANGSD (Korneliussen et al., 2014) and we discovered 108,777 SNPs with a minimum minor allele frequency (MAF) of 0.05 and a significance threshold of P < 10 − 06.
Eastern Türkiye brown bears are characterized by distinct genetic structure, mixed ancestry and high diversity.
To determine genetic structure and ancestry of Eastern Türkiye brown bears together with their conspecifics around the world, we conducted PCA, and calculated admixture proportions based on genotype likelihoods (Skotte et al., 2013; Meisner & Albrechtsen, 2018) using the genome-wide data described above. Eastern (Türkiye, Iran, and Georgia) and Apennine brown bears were separated from the remaining populations (Europe, Russia and Alaska) along PC1, while the Russian and Alaskan samples further separated from the European samples (Greece, Slovakia, Slovenia and Austria) along PC2 (Fig. 2A). Although clustered together relative to other populations and geographically close, the Iranian and Georgian brown bear samples were still highly distinct from the Eastern Türkiye brown bears (Fig. 2A). A separate PCA using only the Eastern Türkiye samples revealed minimal structure within this population and showed no apparent difference between sedentary and migratory individuals (Fig. 2B).
Admixture analysis gave similar results to those of the PCA and determined K = 5 as the optimal number of clusters based on Evanno’s method (Figure S1, Evanno et al., 2005). The three main groups along PC1 (Apennine; Europe + Russia + Alaska; and Eastern brown bears) separated into distinct ancestry classes with only the Georgian and Iranian bears showing slight levels of shared ancestry with European (Greece, Slovakia, Slovenia, and Austria) bears at K = 5 (Fig. 2C). Interestingly, Eastern brown bears (Eastern Türkiye, Iran, Georgia) showed higher amounts of mixed ancestry when compared to the rest of the world populations (Apennine, Europe, Russia, and Alaska) which showed mono ancestry even at higher K values (Fig. 2C). Finally, admixture analysis captured no difference between sedentary and migratory bears within Eastern Türkiye (Fig. 2C).
To summarize the genetic diversity of Eastern Türkiye brown bears in comparison to its conspecifics around the world, we calculated Watterson’s theta (θw) across the genome for the 37 scaffolds that were over 25Mb using the thetaStat module in ANGSD. Excluding the Apennine bear population, genome-wide average nucleotide diversity (θw) was similar between most populations and ranged between 0.001 to 0.002 (Figure S2). Across the 37 scaffolds, the highest diversity was observed in the Eastern Türkiye population, followed by the Iran population (Fig. 3). As expected, the Apennine bear population had significantly lower diversity than all other populations (Fig. 3).
We conclude that Eastern Türkiye brown bears are characterized by distinct genetic structure, mixed ancestry, and higher levels of genetic diversity even when compared to geographically close populations. However, we did not capture any broad scale genetic differences between sedentary and migratory bears in Eastern Türkiye.
High genetic differentiation and distinct genotypes between sedentary and migratory bears at outlier loci.
To determine genomic regions associated with sedentary and migratory behavior we performed an association test using the score statistics (Skotte et al., 2012) implemented in ANGSD. In total we compared 12,207 SNPs between movement categories using a likelihood ratio test (LRT) with a genome-wide significance threshold of 4.09x10− 6 (Bonferroni-corrected α = 0.05; LRT > 21.2). The strongest associated SNP in our panel reached a significance level of 1.4x10− 5 (LRT = 18.8) and we detected no SNP that reached genome-wide significance (Table S3).
To determine if any of the top 20 SNPs (LRT > 13, P < 3.1x10− 4) in our association panel could be candidate regions for behavioral differences between sedentary and migratory bears we performed an outlier test based on the population branch statistics (PBS; Yi et al., 2010) to test if any of these sites showed extreme genetic differentiation in either of the two groups. The mean PBS score across the genome was 0.0205 ± 0.0330 for sedentary bears and 0.0127 ± 0.0277 for migratory bears. Among the 20 top associated SNPs, five fell over the 99.9 percentile range (PBS > 0.3547) in sedentary bears while another five fell over the 99.9 percentile range (PBS > 0.3106) in migratory bears (Fig. 4, Table S3). Moreover, all SNPs with PBS scores over the 99.9 percentile in either sedentary or migratory bears were at least 10 standard deviations above their respective means (Table 1), providing good evidence for the outlier nature of these loci. Out of the 10 highly differentiated SNPs only one, SNP_TM3 (differentiated in migratory bears), fell within a protein-coding gene (within the first exon of CCRL2). For the remaining highly differentiated SNPs, we found 3 genes (Myo5b, RAN, and TRDMT) lying within 1MB distance from SNP_TM4, SNP_TS3 and SNP_TS5, respectively (Table 1).
To further understand genomic patterns of differentiation at the 10 SNPs showing extreme PBS values, we genotyped each individual at these sites and categorized individuals as either heterozygous or homozygous for the ancestral or derived allele (Fig. 5). Genotype calling at these sites revealed three blocks with distinct genotype patterns between sedentary and migratory bears (Fig. 5). Block 1 and 2 constituted regions where the frequency of the derived allele was higher in sedentary bears, with most sedentary bears being heterozygous (block 1) or homozygous for the derived allele (block 2) in contrast to migratory bears which were mostly homozygous for the ancestral allele in block 1 and heterozygous in block 2 (Fig. 5). In contrast, Block 3 had the opposite pattern where the frequency of the derived allele was higher in migratory bears, with migratory bears being mostly heterozygous and sedentary bears being mostly homozygous for the ancestral allele (Fig. 5).
Even though not statistically significant at the genome-wide level, ten of the top 20 SNPs in our association panel showed extreme differentiation in either one of the two behavioral groups and were able to separate sedentary and migratory bears into distinct genotypes. This suggests that the differences in sedentary and migratory behavior might have an important genetic component.
Figure SEQ Figure \* ARABIC 1 The study area in northeastern Türkiye. Blue dot indicated the Sarıkamış garbage dump.
Table 1
Derived allele frequencies of sedentary, migratory and European bears at the ten outlier loci. SNP_TM indicates outlier loci in migratory bears while SNP_TS indicates outlier loci in sedentary bears. The distance to the nearest genes within a 1 MB window of the SNPs are also indicated along with the mean, global, SNP window and 99.9 percentile PBS scores.
SNP Name | Position | SNPs | DAF_Mig | DAF_Sed | DAF_Eur | Mean PBS | Global PBS | 99.9 percentile PBS | SNP window Pbs | Distance to Nearest Gene | Gene |
SNP_TM1 | HiC_scaffold_2-47325857 | SNP14 | 0.0453 | 0.3940 | 0.8072 | 0.0127 | 0.0083 | 0.3108 | 0.6202 | - | LTR |
SNP_TM2 | HiC_scaffold_18-2297497 | SNP11 | 0.3942 | 0.8558 | 0.7709 | 0.0127 | 0.0083 | 0.3108 | 0.3381 | - | LTR |
SNP_TM3 | HiC_scaffold_28-42755955 | SNP6 | 0.7393 | 0.1953 | 0.1835 | 0.0127 | 0.0083 | 0.3108 | 0.4591 | Inside Exon | CCRL2 |
SNP_TM4 | HiC_scaffold_30-24131228 | SNP5 | 0.5368 | 0.1254 | 0.0000 | 0.0127 | 0.0083 | 0.3108 | 0.3175 | 46,501 bp | Similar to Myo5b |
SNP_TM5 | HiC_scaffold_36-51993011 | SNP3 | 0.0000 | 0.3334 | 0.2254 | 0.0127 | 0.0083 | 0.3108 | 0.3799 | 4,935 bp | ELOC |
SNP_TS1 | HiC_scaffold_4-29070918 | SNP1 | 0.5405 | 0.0710 | 0.5115 | 0.0205 | 0.0143 | 0.3547 | 0.4579 | 42,113 bp | RAN |
SNP_TS2 | HiC_scaffold_14-65416581 | SNP17 | 0.4955 | 0.9011 | 0.1928 | 0.0205 | 0.0143 | 0.3547 | 0.5830 | - | LTR |
SNP_TS3 | HiC_scaffold_30-12480299 | SNP4 | 0.4170 | 0.8887 | 0.1927 | 0.0205 | 0.0143 | 0.3547 | 0.6292 | - | LTR |
SNP_TS4 | HiC_scaffold_34-17119378 | SNP19 | 0.4557 | 0.8770 | 0.4691 | 0.0205 | 0.0143 | 0.3547 | 0.3620 | 6,926 bp | TRDMT |
SNP_TS5 | HiC_scaffold_37-70763847 | SNP12 | 0.5032 | 0.9185 | 0.3723 | 0.0205 | 0.0143 | 0.3547 | 0.5442 | - | LTR |
Figure SEQ Figure \* ARABIC 5 SNPs that have the highest log-likelihood ratios at each position for each individual sedentary and migratory bear, along with the genotypes of the associated sites. Blue rectangles represent ancestral, yellow rectangles show heterozygous, and green ones represent derived alleles. Gray rectangles are the ones where the genotype is not known.