Genetic differentiation and asymmetric gene flow among Carpathian brown bear (Ursus arctos) populations—Implications for conservation of transboundary populations

Abstract The abundance and distribution of large carnivores in Europe have been historically reduced. Their recovery requires multilevel coordination, especially regarding transboundary populations. Here, we apply nuclear and mitochondrial genetic markers to test for admixture level and its impact on population genetic structure of contemporary brown bears (Ursus arctos) from the Eastern, Southern, and Western Carpathians. Carpathian Mountains (Europe). Nearly 400 noninvasive brown bear DNA samples from the Western (Poland) and Eastern Carpathians (Bieszczady Mountains in Poland, Slovakia, Ukraine) were collected. Together with DNA isolates from Slovakia and Romania, they were analyzed using the set of eight microsatellite loci and two mtDNA regions (control region and cytochrome b). A set of 113 individuals with complete genotypes was used to investigate genetic differentiation across national boundaries, genetic structuring within and between populations, and movement between populations. Transboundary brown bear subpopulations (Slovakia and Poland) did not show significant internal genetic structure, and thus were treated as cohesive units. All brown bears from the Western Carpathians carried mitochondrial haplotypes from the Eastern lineage, while the Western lineage prevailed in the brown bears from the Bieszczady Mountains. Despite similar levels of microsatellite variability, we documented significant differentiation among the studied populations for nuclear markers and mtDNA. We also detected male‐biased and asymmetrical movement into the Bieszczady Mountains population from the Western Carpathians. Our findings suggest initial colonization of the Western Carpathians by brown bears possessing mtDNA from the Eastern lineage. Genetic structuring among populations at microsatellite loci could be a result of human‐mediated alterations. Detected asymmetric gene flow suggests ongoing expansion from more abundant populations into the Bieszczady Mountains and thus supports a metapopulation model. The knowledge concerning this complex pattern can be implemented in a joint Carpathian brown bear management plan that should allow population mixing by dispersing males.


| INTRODUC TI ON
The brown bear (Ursus arctos L.) is the largest terrestrial carnivore with a wide Holarctic distribution. Its populations formerly occupied all of Europe (Swenson, Gerstl, Dahle, & Zedrosser, 2000;Swenson, Taberlet, & Bellemain, 2011). However, since the 19th century, habitat destruction and human persecution have led to a severe decline among European brown bear populations, resulting in local extirpations and population fragmentation (Servheen, 1990). In western and central Europe, few contemporary populations are isolated from each other and from the western limit of continuous brown bear range located in Scandinavia, Estonia, and western Russia (Zedrosser, Dahle, Swenson, & Gerstl, 2001). The Carpathian population of brown bears spans national boundaries of several countries, consequently, its conservation and management falls under different jurisdictions. Therefore, the challenge of maintaining this wide-ranging carnivore is heightened (Chapron et al., 2014;Swenson et al., 2000Swenson et al., , 2011. The importance of this population is even greater because the Carpathians could have served as a refuge area or as a crucial movement corridor for brown bears, which led to the rise of brown bear populations in eastern and northern Europe during or after the last ice age (Saarma et al., 2007).
In addition to natural processes, human-mediated alterations resulted in severe brown bear population bottlenecks (Taberlet & Bouvet, 1994). Local persecutions, significant habitat fragmentation, followed by intentional translocations could have a disproportional role in shaping genetic diversity of this wide-ranging carnivore (Bray et al., 2013;Crispo, Moore, Lee-Yaw, Gray, & Haller, 2011).
Carpathian Mountains population is the largest brown bear population in Europe outside the continuous distribution range (Figure 1a; Chapron et al., 2014). However, in the early 20th century, it experienced a severe decline due to extensive deforestation and overhunting which led to the isolation of a small population of brown bears in Western Carpathians from the rest of the population (Finďo, Skuban, & Koreň, 2007;Hartl & Hell, 1994). At the end of World War I, the Western Carpathian population of brown bears decreased to roughly 15-75 individuals (after Straka, Paule, Ionescu, Štofík, & Adamec, 2012). In contrast, the East-Carpathian population, resident in Romania, Ukraine, South-Eastern Poland, and North-Eastern Slovakia, never dropped below 800 individuals (Ionescu, 1999).
However, its westernmost portion (Poland and Slovakia) was limited to only several individuals between World War I and World War II (Jakubiec, 2001;Sabadoš & Šimiak, 1981). Changes in management policies and protective legislation have allowed for the restoration of the Carpathian brown bear population (Chapron et al., 2014). In the West-Carpathian population, recent estimates indicate that the population is stable (approx. 1,250 individuals; Kaczensky et al., 2012;Paule et al., 2015), with a majority of bears in Slovakia and approximately a dozen or so in Poland. The East-Carpathian population in the Bieszczady Mountains is much smaller with approximately 83 individuals in Poland (Śmietana et al., 2014) and at least 15 individuals in Slovakia (Straka, Štofík, & Paule, 2013). The Ukrainian part of the East-Carpathian population is occupied by at least 400 brown bears (Zedrosser et al., 2001). The largest population of the Carpathians, which includes portions of the Eastern and Southern Carpathians, harbors more than 6,000 brown bears in Romania (Chapron et al., 2014;Kaczensky et al., 2012).
Population genetic studies using nuclear microsatellite markers revealed the presence of strong genetic structuring between populations from Western Slovakia, Eastern Slovakia, and Romanian Carpathians (Straka et al., 2012). The authors concluded that the present genetic differentiation is a result of nearly 100-year isolation of these geographically close populations; however, they did not estimate levels of gene flow among these populations. In this study, mitochondrial and autosomal markers were used to measure genetic variability within and differentiation among brown bear populations from Western Carpathians (WC; Poland and Slovakia), Bieszczady Mountains (BM; Poland and Slovakia), Ukraine (UKR; mtDNA only), and Romania (ROM). We also estimated the level and direction of admixture that might be indicative for both natural male-biased K E Y W O R D S brown bear, Carpathians, conservation, genetic structure, phylogeography, transboundary populations dispersal and translocations. Moreover, we tested for a possible asymmetry in gene flow, as it could be an indicator of a metapopulation model of genetic connectivity during the recolonization process. Finally, since brown bear populations in Poland (POL) and Slovakia (SVK) are transboundary populations that experience different management practices (the brown bear has been strictly protected in Poland since 1952, while in Slovakia the species is managed by culling), we compared subsamples that originated from these two countries for conservation purposes.

| Samples
In our analysis, DNA was extracted from the following sources: scats (n = 244); hairs found on marking trees and fences crossed by bears (n = 134); dried blood found on a vehicle destroyed by a bear; 14 buccal swabs or hairs collected from bears captured for telemetry studies; an individual rescued from poacher's snare; another one captured due to public safety; and from four dead individuals occasionally found in the field. Samples were collected in Polish (n = 376), Slovakian (n = 20), and Ukrainian (n = 3) Carpathians. Preliminary results of microsatellite DNA analyses of a subsample from Poland were reported by Śmietana, Rutkowski, Ratkiewicz, and Buś-Kicman (2012). We also used DNA isolates obtained from brown bears sampled in Slovakia (n = 35) and Romania (n = 16) that were previously analyzed by Straka et al. (2012). In total, we used 452 samples: 44 from Western Carpathians (POL and SVK), 387 from the Bieszczady Mountains (POL and SVK), 16 from Romanian Carpathians, three samples from Ukrainian Carpathians, and two samples from Greece ( Figure 1). The last two samples were used in a mtDNA survey and their haplotypes were used in network analysis only.

| Laboratory procedures
Two classes of genetic markers were used: nuclear microsatellite Within analyses of the first set of samples, we chose and organized 12 autosomal and unlinked microsatellite markers (Bellemain & Taberlet, 2004;Paetkau, Calvert, Stirling, & Strobeck, 1995;Straka et al., 2012;Taberlet et al., 1997) into three multiplex sets that maximize the number of loci suitable for simultaneous amplification. Set 1 consisted of 5 loci: G10M, G10J, Mu09, Mu61, and Cxx20; Set 2 consisted of 4 loci: G10B, G10C, Mu51, and Mu59; Set 3 amplified 3 loci: Mu10, Mu11, and G10X. In addition, we amplified a ~76 bp fragment of the Sry gene of males (Straka et al., 2012;Taberlet et al., 1997) together with microsatellites in Set 2 in order to genetically identify the sex of the individuals studied. Multiplex PCRs were performed with GeneAmp PCR System 9,700 (Applied Biosystems) in 10 μl reaction volume containing 2 μl of isolated genomic DNA, 4.5 μl Qiagen Multiplex PCR Master Mix (1×), 0.9 μl mix of primers (0.2 μM of each primer), and 2.6 μl RNase-free water. Each multiplex PCR started with an initial activation step at 95°C for 15 min, followed by 42 cycles, with denaturation at 94°C for 30 s, annealing for 90 s (at 49.6°C for Set 1, 58°C for Set 2, and 54°C for Set 3), extension at 72°C for 60 s, and final extension at 60°C for 30 min.
The PCR products were mixed with 10 μl ultragrade formamide and 0.2 μl GeneScan 500-LIZ size standard (Applied Biosystems), denatured at 95°C for 5 min, rapidly cooled and detected using ABI 3,130 Genetic Analyzer (Applied Biosystems). The allele fragment lengths were estimated using the Auto Bin feature in GeneMapper 4.0 software (Applied Biosystems). Four microsatellite loci (Cxx20, Mu10, Mu51, and G10X) produced unreliable results for stool samples (no amplification, irregular stuttering, or artefacts affecting scoring) and therefore were discarded from multiplex sets after initial screening.
PCR reactions were performed from two to four times depending on the consistency of obtained PCR products. During initial screening, samples with <4 genotyped loci were removed and not included in any further analysis. For samples which passed this step but did not produce reliable genotypes at one up to four loci, we attempted to fill in missing or weak data by performing additional PCRs (in multiplex sets or targeting single loci). All pairs of remaining unique geno- Permanent presence is indicated with orange, while sporadic occurrence with yellow. Individuals sampled in Romania and Ukraine are indicated with black dots estimated using GIMLeT 1.3.3 (Valière, 2002). For the remaining eight microsatellite loci, the frequencies of null alleles were low to moderate (Supporting Information Table S1) and thus these loci were used for further analyses. The allelic dropout rate was estimated at 0.049 across loci and depended on allele size of the particular loci (Supporting Information Table S1).

The second set of samples was analyzed by Wildlife Genetics
International (WGI) with the use of the same set of eight selected microsatellite markers, and ZFX/ZFY marker for sex identification (instead of Sry used in Laboratory of Molecular Biology). All microsatellite loci analyzed by WGI were amplified separately. Genotype The reaction conditions followed the same protocol as used for the microsatellites, it ran for 40 cycles with annealing at 57°C.
Sequencing reactions in both directions were performed using the BigDye™ Terminator Cycle Sequencing Kit (Applied Biosystems).
The reaction conditions were as follows: 25 cycles with denaturation at 95°C for 20 s, annealing at 50°C for 15 s, extension at 76°C for 60 s. The detection of sequencing reaction products was carried out on ABI 3,130 Genetic Analyzer. Sequences were aligned manually in the BIoedIT sequence editing program (Hall, 1999

| Phylogenetic analyses
To test phylogenetic relationships among concatenated cytochrome b and control region haplotypes (809-813 bp long), we constructed a neighbor-joining tree with MeGa v.5.05 (Tamura et al., 2011) with 1,000 bootstrap replicates used to assess support for tree nodes. All indels found within mtDNA control region sequences were excluded from phylogenetic and population genetics analyses. We also used brown bear cyt b and mtDNA control region sequences deposited in GenBank (from Bon et al., 2008;Keis et al., 2013;Miller et al., 2012;Taberlet & Bouvet, 1994). Ursus americanus (AF303109) was used as an outgroup. Haplotype network reconstruction was performed in neTwork v4.6.1.0 (Bandelt, Forster, & Röhl, 1999). We also calculated net pairwise divergence (Da) among mtDNA lineages using MeGa v.5.05.

| Population genetics
We used Cervus 3.0.3 (Kalinowski et al., 2007) and For mtDNA analysis, we calculated the number of haplotypes (Nh), haplotype diversity (h), nucleotide diversity (π), number of segregating sites (S), and mean number of pairwise differences (PD) for concatenated mtDNA control region and cyt b using software packages arLequIn (Excoffier & Lischer, 2010) and dnasp v.5 (Librado & Rozas, 2009 Table S2). There was a full consistency of sex identification using Sry and ZFX/ZFY markers for the samples studied in both labs.
For the total sample, the number of alleles per locus (A) ranged from 6 (at 3 loci) to 13 alleles (Mu59); on average 8.63 alleles (Supporting Information Table S1). For the eight microsatellite loci studied, the allele numbers and their size ranges were very similar in three studied brown bear populations (Supporting Information Table   S1), implying genetic affinity among them. The gene diversity at the studied loci ranged from 0.722 to 0.795 and F IS values in brown bear populations were not significantly different from zero (

| mtDNA analyses
The analysis of concatenated mtDNA sequences (809-813 bp) yielded nine haplotypes and the tenth haplotype (H10) was found in Greece (Figure 3a,b). The phylogenetic analysis ( Figure 3) revealed that all mtDNA haplotypes detected in this study belonged to either the Eastern lineage of brown bears (Taberlet & Bouvet, 1994, e.g Table   S3). Molecular diversity indices for the brown bear populations studied and for the whole sample at mtDNA are given in Table 3.   This suggests that the observed admixture level cannot be explained by short-term gene flow between these populations.

| D ISCUSS I ON
Habitat fragmentation is assumed to be a general factor that may increase female philopatry (Henry, Coulon, & Travis, 2016). If we consider that habitat discontinuity and physical obstacles to gene flow might cause genetic differentiation between brown bears even at a relatively small distance in just a few generations of isolation (Straka et al., 2012;Straka, Paule, Štofík, Ionescu, & Adamec, 2011 (Štofík, Bural, Paule, & Straka, 2010), possibly resulting in a genetic admixture. Straka et al. (2012) supposed that migration between Western Carpathians and Bieszczady Mountains was excluded.
However, Paule et al. (2015) identified two males, which migrated from Western Carpathians to the Bieszczady Mountains via straightline (distance about 180-200 km). It is in agreement with the sex-biased dispersal and asymmetric admixture found in our study.

Brown bears in the Bieszczady Mountains (Poland and Slovakia)
after World War II have been a minority (Jakubiec, 2001;Sabadoš & Šimiak, 1981 (Ziółkowska et al., 2016) should be performed, migratory corridors should be identified and incorporated in regional plans of infrastructure and housing development. Fernández et al. (2012), for example, highlighted the need to control unplanned urban sprawl to preserve both brown bear habitats and the connectivity between the Western and Eastern Carpathian populations. To conclude, our genetic survey highlights a pivotal role of transboundary migratory corridors to ensure the long-lasting existence of brown bear populations in Europe.

ACK N OWLED G M ENTS
We would like to thank F. Zięba, K. Bojarska, Z. Wojtas, Regional Directorate of Environmental Protection in Rzeszów, and State Forest Administration in Krosno (Poland) and M. Rosada (NNP "Synevyr," Ukraine) for their help in retrieving samples. We also would like to thank S. Snively for his laboratory work, R. Rutkowski, M. Buś-Kicman, and D. Pietrewicz-Kubicz for their help in DNA isolation from portion of samples, and P. Rode for his help during figure preparation. We thank also D. Paetkau for his insightful comments and for checking English grammar. We also would like to thank C. Abarca and C. Obiesie who line-edited the manuscript for submission. This work was financed by WWF Poland, EEA grants, and Norway grants (project PL 0349) and is also supported by the University of Bialystok (BST-102).

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.