Genetic structure of Eurasian beaver in Romania: insights after two decades from the reintroduction

Once exploited for fur, meat, and extracting the yellowish exudate called castoreum, the Eurasian beaver disappeared from Romania during the eighteenth century. After, the reintroductions carried out two decades ago are currently thriving in the Danube River basin. Using nine nSSR markers, we analysed samples from 98 individuals, and we found no genetic substructure, suggesting high dispersal and gene flow capabilities. The stepwise mutation model (SMM) indicated the existence of a recent genetic bottleneck, though the Eurasian beaver retains high levels of genetic diversity and population growth facilitated variation in nSSR loci. A fine-scale spatial correlation in females was detected, contrasting with males’ dispersal on longer distances. While the movement and establishment of individuals’ new territories were made under natural predation pressure, the mix following natural expansion improved the fitness and could contribute to a higher genetic diversity than the source population. With any reintroduction, a focus on capturing individuals from various geographic origins, as well as securing many and suitable founding individuals (adults, subadults, and juveniles) with mixed origins, could secure the post-genetic bottleneck recovery and higher genetic diversity. Beyond this conservation success, future management strategies should consider building a National Action Plan (NAP) for the species, including a permanent genetic monitoring programme for Eurasian beaver.


Introduction
Eurasian beaver (Castor fiber Linnaeus, 1758) is commonly considered an ecosystem engineer able to shape the physical environment and associated river corridors' hydrologic, geomorphic, biogeochemical, and ecological processes (Wright et al. 2002;Rosell et al. 2005;Wright 2009). Beaver reintroductions are increasingly regarded as an integral component of the broader restoration of river corridors (Denise et al. 2014;Law et al. 2017). However, species activities influence geomorphic heterogeneity and retention of water, sediment, and organic matter influx along river corridors: building dams, digging narrow canals to facilitate their movements, and altering riparian vegetation via herbivory (Laurel and Wohl 2019).
It is well known that the main threat to the beaver populations is in general the low level of genetic diversity Durka et al. 2005), although a limited number of nSSR genetic markers have been developed for beaver population genetics (McEwing et al. 2015). However, small, isolated populations might be exhibiting record losses in genetic diversity, which result in even more significant reductions in their ability to cope with environmental changes (Richard et al. 2017). Besides, monitoring the reintroduced populations is considered crucial (Senn et al. 2014), and the simultaneous optimization of monitoring strategies and the application of adaptive management of the habitat quality is a necessity (Breitenmoser et al. 2001;Andrew et al. 2012).
The refugial populations (such as Bavaria) were strongly bottlenecked and have registered through genetic drift, a loss of genetic diversity (Halley et al. 2020). As the Switzerland example confirmed, the reintroduced Eurasian beaver populations are able to firmly establish, despite a shallow level of genetic diversity, and are rapidly expanding, and inbreeding depression is low (Minnig et al. 2016). Bayesian clustering methods are essential for the assignment of populations to genetically distinct groups for wildlife management and conservation, and these methods contribute to the identification of small-scale substructures (Pearse and Crandall 2004;Crawford et al. 2009;Hobbs et al. 2011).
The provenance of beavers that were reintroduced in Europe has not always been regarded as a significant factor in decision-making, and some mixed genetic forms of beavers have been released in the west (e.g. Germany, Switzerland) and in the east (e.g. Poland, Romania) (Gaywood 2018). However, successful reintroductions generally increase a species' population size and geographic range and restore ecological function to the area where it was extirpated (Doug et al. 2015).
The Eurasian beaver disappeared from Romania at the beginning of the eighteenth century  due to uncontrolled hunting for fur and meat and extracting the yellowish exudate called castoreum. Both sexes use castoreum, a liquid, viscous secretion secreted by the perianal glands to mark their territories, and olfactory signals are essential in identifying individuals and families (Rosell and Schulte 2004). In recent years (1998)(1999)(2000)(2001)(2002)(2003), within a European initiative, 182 individual Eurasian beavers captured in Bavaria (Germany) from various areas (Sinnbronn, Weißenbrunn, Radlig, Wilburgstetten, Mühlhausen, Großaitingen, Marchaney, and so on) were reintroduced into three main river basins: Olt (91 individuals), Mures (56 individuals), and Ialomita (35 individuals). Public opinion consultation, determination of favourable zones, and ecological diagnosis have preceded the reintroductions; nevertheless, no genetic analyses were performed for the reintroduced Eurasian beavers. The 182 beaver individuals were mostly adults; however, there have also been released subadults and juveniles. Currently, over 2250 Eurasian beaver individuals occupy habitats across the country, including in tributary rivers, due to the natural dispersal of beavers occurring first over short distances (Pasca et al. 2016). The CLMAN project report (2015) noted that the Olt River shelters more than 1300 individuals, the Mures River shelters approximately 320 beavers, and nearly 220 individuals were found on the Ialomita River. New presence points were reported in new river basins, such as the Viseu (Maramures County), Somes (Cluj County), Timis (Timis County), Trotus, Uz (Bacau County), Crasna (Buzau County), and Danube Delta.
This study addresses the first time after Eurasian beaver was reintroduced in Romania gaps in knowledge relevant for species conservation: (i) we assessed the genetic diversity of the reintroduced beaver population; (ii) we evaluated the population structuring and spatial autocorrelation, and we tested for the presence of a genetic bottleneck.

Study site and sampling
All samples originate from the following river basins: the Olt, Ialomita, and Mures. Ninety-eight samples were collected within special derogations issued by the Romanian Academy of Science (2016-2017) for conducting scientific studies (Fig. 1). Eurasian beaver individuals were captured using Bailey traps, positioned on movement paths, between 3 to 5 m away from water, well anchored in the ground. Eurasian beavers were tranquillized, and during handling and measuring procedures, we collected hair samples stored in a paper envelope for each individual and kept them at room temperature.

Genetic data analysis
The microsatellite data were examined for the presence of null alleles and microsatellite stutter bands using Micro-Checker (Van Oosterhout et al. 2004). Standard population genetic diversity indices (number of different alleles (N a ), number of effective alleles (N e ), expected heterozygosity (H e ), observed heterozygosity (H o ), and fixation index (F)) were calculated using GenAlEx v.6.5. Furthermore, we used a Bayesian clustering method in Structure v.2.3.4 (Pritchard et al. 2000;Falush et al. 2003) to infer the population structure based on 1,000,000 simulations preceded by 100,000 burn-in replicates conducted using an admixture model and correlated allelic frequencies. The K values varied between 1 and 10; we performed ten repetitions for each K value. However, we have used CLUMPP to build the complied histograms (Jakobsson and Rosenberg 2007). Both the posterior probability of the data with a given value of K (Ln Pr(X|K)) and its rate of change (ΔK) were used to analyse the population structure (Evanno et al. 2005). We estimated the pairwise Euclidean geographic distances between all sampling points, which allowed us to test for the presence of simple isolation by distance (IBD) (Wright 1943), assuming that genetic differentiation is a by-product of geographic distance without taking any landscape features into account. The IBD was calculated in GenAlEx v. 6.5 using 10,000 random permutations. Correlograms showing the pairwise genetic relatedness for the whole dataset and separately for females and males were computed using 9999 permutation and even distance classes. BOTTLENECK V1.2.02 (Cornuet and Luikart 1996;Cornuet 1999;Piry et al. 1999) was used to validate whether the population had undergone through population bottleneck.

Genetic diversity
The number of different alleles ranged from 3 (Cca18) to 7 (CF33) at each locus, with similar results for the number of effective alleles ( Table 2). The mean genetic diversity in the Romanian beaver population (H e = 0.601) was higher than those previously reported in Europe (Frosch et al. 2010;Senn et al. 2014;Minnig et al. 2016) (Table 1).
For all loci, the observed heterozygosity was higher than the expected heterozygosity (Table 2), while the fixation index included negative values. However, this finding may indicate the existence of selection in favour of heterozygotes, which may explain the observation of some types of genetic variability.

Spatial autocorrelation and principal coordinate analysis (PCoA)
Beavers (both sexes) showed significant spatial autocorrelation coefficient (Fig. 4a), r, up to 20 km (r = 0.073, P = 0.001); similar results were obtained for females (r = 0.813, P = 0.009) (Fig. 4b), whereas for males, we detected no significant spatial autocorrelation at any of the distance classes (Fig. 4c). However, the results could suggest that beavers within these distance classes could be related to each other as first cousins.
The principal component analysis (PCoA) offered support for one genetic group (Fig. 5a) when applied on genetic distances matrix, a similar result with the mean of the estimated Ln probability from Structure Harvester. In contrast, the PCoA applied on the geographic distances matrix showed a potential grouping of the samples in four groups (Fig. 5b). The grouping from PCoA on geographic distances matrix represents subdivisions of the basin rivers: (1) Ialomita, (2) Covasna-Brasov, (3) Fagaras-Sibiu-Podisul Hartibaciu-Mures, and (4) Mures-Arad.
Moreover, when analyzing the relationship between geographic distance and genetic divergence using the Mantel test, the results were not statistically significant neither for both sexes (Fig. 6a), females (Fig. 6b) or males (Fig. 6c).
The similarities among individuals with geographic proximity confirm the presence of gene flow, which represents a hypothesis supported by the Mantel test results and the spatial autocorrelation, suggesting that no genetic isolation occurred over time.

Population bottleneck
We used BOTTLENECK to test for population bottleneck in this study, and we applied the infinite allele model (IAM), two-phase model (TPM), and stepwise mutation model (SMM) ( Table 3). While the infinite alleles model (IAM) could be inappropriate for microsatellites due to high microsatellite mutation rates and mutation processes that might retain the memory of ancestral allelic states (Goldstein et al. 1995), we have considered the results only for TPM and SMM. Nevertheless, nSSR loci may closely follow the SMM model, as is the most conservative model for testing for a significant heterozygosity excess caused by genetic bottlenecks (Cornuet and Luikart 1996), yet TPM is more appropriate for nSSR loci (Rienzo et al. 1994;Piry et al. 1999). The sign test indicated the presence of a genetic bottleneck effect for IAM and SMM (p < 0.05), while TPM was not significant (p = 0.200). The sign test, standardized differences test, and Wilcoxon test showed a significant heterozygosity excess under the IAM (p < 0.05). Simultaneously, TPM was significant only for the standardized differences test and Wilcoxon test, and SMM was significant for the sign test and standardized differences test.
The Wilcoxon test results, considered to be more reliable than the sign test and standardized differences test, showed no significant results for population bottleneck under the SMM but has indicated the presence of genetic bottleneck effect for IAM (p = 0.00195) and TPM (p = 0.01367).

Discussion
Beavers are ecosystem engineers whose presence has multiple effects on riparian habitats (Halley et al. 2020); however, the low genetic variability among populations represents a threat in a changing environment, especially if the number of founders is low (Iso-Touru et al. 2020). Whereas the Castor fiber population in Finland (Iso-Touru  To our knowledge, we provide the first population genetic assessment on Eurasian beaver in Romania. Three key Fig. 4 Spatial autocorrelation correlograms of the entire sample a), for females b), for males c). The pairwise genetic relatedness for nine distance classes is ranging from 1 to 500 km. r is the coefficient of similarity among the individual genotypes within a distance class (solid line). Bootstrap error bars are also shown for each distance class. The red lines (U and L) represent the upper and lower confidence intervals for the null hypothesis (no autocorrelation among the genotypes within the distance classes (r = 0) European Journal of Wildlife Research (2021) 67: 104 insights emerge from this work. First is that the Eurasian beaver population in Romania has a higher degree of genetic diversity than other populations in Central Europe Minnig et al. 2016), Finland, Russia, Estonia, andLithuania (Iso-Touru et al. 2020), suggesting that species is firmly established ( Table 1, Table 2). The founding population itself in Bavaria has a high level of genetic diversity, a likely consequence of an admixed descent (Campbell-Palmer et al. 2020). However, in our case, individuals/groups were repeatedly removed under the pressure of natural and human-introduced predators (such as brown bears and stray dogs; Anonymous 2021), and mixing following natural expansion could contribute to a higher genetic diversity than the source population, improving the fitness.
Overall, model-based clustering inferred one Eurasian beaver population in Romania (Fig. 2, Fig. 3). The lack of genetic sub-structuring (Fig. 5a, b) suggests high dispersal and gene flow capabilities and implies that suitable habitats facilitated the movement and establishment of individuals' new territories. Habitat selection in recently documented beaver recolonizations is affected by different scale-related phenomena such as food availability (woody plants), densities, and vegetation cover, or could be an independent selection at the largest scale (between watersheds) (Zwolicki et al. 2018).
The third insight from our study is the pattern detected of fine-scale spatial correlation in females (Fig. 4b) which could be attributed to the benefit from easy access to resources, an adaptive strategy depending on population density (Mayer et al. 2017b). However, senescing individuals become more tolerant towards philopatric offspring, allowing them in their natal territory to help in protecting boundaries (Ekman et al. 2001). The takeover of the breeding position could benefit from remaining in the natal territory (Mayer et al. 2017a). In contrast, the Eurasian beaver's male's dispersal was documented on longer distances (Saveljev et al. 2002;Mayer et al. 2017b), which was confirmed also here (Fig. 4c). These patterns could have significantly impacted dispersal after reintroducing Eurasian beavers, gene flow, and the ability to colonize new and formerly occupied areas. The existence of similarities among individuals with geographic proximity confirms the presence of gene flow (Fig. 6).
Moreover, after reintroduction, population growth facilitated variation in nSSR loci in the Eurasian beaver population, even though passing the genetic bottleneck (Table 3), the population recovered post-genetic bottleneck and retained genetic diversity. However, dedicated studies using SNP and mitochondrial DNA and expanding the sample size to incorporate transborder areas that could be a permanent source of new genes are needed. While our results provide molecular evidence for the effectiveness of conservation strategies, future management strategies should consider building a National Action Plan (NAP) for the species, which includes a permanent monitoring programme for the Eurasian beaver. This NAP could incorporate digital models if restocking is needed, easily accessible for the decision-makers and stakeholders. A permanent genetic monitoring programme should be made available by NAP using microsatellites, SNPs, and mitochondrial DNA. Altogether, the Eurasian beaver represents a conservation success in Romania, and no restocking is needed at this time. However, if restoking will be considered for various reasons in the future, we recommend that the Eurasian beavers be located outside the area where the reintroduction nucleus was taken initially.
Foremost, the principles to be drawn from study results relating to mammalian reintroductions, in general, are to secure for relocation many individuals captured from various areas with source population large enough and with mixed origins, by using most adults, but also subadults and juveniles. Nevertheless, during reintroductions from mixed origins, populations could pass through genetic bottleneck due to founder effects associated with population subsampling and post-release mortality, even though securing many founding individuals could have ameliorated these effects in our case. Other recent studies on mammals (Gajdárová et al. 2021) have demonstrated the importance of protecting migratory corridors and landscape permeability and protecting the dispersing individuals along with securing suitable areas, which should be considered in every reintroduction. The main conclusion is that we ensured proper founder composition and high numbers without pre-release genetic screening, which was hardly accessible at the same time; however, they should represent vital steps in any contemporaneous reintroductions.
Author contribution A.F. and G.I. conceived of the presented work, G.I. ensured the funding, C.P., M.F., and A.G. gathered most samples. E.C. conducted the laboratory work and the primary analysis; A.F. analysed and interpreted the data and wrote the manuscript draft together with G.I. M.F. developed tables and figures. All the authors critically reviewed the manuscript and approved it.
Funding This work was supported by the projects PN16330207 and PN19070603 financed by the ANCSI (National Authority for Scientific Research) Nucleu Programme.
Data availability All relevant data are within the paper (Please, see tables and figures).
Code availability Not the case.