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 (Babik et al. 2005; 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 (Ionescu et al. 2010) 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–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.

Materials and methods

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.

Fig. 1
figure 1

GPS locations of the genetic samples: males are represented in blue, and females are represented in red

Full size image

DNA extraction, amplification, and sizing

DNA was extracted from hair samples using the Maxwell® 16 Tissue D.N.A. Purification Kit (Promega, USA) (Cotovelea et al. 2015) and optimized according to Fedorca et al. (2018). Nine fluorescently labelled microsatellite markers (Crawford et al. 2009; Pelz-Serrano et al. 2009; Frosch et al. 2010; Minnig et al. 2016) were used and divided into two multiplexes: Multiplex I (Cca13, CF33, CF07, CF06) and Multiplex II (Cca18, CF31, CF05, CF44, CF19). The PCR reactions were performed in a 15 μl total volume, containing 7.5 μl of Qiagen Multiplex PCR mix, 50–100 ng/μl DNA and fluorescently labelled primers (0.2 and 0.4 μM forward and reverse primers, respectively). Fragment analysis was performed using the GenomeLab™ GeXP Genetic Analysis System (Beckman Coulter), and alleles were scored using GenomeLab™ Software (Beckman Coulter, Inc). Furthermore, we proceeded on amplifying the non-labelled SRY primer (forward: 5′-TGA ACG CAT TCA TGG TGT GGT-3′; reverse: 5′-AAT CTC TGT GCC TCC TGG AA-3′) in 5 μl total volume; the amplification was carried out using Qiagen Multiplex PCR and 50–100 ng/μl DNA after pre-denaturation (95 °C for 15 min) in 35 cycles (94 °C denaturation, 60 °C annealing, and 72 °C elongation; 45 s for each step) and final extension at 60 °C for 30 min.

Genetic data analysis

The microsatellite data were examined for the presence of null alleles and microsatellite stutter bands using MicroChecker (Van Oosterhout et al. 2004). Standard population genetic diversity indices (number of different alleles (Na), number of effective alleles (Ne), expected heterozygosity (He), observed heterozygosity (Ho), 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.

Results

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 (He = 0.601) was higher than those previously reported in Europe (Frosch et al. 2010; Senn et al. 2014; Minnig et al. 2016) (Table 1).

Table 1 Eurasian populations and comparable measurements of genetic diversity
Full size table

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.

Table 2 Genetic diversity within beaver population in the study
Full size table

N sample size, Na numbers of alleles, Ne the effective number of alleles, Ho observed heterozygosity, He expected heterozygosity, and F fixation index.

Population structuring

We inferred no presence of genetic clusters when using Bayesian clustering, employed in STRUCTURE, and plotted in CLUMPP (Fig. 2; k = 2, 3, 4, 5, 6, 7). Structure Harvester indicated the presence of 6 genetic groups (Dk) and one genetic group when using LnPr (X|K) method (Fig. 3).

Fig. 2
figure 2

Maximum likelihood analyses of the population structure resulted from STRUCTURE and plotted in CLUMPP for k = 2, 3, 4, 5, 6, 7 with the vertical bars representing individuals and their respective cluster memberships

Full size image
Fig. 3
figure 3

Distribution of ΔK and the mean of the estimated Ln probability for the data, as determined using Structure Harvester

Full size image

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.

Fig. 4
figure 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)

Full size image

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.

Fig. 5
figure 5

Principal coordinate analysis (PCoA) of all samples

Full size image

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

Fig. 6
figure 6

Mantel test results for the entire database a), females b) and males c)

Full size image

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

Table 3 Results (P-values) of bottleneck testing under three models
Full size table

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

Analyses with BOTTLENECK used three microsatellite mutation models: infinite allele model (IAM), two-phase model (TPM), and stepwise mutation model (SMM), statistically significant (p < 0.05).

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 et al. 2020) originated from one female and two males and went through a severe population bottleneck, these could not be compared with our case, where we have reintroduced a large number of individuals of mixed origins (Russia, Norway, and France) from various geographic areas in Bavaria.

To our knowledge, we provide the first population genetic assessment on Eurasian beaver in Romania. Three key 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 (Frosch et al. 2014; Minnig et al. 2016), Finland, Russia, Estonia, and Lithuania (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.