Modelling of competitive interactions between native Eurasian (Castor fiber) and alien North American (Castor сanadensis) beavers based on long-term monitoring data (1934–2015)

Modelling of competitive interactions between native Eurasian (Castor fiber) and alien North American (Castor сanadensis) beavers based on long-term monitoring data (1934–2015)

Ecological Modelling 409 (2019) 108763 Contents lists available at ScienceDirect Ecological Modelling journal homepage: www.elsevier.com/locate/ecol...

6MB Sizes 0 Downloads 0 Views

Ecological Modelling 409 (2019) 108763

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

Modelling of competitive interactions between native Eurasian (Castor fiber) and alien North American (Castor сanadensis) beavers based on long-term monitoring data (1934–2015)

T

Varos G. Petrosyana, , Victor V. Golubkovb, Nikolay A. Zavyalovc, Lyudmila A. Khlyapa, Natalia N. Dergunovaa, Fedor A. Osipova ⁎

a b c

A.N. Severtsov Institute of Ecology and Evolution, 33 Leninskii Prospekt, Moscow, 119071, Russia Federal Research Center “Informatics and Management”, 44 Vavilov St., Moscow, 119333, Russia Rdeysky Nature Reserve, Kholm, Novgorod Oblast, 175271, Russia

ARTICLE INFO

ABSTRACT

Keywords: Competition population dynamics Native species invasion nature reserves

Our goal is to assess consequences of the introduction of alien North American (Castor öanadensis (Cc)) beaver into the Nature Reserves inhabited by Eurasian beaver (Castor fiber (Cf)) in European Russia using a mathematical model. For this reason, we have developed a two-species model of population dynamics. Long-term (1934–2015) monitoring data on Cf population dynamics in six Nature Reserves are used in computer modelling of competitive interactions between native Cf and alien Cö. The Reserves are located in the European part of Russia in the north, south, and central part of Cf range. We have simulated the dynamics of both species populations after the introduction of Cö into the habitats occupied by Cf. The model demonstrates that Cf is displaced by Cc in all the Reserves after the introduction of 2 – 24 individuals of Cö. However, the duration of exclusion of one species by the other varies as a function of ecological conditions, initial number of individuals, and fecundities. Our model shows that, in case of introduction of 12 Cö beavers, the size of Cf population starts to decrease after 31–146 years as a result of competition. We study the conditions providing the coexistence of both species and find that Cö population dynamics after Cf exclusion can be described by four patterns: irruptive (Lapland Reserve), single-stage (Prioksko-Terrasny Reserve), multi-stage (Darwin, Central-Forest, and Khoper Reserves) and logistic population growth (Oka Reserve). Species biology in terms of the fecundity, family size, rate of individual development until sexual maturity, life-span, and the age structure of populations are compared between species to detect the mechanisms providing the competitive advantage of Cö over Cf.

Introduction By the end of the 19th century, there were approximately 1,500 Eurasian beavers (Castor fiber L., 1758) (Cf) in 11 isolated populations in Eurasia (Portenko 1926, Grave 1931, Lavrov 1981, Nolet, Rosell 1998) (Table 1, Fig. 1). Subsequently, those (mainly European) populations served as a source for Cf reintroduction and restoration of its range. Besides, the North American beavers (Castor canadensis Kuhl, 1820) (Cc) were brought and released in the north-eastern parts of Poland in 1926 (Masuria), into the river Goryn’ (former Poland territory and now Ukraine) in 1933–1934, and in Finland in 1937 (Lathi, Helminen 1974, Ermala et al. 1989). The purpose of those introductions was to restore beaver’s populations in Europe because, at that time, many zoologists



considered the beavers of the Old and New Worlds as the same species. Only in 1973, Lavrov and Orlov (1973) showed that Eurasian beaver (C. fiber) had 48 chromosomes, whereas North American beaver (C. canadensis) 40 chromosomes (Lavrov, Orlov 1973), so Cc turned out an alien species for Eurasia. In Ukraine, Cс became extinct in 1957 and in Poland in 1979 (Zurowski 1980, Parker et al. 2012). In Finland, it has spread throughout its southern part of the country concentrating near the eastern boundaries, while it was rarely met in the north-west, near the borders with Norway (Kauhala, Turkia 2013). In the 50 s of the 20th century, Cс crossed the Russian-Finland border and invaded habitats in Karelia Republic and Karelian Isthmus (Zaikin 1959, Segal, Orlova 1961, Ivanov 1975, Monitoring and conservation … 2010). In Russia, alien Cс beavers were repeatedly caught and introduced to beaver-free

Corresponding author. E-mail address: [email protected] (V.G. Petrosyan).

https://doi.org/10.1016/j.ecolmodel.2019.108763 Received 26 March 2019; Received in revised form 12 July 2019; Accepted 22 July 2019 Available online 05 August 2019 0304-3800/ © 2019 Elsevier B.V. All rights reserved.

Ecological Modelling 409 (2019) 108763

300 30–40 ˜ 300

Serebrennikov 1929, Lavrov, Lavrov, 1986, Saveljev 2002 Lavrov 1969, Lavrov and Lavrov, 1986 Serebrennikov 1929, Lavrov, Lu 1961, Stubbe, Dawaa 1986, Stubbe et al. 1991

areas in Karelia and the Russian Far East from the middle of 1960s to the end of 1980s. Several cases of unintentional Cc introduction in Western Europe are known. In France, in 1977, three Cc beavers escaped from a private zoo and reached the Bourdon Reservoir on the Yonna River. According to Dewas et al. (2012), in 1984, the local Cс population consisted of 15–20 beavers. From the initial sites of introduction, the beavers were spreading and creating settlements in the basin of the Loire River. Because of the simultaneous presence of Cf and Cс in the Loire, it was decided to eradicate Cс. In 1984–1985, 24 North American beavers were captured and exterminated. In Austria (1978–1979 and 1986) and in Hungary (1990s), Cс also escaped from captivity. The beavers in Austria did likely not reproduce and soon became extinct (Sieber, Bauer 2001). The outcome of introduction in Hungary remained unclear (Parker et al. 2012). Cс was introduced to Germany several times and at different locations. In 1981 and 1989, Cс beavers were released in the North Rhine-Westphalia (NRW) which served as a source for further Cс distribution (Dewas et al. 2012). In 1994, Cс was released in Rheinland-Palatinate (RLP), which is located in the south-west of Germany, near the border with France, Luxemburg, and Belgium. In that region, the Cc population was represented both by reproductive beavers and by sterilized ones, which had been released to restrict the total population size (Dewas et al. 2012, Parker et al. 2012). A more tangled situation emerged on the territory of Bavaria, Germany, where beavers were introduced from Poland (Beaver farm in Popel’no) in 1966 (Parker et al. 2012), and Finland in 1972 (Lathi, Helminen 1974), i.e., before Cf and Cc were recognized as different species. Parker et al. (2012) indicated that Cс beavers could also be among those brought from the Polish beaver nursery, though the identification of 1,500 beavers did not reveal any. Another introduction of Cс from Finland was little discussed in literature sources, though there are data confirming the capture of naturalized Cс individuals in Finland and their transportation in 1972 for introduction in Bavaria (Lathi, Helminen 1974, Anonymous 1972). This fact is very important for understanding the species and subspecies structure of beaver populations in Bavaria, where, since 1970, individuals of different origins were repeatedly released (C. f. albicus, fiber, belorussicus, galliae, orientoeuropaeus). Later on, the beavers from Bavaria were regularly captured for distribution to the European countries: Hungary (2000–2002), Belgium (1998), Croatia (1996–2001), and Bulgaria (2001–2002) (Halley, Rosell 2002). Our analysis shows that Bavaria could be a hidden source of Cс in Central and Western Europe. The availability of the sources for Cc distribution in three European countries, namely, Germany, Luxemburg, and Belgium, was affirmed in 2006, when a dead Cc beaver was found at the border of Luxemburg and Germany (Dewas et al. 2012). Afterwards, Cc was met in four locations of Luxemburg (Herr, Schley 2009, Dewas et al. 2012). During 2009–2010 winter studies, 13 Cc beavers were removed from six locations in Luxemburg and Belgium, while eight beavers were withdrawn from three habitats in Germany (RLP) (Dewas et al. 2012). Currently, the largest populations of Cc in Eurasia are found in Karelia (Russia) and Finland. The Cc population in Karelia consists of 13,000 beavers (Monitoring and conservation … 2010). In Finland, there were 6,100 settlements, i.e., over 24,000 animals in 2013 (Brommer et al. 2017). The population size of Cf, based on monitoring counts in the autumn of 2017, was estimated at 3,300–4,500 animals, and Cc at 10,300‒19,100 (Beaver as a renewable resource 2019). At present, the two species are partly sympatric in Pirkanmaa and close to each other in Etelä-Pohjanmaa regions in Finland. The Cf is considered as ‘near threatened’ in the 2015 Finnish Red Book (Liukko et al. 2016). The issues of coexistence and competitive relations between these two species arise with regard to the recent findings of Cc in Central and Western Europe (Dewas et al. 2012). In Karelia and Finland, the creation of sympatric zones took from 67 to 80 years after the Cc introduction. The processes of mutual

9 10 11

Usman’ and Ivnitsa Rivers, the left tributary of the Voronezh River, the upper Don basin Rivers Konda, Bolshaya and Malaya Sos’va, tributaries of the Ob, Western Siberia River Azas, the upper Yenisei basin (Respublika Tyva) Bulgan River, the border of the south-western Mongolia and China (Xinjiang)

C. f. pohlei Serebrennikov, 1929 C. f. tuvinicus Lavrov, 1969 C. f. birulai Serebrennikov, 1929

Lavrov 1981, Stubbe, Romashov 1992

Lavrov 1981, 1983, Lavrov and Lavrov, 1986

˜ 400 (to the end of 1930s) 70 7

5 6

8

C. f. belorussicus Lavrov, 1981 C. f. belorussicus Lavrov, 1981

C. f. vistulanus Matschie, 1907 (?), C. f. belorussicus Lavrov, 1981 (?) C. f. orientoeuropaeus Lavrov, 1981

Grave 1931, Lavrov 1981 Portenko 1926, Grave 1931, Lavrov 1981

Halley, Rosell 2002 Halley, Rosell 2002 Heidecke, Hörig 1986 Fedushin 1935, Lavrov 1981

100 30 200 80-100 (in 1926) < 30 < 30 fiber L., 1758 galliae Geoffroy, 1803 albicus Matschie, 1907 belorussicus Lavrov, 1981 f. f. f. f. C. C. C. C.

Nid River, the south of Norway (Telemark) Delta of the Rona River (France) Elba River (Germany) Berezina River, the right tributary of the Dnepr River (Berezinskaya population), (Belarus) Sozh River, the left tributary of the Dnepr River (Sozhskaya population) (Belarus) Teterev River and the Pripyat’ River with its tributaries Ubort’ and Uzh (Dneprovskaya population) (north of Ukraine) Neman River basin (the border of Poland and Belarus) (Nemanskaya popluations), 1 2 3 4

Subspecies Location Ordinal number

Table 1 Characteristics of the Eurasian beaver isolated populations preserved in Eurasia by the end of the 19th century.

Population size, ind.

References

V.G. Petrosyan, et al.

2

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

Fig. 1. Geographical locations of 11 isolated Eurasian beaver populations at the 19th century (see Table 1 for the characteristics).

Russia by means of a mathematical model based on long-term monitoring data on beaver population dynamics. For this reason, we have developed a population dynamics model of the resource–consumer type with two species of competing consumers. This model has originated from the previously developed single-species model (Petrosyan et al. 2013), which was highly efficient for prediction of Cf population dynamics in Nature Reserves in the European part of Russia (Petrosyan et al. 2016).

penetration and exclusion of one beaver species by another continue in Karelian Isthmus (Leningrad Region) and in the north of Arkhangelsk Region (Monitoring and conservation … 2010). According to the monitoring studies in Karelia, some regions occupied by Cc over 20 years ago are currently inhabited by Cf and vice versa (Danilov, Fedorov 2015). Data on ecology and behavior of Cc and Cf in South Karelia showed that, in identical orographic, edaphic and hydrologic conditions, Cc and Cf beavers build lodges and dams with the same frequency (Danilov, Fedorov 2015). Parker et al. (2012) also indicated that differences in the life history, ecology, and behavior were insignificant and suggested that ecological niches of the two species completely overlap. However, the fecundity and family sizes are significantly greater in Cc than in Cf (Rosell, Parker 1995, Parker et al. 2012). Thus, there is no unified view on the competitive interactions between the two species, although they are in the focus of scientific attention (Monitoring and conservation … 2010, Parker et al. 2012, Danilov, Fedorov 2015, Dewas et al. 2012, Frosch et al. 2014). For example, Parker et al. (2012) analyzed the history of species restoration, their biology, and ecology, then suggested a strategic plan of Cc eradication from Eurasia. However, these publications do not give any specific estimates of the rate of competitive exclusion of one species by the other. Some researchers considered Cf as inferior relative to Cc (Müller-Schwarze, Sun 2003), but according to the monitoring data in Karelia, Cf has an advantage over Cc at the initial stages of competition (Danilov et al. 2007, p. 54). Cc was long believed to be a superior competitor as it achieves sexual maturity faster, has greater fecundity, and its building activity is more intense (Müller-Schwarze, Sun 2003). However, data on ecology and behavior of two beaver species in the same conditions of the south Karelia showed that Cc and Cf built the lodges and dams with the same frequency (Danilov, Fedorov 2015). Consequently, the rate of population growth is the only competitive advantage of Cc. So, it is not clear whether two beaver species would coexist on the same territory or one of them would be eventually excluded by the competitor. Therefore, our goal is to assess the consequences of the Cc introduction into Nature Reserves inhabited by Cf in the European

Materials and Methods Long-term monitoring data To establish demographic and resource parameters characterizing the state of beaver populations, we use the long-term monitoring data on Cf population dynamics in six Nature Reserves: Lapland Nature Reserve (LNR), Darwin Nature Reserve (DNR), Central-Forest Nature Reserve (CFNR), Prioksko-Terrasny Nature Reserve (PTNR), Oka Nature Reserve (ONR), and Khoper Nature Reserve (KNR) (Petrosyan et al. 2016). These Nature Reserves are located in European Russia in the northern, southern and central parts of the Cf range (Fig. 2, Table 2). Cf was introduced twice in the LNR, in 1934 and 1937 (14 beavers in total) (Kataev 2015); in the DNR in 1976 (Zavyalov et al. 2005); in CFNR in 1936 and 1937 (four couples) (Korablev et al. 2011). In 1948, after an almost 400-year absence from the territory of the PTNR (Albov, Khlyap 2015), two couples of Cf beavers were released into the Tadenka river. In 1937–1940, 23 beavers of Cf were released in the ONR (Kudryashov 1975). In 1937–1939, 23 Cf beavers were released in KNR (D’yakov 1975). Demographic and resource parameters of the model were derived from the long-term monitoring data on the beaver population dynamics for LNR from 1934 to 1985, for DNR from 1980 to 2001, for CFNR from 1936 to 2010, for PTNR from 1948 to 2015, for ONR from 1949 to 2010, and for KNR from 1957 to 2012 (Petrosyan et al. 2016).

3

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

Fig. 2. Locations of six studied Nature Reserves (shown as black spots). The names and boundaries of vegetation subzones in the territory of Russia and neighboring states (the former USSR) follow Ogureeva et al. (1999). Latin letters indicate the subzones populated by beavers: a, the hypoarctic (taiga) type of altitude zonation; b, northern taiga subzone; c, medial taiga subzone; d, southern taiga subzone; e, subtaiga subzone; f, broad-leaved forest subzone; g, forest-steppe subzone; h, northern steppe subzone.

The model

Potential resource is a resource that can be restored from degraded resource and turn into the active form. Degraded resource cannot be used by beavers; it is a product of consumption process by beavers. In nature, the state of these resources is dependent on the habitat conditions: interchanges of periods with beavers and those without them (Logofet et al. 2015, Zavyalov et al. 2016).

The two-species model is built on the following assumptions. Assumption 1. During the study period, the common resource R (k ) does not change with time k, i.e., the favorable common resource of two species is limited and characterized by R: R (k ) = R = const . It means that potential changes of common resources caused by other reasons are significantly less than the impact of beavers.

Assumption 3. In the model, the rate, F (x ) , of beaver population increase depends on the amount, x, of active resource available for one individual, x = Rk(a) Pk , where Pk = Pk(1) + Pk(2) , Pk is the total of two beaver population sizes, and Pk(1) , Pk(2) are population sizes of species 1 (native) and species 2 (alien) at time k, respectively. The rate F (x ) is constrained in such a way that the function 1 + F(x) (the population growth rate) is bounded, continuous and monotone increasing for all x ≥ 0, guaranteeing zero or positive values. Fig. 3 shows a qualitative pattern of the beaver population increase coefficient F (x ) , as a function of x.

Assumption 2. According to published data (Gurney, Lawton 1996, Wright et al. 2004, Petrosyan et al. 2013), the common resource comprises three components Rk(a) , Rk(p) , and Rk(d) ; it can be represented as follows: R = Rk(a) + Rk(p) + Rk(d). Where, Rk(a) , Rk(p) and Rk(d) are active, potential and degraded resources, respectively. Active resource is a resource available for consumption. 4

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

Table 2 General description of the six Nature Reserves. Nature Reserve

Location, with geographical coordinates (Fig. 2)

Description of the land cover, references

Lapland (LNR)

Mountain areas of the Kola Peninsula to the north of the Arctic Circle, 67°39'-68°15′N, 31°44'-32°44′E

Darwin (DNR)

Shore of the Rybinsk Reservoir, 58°28'-58°50'N, 37°29'-38°10'E

Central-Forest (CFNR)

Center of the European part of Russia, 56°26'-56°39'N, 32°39'-33°01'E

Prioksko-Terrasny (PTNR)

Left side of the Oka River, 54°51'-54°55'N, 37°34'37°41'E

Oka (ONR)

Left side of the Oka River, 54°40'-55°0'N, 40°35'-41°0'E

Khoper (KNR)

Valley of the Khoper River, 51°14'-51°24'N, 41°58'-41°94'E

The hypoarctic (taiga) type of altitude zonation (Ogureeva et al. 1999). Forests, mountain tundra, lakes and rivers, and wetlands occupy 61%, 28%, 3%, and 8% of the Reserve’s territory, respectively. There are 8 lake-river systems in the Reserve, with several rivers flowing from the source to the mouth through the Reserve territory (Beavers in the Reserves… 2018). The southern taiga subzone (Ogureeva et al. 1999). Pine (Pinus sylvestris L.) forests occupy 73.5% of the area, birch (Betula pendula Roth), spruce (Picea abies (L.) Karst.), aspen (Populus tremula L.) and black alder (Alnus glutinosa (L.) Gaertn.) occupy 19.6%, 5.2%, 1.3%, and 0.4%, respectively. The area of temporary flooding of the Rybinsky Reservoir totals to 20,000 ha. Wetlands occupy 27% of the Reserve area (Zavyalov 2015). The sub-taiga subzone (Ogureeva et al. 1999). The territory serves as a watershed of two large rivers, the Volga and the Western Dvina. Spruce forests occupy 46% of the area, birch and aspen forests 40%, pine forests 9%, black alder forests 1%. About 4% of the area is occupied by the peat bogs (Korablev et al. 2011; Zavyalov 2015). The ecotone area between sub-taiga and broad-leaved forests (Ogureeva et al. 1999). Almost the entire territory is covered with forests. Pine (40%) and birch (35%) forests prevail. Among other abundant forest species, there are aspen, spruce, oak (Quercus robur L.), linden (Tilia cordata Mill.), and black alder. The hydrological network is characterized by small rivers, lakes, and wetlands (Beavers in the Reserves… 2018). The subtaiga subzone (Ogureeva et al. 1999). Forests mostly consist of pines with birches mixed with broad-leaved tree species. Hydrological network involves the Oka River with tributaries, numerous oxbows, and lakes. Large wetland areas occupied the central and northwestern parts of the Reserve (Beavers in the Reserves… 2018). The forest-steppe subzone (Ogureeva et al. 1999). About 80% of the area is covered mostly with floodplains and upland oak forests. Its smaller part is occupied by steppes and meadows. During the high water, about 80% of the area is usually flooded. There are about 400 lakes and oxbows in the Reserve (Beavers in the Reserves… 2018).

of model parameters q1, q2 , q3, q4 (Petrosyan et al. 2013). Model (1) includes the following parameters: q1 is the maximal value of the population increase rate when the relative active resource tends to 1; q2 is the population increase rate when the level of active resource per individual coincides with the optimal level under particular environmental conditions, where always q2 < q1; q3 is the threshold of per capita active resources at which the pattern of population increase rate alters; q4 is the parameter determining the steepness of function F (x , q1, q2, q3, q4 ) at the inflection point (the value of the derivative w.r.t. x), which is reached at x = q3 (the value of the derivative at this point equals (q1 q2) q4 ) ; q5 is the fraction of the potential resource that transforms into the active form for a year, i.е., the intensity of resource restoration; q6 is the fraction of the active resource utilized by an individual and thus transformed into the degraded recourse; q7 is the fraction of the degraded resource transforming into the potential resource for a year. All the above parameters of the model can take only positive values except for the parameter q4 , due to their physical meanings. Further, we can assume that the parameters q4(i) (i = 1, 2) take nonnegative values since function F(x , q1(i) , q2(i) , q3(i) , q4(i) ) depends on them. However, parameter q4(i) cannot be equal to 0. Indeed, according to Assumption 3, function F (x , q1(i) , q2(i) , q3(i) , q4(i) ) should be continuous for all x, while it displays a jump discontinuity at x = q3(i) when q4(i) = 0 . Thus, we can consider all the parameters of system (1) to take positive values. Also, parameters, q5 , q6 , q7 , of the beaver environment should not exceed 1 as fractions. To have F (x , q1(i) , q2(i) , q3(i) , q4(i) ) monotone increasing, it is necessary and sufficient to meet the following condition: q1(i) > q2(i ) . According to Assumption 3, the following conditions should be met:

Fig. 3. Form of the population rate of increase, F (x ) , as a function of the active resource amount per individual, x = R (a) P .

Assumption 4. Since the absolute values of resources are commonly difficult to obtain in the field, the proposed model uses the relative values measured as fractions of the total resource. It means that R = 1. Ensuing from Assumptions 1–4, the parametric discrete-time model describing the development of competitive relations between two beaver species takes on the following form: (1) (a ) Pk(1) Pk , Q(1) ) Pk(1) + 1 = Pk + F (Rk (2) (a ) Pk(2) Pk , Q(2) ) Pk(2) + 1 = Pk + F (Rk

q6(2) Pk(2) + q5 Rk(p)

Rk(a+)1 = Rk(a)

q6(1) Pk(1)

Rk(p+)1 = Rk(p)

q5 Rk(p) + q7 Rk(d) , Rk(d) = 1

k = 0, 1, 2, ..., F (x , Q) = q2 + (q1

q2)

x q42

+ (x

q3 q3)2

Rk(a)

Rk(p)

,

Q(1) = (q1(1) , q2(1) , q3(1) , q4(1) ), Q(2) = (q1(2) , q2(2) , q3(2) , q4(2) ), Q = (q1, q2 , q3, q4).

(1)

Here, and Rk(d) are the numbers of beavers of the 1 st (Cf) and 2nd (Cc) species, the levels of active, potential, and degraded resources at the kth year, respectively. F (x , Q) = F (x , q1, q2 , q3, q4 ) is a parametric model of the rate of population increase depending on the level of the active resource per beaver, where x = Rk(a) Pk ; Q is a vector

Pk(1) , Pk(2) , Rk(a) , Rk(p) ,

q2(i)

q1(i)

q3(i) q3(i) +

(q3(i) )2 + (q4(i) )2

q2(i) +

(q3(i) ) 2 + (q4(i) ) 2 q3(i) +

(q3(i) )2 + (q4(i) )2

; i = 1, 2

(2)

5

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

Analysis of competitive interactions

between the species, which is practically improbable (see Appendix). Case 3. Pc(1) > 0, Pc(2) = 0 . If the two species have identical death rates but their fecundities are different, the positive equilibrium population size cannot be realized for the species (Cf) with lesser fecundity.

We analyze the population dynamics of competing Cc and Cf species via stationary solutions of system (1). Stationary solutions are obtained in the traditional way, i.e., by looking for the solutions that are constant over time, then investigating their asymptotic stability by the Lyapunov method in the first approximation. Additionally, the two-species model (1) is used to carry out computer experiments. The population (q1, q2 , q3, q4 ) and resource (q5 , q6 , q7 ) parameters for Cf are taken from the single-species model based on the long-term monitoring data on Cf population dynamics in 6 Nature Reserves (Petrosyan et al. 2016). The parameters for Cc are determined from conditions (2) and (A.11) (see Appendix). Parameters q1 and q2 for Cc are calculated from the fecundity ratio (Frat) of Cc to Cf averaged to 1.67 ( ± 2SD = 0.346) (Lavrov 1960, Lathi, Helminen 1974, Kanshiev 1983, Rosell, Parker 1995, Danilov 2009, Danilov et al. 2011, Parker et al. 2012). The other parameters of the model, q3 , q4 , q5 , q6, q7 do not differ between Cf and Cc in our computer experiments. In all the computer experiments, the initial abundance of Cf is equal to the equilibrium level, Cf*, obtained in a single-species model (Petrosyan et al. 2016) or to the maximum value of the population size (Cf**) during the monitoring period if Cf** < Cf*. Such an approach is necessary to simulate conditions for the introduction of alien Cc into the habitats of Cf. Cc beavers are introduced by groups varying from 2 to 24 beavers at a step of two individuals (2, 4, 6, … 24). For each scenario of introduction, we determine the duration from the moment of Cc introduction until the phase of competitive Cf exclusion by Cc. For the odd numbers of 3, 5, 7, … 23 of introduced Cс beavers, the duration of the coexistence phase for all reserves are established based on nonlinear regression models (RM) using the estimates for even Cc numbers. We evaluate the suitability of these models using residual analysis. It means that the fitted curves are considered acceptable if we cannot reject the following statistical hypotheses: 1) that the average values equal zero, and 2) that the residuals obey a normal distribution according to the Kolmogorov – Smirnov, Shapiro – Wilk, and Watson tests. These curves adequately describe the relationships between the number of introduced Cc beavers and the duration of the two-species coexistence in Nature Reserves.

Case 4. Pc(1) = 0, Pc(2) > 0 . The zero equilibrium population size of the first species (Cf) is asymptotically stable for all studied Reserves (see Appendix), i.e., if the two beavers species have similar death rates, then the positive stationary population size of the species (Cc) with higher growth rate (fecundity) can be realized. Sufficient conditions for the asymptotic stability of the unique stationary solution in the above form are presented in the Appendix. The analysis of stationary solutions has shown that, in general, for any values of demographic and resource parameters that meet condition (2), the competitive exclusion does occur, which illustrates the Gause (1934) exclusion principle. In other words, if two species use the same food resources and occupy the same habitat, the Cc population grows faster, thus, constraining the Cf population growth. Therefore, if environmental conditions are stable, the Cf population may get extinct. However, the duration of coexistence and the rate of Cf exclusion depend on the initial abundances of both species, their demographic parameters and resource parameters. The rates of Cf exclusion have been established under various environmental conditions in the computer experiments. The initial numbers of introduced Cc individuals vary from 2 to 24 for each Reserve. Fig. 4 shows the duration of species coexistence phase until the onset of the Cf exclusion phase in 6 Reserves at different initial abundances of beaver populations. Model parameters for different Reserves are presented in Table 3; they have been used in all the subsequent computer experiments. The analysis has shown that there are no causes for rejecting the hypotheses about the normal distribution of residuals (with zero arithmetic average) for all models of regression (Table 4), i.e., these curves accurately illustrate the coexistence duration depending on the number of introduced Cc individuals in different Reserves. The dynamics of beaver population sizes for the case when Cc = 12 is given in Fig. 5. In LNR, as shown in Fig. 4a, the duration of Cf exclusion depends on the initial abundance and varies within a wide range, from 5 to 63 years. In particular, at the initial abundance of Cf equal to 25 ind. and that of Cc equal to 12 ind., the populations of both species increase at first, following the irruptive pattern of development (Fig. 5а). In this case, Cc begins to suppress Cf after 31 years. Such a pattern of both species development of has been observed at all tested initial population sizes in the computer experiments, i.e., Cf = 25, Cc = 2, 4, … 24 ind., and for the set of q parameters (see Table 3). Therefore, the pattern of the alien Cс population development in LNR does not significantly differ from that of Cf in the single-species case (Petrosyan et al. 2016). When the initial Cc population consists of 24 ind. (carrying capacity level), Cf is excluded faster, for 5 years. The habitat capacity measured by the q parameters can support 26 Cc ind. on the average, with insignificant fluctuations from 23 to 31 individuals. Hereafter, the carrying capacity equals the equilibrium number of beavers, which has been determined from a single-species model (Petrosyan et al. 2016). In DNR, the duration of Cf exclusion phase depends on the initial abundance of both species as a nonlinear function (R2 = 99.2 %, Fig. 4b). When the initial number of Cf beavers equals 250, of Cc beavers equals 2, 4, …, 24, the duration of coexistence phase until the beginning of exclusion phase equals 133, 116, …, 66 years, respectively. If the initial number of Cf = 250 ind. and that of Cc = 12 ind., the coexistence phase duration until the beginning of exclusion phase equals 83 years (Fig. 5b). Petrosyan et al. (2016) indicated that the Cf population growth rate in the DNR exhibits a multistep pattern with quasi-periodic oscillations. Fig. 5b illustrates the growth pattern of the Cc population that pertains for all the tested initial Cc population sizes. The habitat capacity measured by the q parameters for Cf can support

Results The analysis has revealed that the system of Eq. (1) has a stationary (time-independent) solution: Pc(1) , Pc(2) , Rc(a) , Rc(p) (d ) (a ) (p) Rc ). This solution must be consistent with the fol(R d = 1 R c lowing system of algebraic equations:

q2(1) + (q1(1)

q2(1) )

q2(2) + (q1(2)

q2(2) )

q6(1) Pc(1)

Rc(a)

q3(1) Pc

(q4(1) Pc )2 + (Rc(a) Rc(a)

q3(1) Pc ) 2

q3(2) Pc

(q4(2) Pc )2 + (Rc(a)

q3(2) Pc ) 2

Pc(1) = 0

Pc(2) = 0

q6(2) Pc(2) + q5 Rc(p) = 0

q5 Rc(p) + q7 Rc(d) = 0 Rc(d) = 1

Rc(a)

Rc(p) , Pc = Pc(1) + Pc(2).

(3)

The existence and stability analysis of stationary solution is presented in Appendix. The characteristics of all possible stationary solutions of system (1) are given below. Case 1. Pc(1) = 0, Pc(2) = 0 . Such a stationary solution is unstable. It means that the populations of two beaver species of cannot go extinct simultaneously. Case 2. Pc(1) > 0, Pc(2) > 0 . In order that stationary population sizes exist for both species, it is necessary that all the q parameters be equal 6

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

Fig. 4. Model estimates of the Cc-Cf coexistence phase duration in 6 Reserves as a function of the initial number of introduced Cc beavers: X axis shows the number of introduced Cc individuals; Y axis, years after the introduction; dots represent the model estimates; the bold curve is a regression model of this dependence (with a 95% predicted confidence interval). Panel a represents LNR; a, DNR; b, CFNR; c, PTNR; d, ONR; e, KNR; f. The initial number of Cf (ind.) in the Reserve: a, 25; b, 250; с, 350; d, 45; e, 300; f, 350 (parameter values as in Table 3).

beaver population sizes is also nonlinear (R2 = 98.9 %) (Fig. 4e). If the initial population sizes are 300 ind. for Cf and 2, 4, …, 24 ind. for Cc, the exclusion phase will begin after 132, 115, …, 68 years, respectively. In contrast to the other Reserves, the Cc population growth in ONR is properly approximated by a logistic curve (Fig. 5e), the pattern at the other initial population sizes being similar. The ONR habitat capacity measured by the q parameters for Cf population can support 474 ± 127 beavers. The KNR data also give evidences of a nonlinear dependence of the coexistence phase duration on the initial beaver abundances (R2 = 99.3 %, Fig. 4f). If the initial Cf population size equals 350 ind., the initial Cc population size equals 2, 4, …, 24, the exclusion phase will begin after 108, 92, …, 60 years in accordance with Cc initial numbers, respectively. Fig. 5f demonstrates that the phase of two-species coexistence lasts about 70 years. Besides, the dynamics of Cc population growth rate follows a multistep pattern at the other initial population sizes, too. The KNR carrying capacity measured by the q parameters can reach 781 ± 124 beavers. We have analyzed the stationary solutions (Cases 3, 4; see Appendix) assuming that the death rates of both species are equal. Unfortunately, there is little published data for comparison. The annual death rate in a sustainable Cc population is 26–30%. The death rate of adults above 2.5 years is 30–32%. During the first year of life, the death

on average 713 ± 22 Cc ind. at the end of the computer experiment. In CFNR, how the duration of coexistence phase depends on the initial population sizes is also described by a nonlinear function (R2 = 99.9%, Fig. 4c). Yet, in this Reserve, Cf is displaced markedly slower (Fig. 4c). If, in particular, the initial Cf population size is 350 ind. (Fig. 5c), that of Cc equals 2, 4, …, 24 ind., then the Eurasian beaver extinction will begin after 152, 133, …, 80 years, respectively. Fig. 5c shows that the Cc population dynamics starting with 350 ind. of Cf and 12 ind. of Cc is following a multistep pattern. Petrosyan et al. (2016) indicated that Cf population dynamics was similar in a single-species case for this Reserve. The capacity of this reserve habitat measured by the q parameters can support 316 ± 120 beavers. In PTNR, if the Cf initial population size equals 45 ind. (Fig. 5d) and that of Cс equals 2, 4,…, 24, then the duration of coexistence phase depends on the Cc initial population size as a nonlinear function (R2 = 99.7 %, Fig. 4d). Cf would be excluded after 389, 300, …, 70 years in accordance with the introduced numbers of Cc beavers. If 12 Cc individuals are introduced, the Cf exclusion will start after 146 years (Fig. 5d). Petrosyan et al. (2013, 2016) showed that Cf population dynamics in PTNR followed a one-step pattern with quasi-periodic oscillations and the capacity of its habitat measured by the q parameters (Table 3) can support 42 ± 3 Cc ind. in this Reserve. In ONR, the dependence of the coexistence duration on the initial Table 3 Parameter values for the six Nature Reserves. Nature Reserve LNR DNR CFNR PTNR ONR KNR

Species Cf Cc Cf Cc Cf Cc Cf Cc Cf Cc Cf Cc

q1 0.112287 0.142287 0.279701 0.304701 0.584466 0.607466 0.29157 0.29837 0.192945 0.215945 0.144297 0.215297

q2 0.010014 0.040014 0.001967 0.026967 0.209638 0.232638 0.12684 0.133640 0.000235 0.023235 0.035 0.106

q4

q3

-06

q5

q6

q7

0.000824

9.88×10

0.004985

0.000166

0.3427

1.08×10-05

0.000135

0.013125

8.3×10-06

0.010938

0.000227

2.5×10-04

0.007512

8.68×10-06

0.004352

0.000438

-05

8.49×10

0.023921

0.000338

0.036723

2.1×10-05

2.41×10-05

0.017014

1.86×10-05

0.01835

2.78×10-05

1.85×10-06

0.014714

1.98×10-05

0.011488

7

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

Table 4 Results of residual analysis (deviations of the species coexistence duration from the model values). Nature Reserve LNR DNR CFNR PTNR ONR KNR

Mean ± SE -6

0.5 × 10 ± 1.28 0.2 × 10-6 ± 1.71 –0.2 × 10-4 ± 0.67 –0.5 × 10-4 ± 5.04 0.1 × 10-4 ± 1.97 –0.4 × 10-4 ± 1.16

Kolmogorov–Smirnov

Shapiro–Wilk

Watson

0.19 (P = 0.27) 0.17 (P = 0.49) 0.21 (P = 0.15) 0.2 (P = 0.19) 0.11 (P > 0.87) 0.19 -5-5(P = 0.27)

0.93 0.94 0.91 0.88 0.96 0.95

0.058 0.062 0.061 0.082 0.029 0.063

rate is the least owing to the protection of young individuals in the parental colony. During the second year of life, it is about 40% since the yearlings begin to disperse, leaving the parental colony. During the following 2 years, when the beavers spread along new habitats and their individual mass continuously increases, the death rate also amounts to about 40%. When the individual growth ceases at the age of 4.5–10.5 years, the death rate decreases down to about 20%. The annual death rate is similar for males and females. The main causes of natural death for Cc are wolf predation, tularemia, shortage of food, and intraspecific fights (Payne 1989). There are few data on the Cf death rates. In Voronezh Reserve, during the period of 1937–1963, the age composition of dead beavers (n = 192) distributed as follows: 57% of adults, 26% of yearlings and two-year-old beavers, and 17% of under-yearlings. In KNR for the same period (n = 55), these estimates constituted 38%, 42% and 20%,

(P = 0.38) (P = 0.49) (P = 0.3) (P = 0.09) (P = 0.69) (P = 0.53)

(P = 0.34) (P = 0.29) (P = 0.31) (P = 0.15) (P > 0.84) (P = 0.29)

respectively (D’yakov 1975). During the same period in Voronezh Reserve, approximately 2% of the average annual number of beavers died (Barabash–Nikiforov et al., 1961). The main causes of beaver death in Voronezh Reserve were beaver bites during their fights, predator attacks and natural disasters. In KNR, they were natural disasters, poaching, predators, and beaver bites (D’yakov 1975). In ONR, 47 deaths of beavers were recorded in 1967–1972. The main natural causes of those deaths were beaver bites, natural disasters, asphyxiation under ice, hunger, and predator attacks (Kudryashov 1975). The main difference in death rates between the two species is likely their different resistances to tularemia. While the high Cf resistance to tularemia was shown both in field studies (Borodina 1966) and in experimental research (Barabash–Nikiforov et al. 1961), Cc is, on the contrary, known to be vulnerable to this disease. Cases of its mass death were recorded repeatedly (Novak et al. 1987).

Fig. 5. Simulated population dynamics after the introduction of 12 Cc ind. in 6 Reserves: X axis, years after the introduction; Y axis, the absolute numbers of individuals; black line for Cf population; red line for Cc; blue line is a logistic curve regressed for the simulation data. Initial numbers of Cf are given on the corresponding panels (parameter values as in Table 3). The same panels as in Fig. 4gr5 8

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

Fig. 6. Simulated population dynamics of Cf (black line) and Cc (red line) beavers after the introduction of 12 Cc ind. in 6 Reserves populated by Cf for the case where the death rate of Cc is higher than that of Cf, with the same q parameters (Table 3). The same panels and axes as in Fig. 5.

If the death rates of the species are different, the higher fecundity of Cc can be counter-balanced by its greater death rate providing a longterm coexistence of two species. Model prediction of this case is shown in Fig. 6. Figs. 6a–f indicate that the introduction of 12 Cc individuals in LNR, DNR, PTNR, CFNR, ONR, and KNR, provides eventually the equilibrium population size of Cf equal to 10, 36, 17, 14, 25 and 16 ind., respectively. Our computer experiments show that Cc can long survive in all the Reserves if the death rate is sufficiently high.

degradation, and regeneration resources. 2. Our long-term monitoring data (50–70 years) on the beaver populations in six Natural Reserves and the parameters of the single- and two-species models allow us to describe the main patterns of Cf and Cc population dynamics under competitive interactions between two species. 3. While our single-species model can assess the demographic (q1, q2 , q3, q4 ) and environmental (q5 , q6 , q7 ) parameters on the basis of the long-term monitoring data on the Cf populations, the two-species model can assess these parameters for Cf and Cc. The parameters q1, q2 and q1*, q2* are measured by means of available data on the fecundity and family sizes of both species. For example, the average family size of Cf is 3.8 (SE = 1.0, N = 13 studies; range: 2.4–5.5 ind.), and that of Cс is 5.2 (SE = 1.4, N = 51 studies; range: 2.7–9.2 ind. (Rosell, Parker 1995, Parker et al. 2012). Parker et al (2012) concluded in their review that the mean foetus number was less for Cf ( 2.5) than for Cc ( 4.0), and it was supported by Pyotr I. Danilov’s data (Danilov 1995, Ruusila et al. 2000), who found that the litter size of Cf = 1.9 and that of Cc = 3.2. The average size of Cc kit was reported 4.7 (n = 9) in Finland (Lathi, Helminen 1974), and 3.3 (n = 87) in Karelia (Danilov 2009). The analysis of a large set of data on the fecundity of beavers (Cf: n = 896, Cc: n = 1080) evidenced that the Cc fecundity, on average = 3.27 (min–max: 1–12, embryos per female) is higher than that in Cf, on average = 2.92 (min–max:1–9, Kanshiev 1983; Danilov 2009). In addition, Danilov et al. (2011) showed that fecundity is one of the greatest distinctions in the reproductive biology of the Cc. It is much higher in Cc than in Cf (the fecundity of Cc = 4.0, maximum

Discussion The long-term data on Cf population dynamics in Lapland, Darwin, Central-Forest, Prioksko-Terrasny, Oka, and Khoper Reserves located in the northern, central and southern parts of this species range in European Russia were used earlier in a single-species population model for the Eurasian beaver (Zavyalov et al. 2015, European beaver 2012, Petrosyan et al. 2016). The single-species model revealed four patterns of the beaver population dynamics, namely, irruptive (LNR), single-step with quasi-periodic oscillations (PTNR), multistep with quasi-periodic oscillations (DNR, CFNR, KNR), and periodic oscillations around a logistic growth curve single-species model could predict the main patterns of Cf populations dynamics in different habitats, we have modified it to analyze the competitive interactions between two beaver species during their coexistence. The advantages of the two-species model are the following: 1. This model involves interactions of both species with environmental conditions including the intensity of exploitation (usage), 9

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

foetuses = 12, the fecundity of Cf = 2.2, maximum foetuses = 6). The higher fecundity of Cс (2.86 kits per female) than that of Cf (1.93 kits per female) was also documented for beavers in captivity of Voronezh Nature Reserve (Lavrov 1960). 4. In the model, two species have identical environmental parameters (q5 , q6 , q7 ) in accordance with published data. Reserchers formerly argued that the Cc builder activity is more effective and that Cc is better adapted to critical environmental conditions, and more often feeds on grey alder (Alnus incana) (Ruusila 1997, Danilov, Kan’shiev 1983). However, the current data demonstrate no differences in their building activity and diet. Both species likely occupy the same ecological niche (Danilov et al. 2011, Danilov, Fedorov 2015). We divide population dynamics into two phases after Cc introduction into habitats occupied by Cf. During the first phase, the Cc population size increases while the number Cf beavers decreases. The abundance of Cc is less than that of Cf throughout the first phase. Hence, this stage can be called the species coexistence phase. After this phase, the Cc population size becomes equal to that of Cf, afterwards Cc begins to suppress Cf. While Cf is dying out, the population size of Cc is tending to a level in accordance with habitat capacity measured by the q parameters. By the end of the second stage, Cc excludes Cf completely. Cf suffers extinction, while the abundance of the other species establishes at the maximal level determined by the habitat carrying capacity that is described by the q parameters. We call the second period as competitive exclusion phase. In simulation experiments, the number of introduced Cc individuals has varied from 2 to 24. Such a size of the introduced beaver group corresponds to the number of beavers reintroduced in Eurasia last century. In particular, 14 (Kataev 2015), 8 (Korablev et al. 2011), 4 (Zavyalov et al. 2010), 23 (Kudryashov 1975), and 22 (D’yakov 1975) Cf ind. were released into LNR, CFNR, PTNR, ONR, and KNR, respectively. Seven Cc beavers and 19 Cf ones (Lathi, Helminen 1974) were introduced in Finland, in 1935–1937. In 1998 and 1999, 40 and 12 Cf ind. were released in Belgium, respectively (Halley, Rosell 2002). In France, Cf were introduced in groups including 3 to 21 beavers, totally 273 Cf ind., between 1957 and 1988 (Dewas et al. 2012). The population dynamics predicted by our model in LNR, which is characterized by extreme environmental conditions, demonstrates that after the introduction of 2 – 24 Cc individuals, the phase of coexistence lasts from 63 to 5 years. The maximum number of Cc individuals during the phase of competitive exclusion does not exceed 30 ind. (Fig. 5a). These estimates confirm the results of repeated Cc introductions into similar ecological conditions in the north-west of Finland, where the current population of Cc did not exceed 50 ind. (Lathi, Helminen 1974, Parker et al. 2012). The outcome of our computational experiments shows that, for different initial numbers of Cc, the predicted population growth is always irruptive. Similar patterns of population dynamics were observed for Cf in two regions of Sweden (Hartman 2003) and in LNR (Petrosyan et al. 2016). Noteworthy is that, when 24 ind. of Cc were introduced, the phase of coexistence lasted 5 years, i.e., this phase involved 2-3 generations. Published data indicate that Cс more frequently develop to maturity at the end of their second year of life as compared to Cf (Henry, Bookhout 1969, Danilov 2009). Since environmental conditions in PTNR are unfavorable, the duration of Cf exclusion phase should last longer than in the other study Reserves (Fig. 4d). In this Reserve, the model predicts that Cf population dynamics for single-species community tends to equilibrium state with quasi-periodic oscillations with periods of 14–26 years (Petrosyan et al. 2013). Fig. 5d shows that after the introduction of 12 Cc individuals, the Cf population size decreases according to quasi-periodical oscillations. In this case, Cc can exclude Cf for 400 years. The population of Cc, on the contrary, grows with a linear trend and small oscillations, yet the quasi-periodicity in Cc population dynamics is not well expressed. The relations between two species during their coexistence phase depend on the q parameters that describe the habitat carrying capacity.

In the suboptimal habitats of DNR, two species coexisted for 136–66 years after the introduction of 2 to 24 of Cс individuals (Fig. 4b). The model shows that Cf gets extinct when Cc reaches the first step of its population growth. A multistep pattern of population dynamics was first found for Cf in this Reserve (Petrosyan et al. 2016); yet this pattern is better expressed in Cc population dynamics. In CFNR and KNR suboptimal habitats, Cf is displaced by Cc at all the tested initial numbers of introduced Cc individuals (Fig. 4с, f) much faster than in PTNR (Fig. 4d) due to the high rate of food resource renewal. In particular, when 12 ind. of Cc are introduced, the two beaver species coexist for 100 years in CFNR, and 71 years in KNR. In CFNR and KNR, the model also predicts the multistep population dynamics similar to that in DNR (Fig. 5c, f). In KNR, the multistep growth of Cc population is very distict: the population reaches the first step for 25 years after the introduction while in CFNR it takes 150 years. Noteworthy is that the multistep pattern of Cc population dynamics is more expressed than that of Cf. In favorable conditions of ONR, Cf is displaced fairly quickly too at all initial numbers of Cc individuals varying from 2 to 24 (Fig. 4e). In this Reserve, Cc population dynamics, in contrast to other studied Reserves, demonstrates logistic oscillatory growth (see Fig. 5e, logistic curve), and the pattern is similar for other initial population sizes. The phase of coexistence lasts about 85 years in ONR at 12 initial Cc individuals. The model shows (Fig. 5) that if 12 beavers of Cc are introduced, the phase of Cf and Cc coexistence would take 31, 71, 83, 85, 100, and 146 years in LNR, KNR, DNR, ONR, CFNR, and PTNR, respectively. We relate the high rate of competitive exclusion in LNR to its extremely low carrying capacity for beavers because of the long regeneration time of trees and shrubs in the Far North. The moderate rates of competitive exclusion in the last four Reserves are provided by sufficient availability of food resources. However, in PTNR, the food resources for beavers are rather scarce (Zavyalov et al. 2010) yet the predicted rate of the competitive exclusion is quite low, about 146 years. There may be several reasons for this phenomenon. First, beaver food resources restore faster in PTNR compared to LNR and previously abandoned habitats are repopulated after four years (Zavyalov et al. 2016). In contrast, the restoration of woody vegetation in LNR lasts at least 50 years (Kataev 2011). Second, beavers in PTNR frequently change their habitats, sometimes several times a year (Zavyalov et al. 2016). Therefore, there will be hard for Cc to settle if all suitable habitats are already occupied. Third, the beavers demonstrate high building activity in PTNR, they continuously create new dams and repair old ones; numerous dams remain long after beavers abandon the habitat (Goryainova et al. 2014, Zavyalov et al. 2016). It is known that all these factors can ultimately enlarge carrying capacity for beavers (Logofet et al. 2015). Thus, even under the current scarcity of food resources owing to long-term exploitation by beavers, the carrying capacity in PTNR is higher than in LNR and the environmental conditions for beaver development are much better. The model results indicate that Cc would penetrate the Cf range diffusively. The rate of Cc penetration into the central part of the Cf range is still unknown. The rate of dispersal among vacant habitats is similar in both species (Danilov et al. 2007). The following patterns of beaver migrations are distinguished: short-term migrations of yearlings; dispersals of animals at the age of 2 creating new colonies; migrations of the whole colonies to the other water bodies within their ranges; individual migrations of adults that lost their partners. The distance of beaver migrations varies greatly, and it frequently depends on the availability of favorable and free habitats. The direction of migration can be either upstream or downstream. Beavers are also able to overcome watersheds overland. The distances of migrations in uninhabited areas are commonly greater than within habitats with established populations (Saveliev et al. 2010). Since the Cc dispersals were analyzed in inhabited areas, their pattern can differ from that among free habitats. Thefore, Cc would likely migrate from the north to the south of 10

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

frontiers (Khlyap et al. 2011, Bonesi, Tom 2012). The American mink has replaced the European mink (Mustela lutreola L. 1761) almost everywhere (Maran, Henttonen 1995, Tumanov 2009). Its further population growth in the north-west of Russia might negatively affect the state of small populations of the native European mink (Bobrov et al. 2008). There are some more alien species that have caused similar negative effects in Europe. In particular, the Eastern grey squirrel (Sciurus carolinensis Gmelin, 1788) has gradually excluded the native red squirrel (S. vulgaris L., 1758) in the UK and Italy (Wauters et al. 2002, Martinoli et al. 2010).

European Russia. In general, we suggest that vacant habitats are the most vulnerable to Cc colonization, while favourable areas are more resistent. Cf populations might remain as a number of isolated populations within the current range of this species. The model predicts that Cc would displace Cf in all the study Reserves located in different ecological conditions in the northern, central, and southern parts of the Cf current range. However, if the Cc death rate is high enough to counter-balance the high Cc birth rate, then a long-term coexistence of both species is likely. Nevertheless, the equality of death rate and birth rate may not be sufficient for their coexistence in some cases since Cc mature earlier than Cf, which gives another advantage of Cc over Cf. Although Cf sexually matures at the end of the second year of life, only a small share of young animals participate in reproduction. In Voronezh Reserve, this share did not exceed 8–10% (Lavrov 1956, 1981); it is 12% in KNR (D’yakov 1975); 7% in ONR (Kudryashov 1975); 6.6% in the Leningrad Oblast and the south of Karelia (Kanshiev 1983). On the contrary, Osborn (1953) indicated that 21.2% of Cc females start reproduction at the age of 1.5–2 years. To conclude, complex models are very useful to predict the consequences of alien species invasions under different environmental conditions and at various initial numbers of individuals. The suggested two-species model can predict population dynamics of the two beaver species if we know environmental conditions, the initial population sizes and the duration of their coexistence. For example, 2–8 ind. of Cc from Finland were introduced into the vacant habitats in Bavaria in 1972 (see Introduction). Then, by 1997, Cc population could have reached 266 ( ± 34) ind. if the ecological conditions for beavers in Bavaria had been suboptimal, i.e., similar to those in DNR, CFNR, and KNR. If the conditions had been optimal, the population size could be 577 ( ± 87) ind. If the minimum Cc population growth rate in Bavaria had been equal to that in Finland (known from the real monitoring data, Lathi, Helminen 1974) for 25 years, then Cc population size could reach 550 ( ± 50) individuals. These estimates enable us to predict a likely secondary source of Cc invasion that could have formed in Bavaria by 1997. It was likely that the beavers from such a source were removed and introduced into other European countries in 1998–2002. The 66 Cc individuals that were detected, removed, sterilized or transferred to zoos in Germany (RLP, NRW), France, Belgium, and Luxembourg (Dewas et al. 2012, Frosch et al. 2014) constituted about 12 % of the population size predicted by the model. Our model predicts that if strategic plans for the eradication of the North American beaver including recommendations of Parker et al. (2012) are not accepted, the Eurasian beaver will likely get extinct in a significant part of its restored range in the 21 st century. The history of biological invasions has well demonstrated that, in spite of certain success in reintroduction and restoration of some species, there can be serious negative consequences for native species, especially those that are ecologically close to invaders. For example, the American mink (Neovison vison (Schreber 1777)) intentionally introduced into various regions of Western Europe in 1920–1930 and Russia in 1928 was successfully naturalized. Currently, its range extends across Eurasia from the Atlantic to the Pacific and nearly reaches the maximum possible

Conclusions Ecological and biological similarities between the two beaver species can result in the competitive exclusion of one of them. Since Cc was shown to have higher fecundity, larger family size, and early sexual maturation, it is likely to be superior over Cf and to exclude Cf when they compete despite the similarities in the life span and population age structure. One of the top priorities for the conservation of the Eurasian beaver at the continental, national and regional levels is to develop an international program for monitoring the status of Cc populations in Eurasia. Such a program should comprise the strategies developed by the international scientific communities, such as Aichi Biodiversity Target 9 (Convention on Biological Diversity 2001) and EU Biodiversity Strategy 2020 Target 5 (European Commission 2011), including the measures to constrain population sizes of the North American beaver, to reduce its ranges, and, at the slightest threat, to completely eradicate it. The necessity of such a program is dictated by the obligations arisen from the international treaties on biodiversity conservation, and, especially, from the Aichi Biodiversity Target 9: “By 2020, invasive alien species and pathways are identified and prioritized, priority species are controlled or eradicated, and measures are in place to manage pathways to prevent their introduction and establishment" (Convention on Biological Diversity 2001). Acknowledgements We are thankful to all the staff of Nature Reserves for their assistance in creating the databases of long-term (1934–2015) monitoring data on Eurasian beaver (Castor fiber). We would like to thank I. Y. Feniova (A.N. Severtsov Institute of Ecology and Evolution, Russian Academy of Sciences) for her help with English translation. The original manuscript was greatly improved by comments from anonymous reviewers. The authors are also grateful to ESRI (USA) for providing a free-of-charge licensed version of ArcGIS Desktop Pro 10.4.1 (Esri Sales Order number 3128913; Esri Delivery number 81833751, User customer number 535452). Development of mathematical models of species population dynamics, statistical analysis, data interpretation, literature review and material preparation for publication are supported by the Russian Science Foundation (grant no. 16-14-10323). Data collection on the territory of the Reserves was supported by the programs of the Presidium of the Russian Academy of Sciences № 41 "Biodiversity of natural systems and biological resources of Russia".

Appendix A. Asymptotic stability analysis of stationary solutions Asymptotic stability of the stationary solutions is studied by the Lyapunov method in the first approximation. Stability theory of finite-difference equations (e.g., Romanko 2012) states that for stability of a stationary solution, it is necessary that the eigenvalues of the Jacobi matrix derived for the stationary solution do not exceed 1 in absolute value. Moreover, if all the eigenvalues are strictly less than one, the stationary solution is asymptotically stable. If at least one eigenvalue is strictly greater than one, the stationary solution is unstable. Case 1. Pc(1) = 0, Pc(2) = 0 . Since q5 0 and q7 0 , using Eq. (3), we can derive thatRc(a) = 1, Rc(p) = Rc(d) = 0 . With respect to the system of Eq. (1), the Jacobi matrix derived for the stationary solution can be written as follows:

11

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( )( ) P1

(A.1)A =

P1

P1

P1

P (1)

P (2)

R(a)

R(p)

P2 P (1)

P2 P (2)

P2 R(a)

P2 R(p)

Ra P (1)

Ra P (2)

Ra R(a)

Ra R(p)

( )( )( )( ) Rp

Rp

Rp

Rp

P (1)

P (2)

R(a)

R(p)

1 + q1(1)

0

0

0

0

1 + q1(2)

0

0

=

q6(1)

q6(2)

0

0

1 q7 1

q5 q5

q7

c

Where, P1, P 2 , Ra , Rp are the right-hand sides of the first, second, third, and fourth equations of system (1), respectively, c means that the partial derivatives of the right-hand sides w.r.t its variables are calculated at the stationary solution. The characteristic equation to determine the eigenvalues of matrix (A.1) takes on the following form:

E) = (1 + q1(1)

det(A

)(1 + q1(2))

)[(1

)(1

q5

q7

(A.2)

) + q5 q7] = 0

The eigenvalues are found as follows: 1

= 1 + q1(1) ,

2

= 1 + q1(2) ,

3

=1

q5,

4

=1

(A.3)

q7

and q1(2)

>0 > 0 (see model) that 1 > 1 and Hence, it follows from the populations of two species cannot get extinct simultaneously. q1(1)

2

> 1 and, consequently, the stationary solution is unstable. It means that

Case 2. Pc(1) > 0, Pc(2) > 0 . In this case, system (3) reduces to

(A.4)

q2(1) + (q1(1)

q2(1) )

q2(2) + (q1(2)

q2(2) )

q6(1) Pc(1)

q3(1)

xc

(q4(1) )2 + (x c

q3(1) )2

q3(2)

xc

(q4(2) )2 + (x c

q3(2) )2

=0 =0

q6(2) Pc(2) + q5 Rc(p) = 0

q5 Rc(p) + q7 Rc(d) = 0 Rc(d)

Rc(a)

=1

Rc(p) , Pc = Pc(1) + Pc(2), xc = Rc(a) Pc . (a)

It follows from (A.4) that the first two equations depend only on the value of x c = Rc (P1c + P2c ) to be found. It means that x c must simultaneously obey the first two equations differing from each other only by the set of q parameters for different beaver species. After solving the first and second equations of (A. 4) we obtain

q2(1) q4(1)

x c = q3(1)

q1(1) (q1(1)

2q2(1) )

q2(2) q4(2)

= q3(2)

q1(2) (q1(2)

2q2(2) )

(A.5)

According to (A. 5), the existence of stationary solution beaver species 1 and 2:

q2(1) q4(1)

q3(1)

q1(1) (q1(1)

2q2(1) )

Pc(1)

> 0,

Pc(2)

> 0 implies the following condition (equality) for the q parameters of

q2(2) q4(2)

= q3(2)

q1(2) (q1(2)

2q2(2) )

(A.6)

Since the restriction on the parameters of the model in the form of equality (A. 6) for arbitrary parameters of the two beaver species is unlikely in nature, the stationary solution Pc(1) > 0, Pc(2) > 0 can hardly be realized. Case 3. Pc(1) > 0, Pc(2) = 0 . It is necessary and sufficient for the stationary solution to exist if the following condition is satisfied (Petrosyan et al. 2013): (A.7) q3(1) +

q6(1) q5

+

q6(1)

q1(1) (q1(1)

q7

2q2(1) )

q2(1) q4(1) > 0.

Moreover, Petrosyan et al. (2013) showed that there is only one stationary solution Pc(1) > 0 . If condition (A.7) is satisfied, the stationary values of the beaver population sizes and resources are found as follows:

Pc(1) =

q1(1) (q1(1) (q3(1)

Rc(a) = 1

+

q6(1)

q5 +

q6(1)

q7 )

q1(1)

2q2(1) ) (q1(1)

2q2(1) )

q2(1) q4(1)

,

(q6(1) q5 + q6(1) q7 ) Pc(1), Rc(p) = (q6(1) q5) Pc(1), Rc(d) = (q6(1) q7 ) Pc(1).

(A.8)

In this case (see Case 1), the following matrix can be obtained by rearranging the two first rows and columns of the Jacobi matrix (since rearranging rows and columns does not change the eigenvalues):

12

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

1 + F(Rc(a) Pc(1), Q(2) ) f(Q(1) ) Rc(a) q6(2)

A=

0

Pc(1)

1

0 Q(1)

=

(q1(1) ,

0

f(Q(1) ) Rc(a) q6(1)

Pc(1)

0

1 q7

q5 q5

0

q2(1) ,

q3(1) ,

F(x , Q) = q2 + (q1

q4(1) ),

q2 )

q42

Q(2)

=

x

q3

(q1(2) ,

q2(2) ,

0

f(Q(1) )

q3(2) ,

1

q4(2)),

q7

Q = (q1, q2 , q3, q4 ) 3

q3) 2

+ (x

,

f(Q) =

[q1 (q1 2q2)] q4 (q1 q2)2

2

. (A.9)

Matrix (A. 9) has one real eigenvalue equal to 1

= 1 + F(Rc(a) Pc(1), Q(2) ) = 1 + q2(2) + (q1(2)

q3(2)

xc

q2(2) )

(q4(2) ) 2 + (xc

q3(2) )2

, xc = Rc(a) Pc(1)

(A. 10)

Here, as before and further, we assume that the death rates of both beaver species are identical and their fecundity rates do not change with time. This means that the q parameters of model (1) are related as

q3(2) = q3(1) , q4(2) = q4(1) , q1(2) = q1(1) + b, q2(2) = q2(1) + b , Since we have Pc = 1

Pc(1)

= 1 + b + q2(1) + (q1(1)

> 0,

Pc(2)

q2(1) )

(A. 11)

= 0 , it follows from (3) and (A.10) that q3(1)

xc

(q4(1) ) 2 + (xc

q3(1) )2

= 1 + b , xc = Rc(a) Pc(1)

(A. 12)

> > Since b > 0 here as the fecundity of Cc is greater than Cf matrix (A.9) has one real eigenvalue strictly greater than 1. Consequently, the stationary solution is asymptotically unstable according to Lyapunov’s theorem (Romanko 2012). In other words, if two species have similar death rates, the positive stationary population size of the species with lower growth rate (fecundity) cannot be realized. (q1(2)

q1(1),

q2(1) ),

q2(2)

Case 4. Pc(1) = 0 , Pc(2) > 0 . In this case, all the assumptions and restrictions are the same as in Case 3. Therefore, all the formulas and results obtained in Case 3 are applicable in this case with only small change, i.e., the upper scripts 1 and 2 must be permuted in these formulas and results. The Jacobi matrix derived for the stationary solution, the q parameters of model (1), and the eigenvalue are determined as follows (see A.13, A.14, A.15):

1 + F (Rc(a) Pc(2), Q(1) )

0

(Q(2) ) Rc(a) q6(1)

(Q(2) ) Rc(a) q6(2)

f

A= (A.13)

Pc(2)

1

f

0 Q(1)

=

(q1(1) ,

q2(1) ,

0 Pc(2)

f

0

1 q7

q5 q5

0 q3(1) ,

F (x , Q) = q2 + (q1

q4(1) ), q2)

Q(2) x

=

q3

q42 + (x

q3)2

(q1(2) ,

0

(Q(2) )

q2(2) ,

q3(2) ,

, f (Q) =

[q1 (q1

1

q4(2)),

q4 (q1

, q7

Q = (q1, q2 , q3, q4 ),

3

2q2)] 2 q2)2

;

(A.14)q3(1) = q3(2) , q4(1) = q4(2) , q1(1) = q1(2) + b, q2(1) = q2(2) + b , (A.15)

1

= 1 + b + q2(2) + (q1(2)

q2(2) )

x c q3(2) (2) 2 (q4 ) + (x c q3(2) )2

= 1 + b, x c = Rc(a) Pc(2)

Since b < 0 here as the growth rate (fecundity) of Cf is lower than Cc (q1(1) < q1(2), q2(1) < q2(2) ), the eigenvalue is strictly less than 1. It means that the stationary solution with Pc(1) = 0 may be asymptotically stable. The stationary solution Pc(2) > 0 , according to Petrosyan et al. (2013), does exist if and only if the following condition holds true:

q3(2) +

q6(2) q5

+

q6(2) q7

q1(2) (q1(2)

2q2(2) )

q2(2) q4(2) > 0

(A.16)

If condition (A.16) is satisfied, the stationary values of the beaver population sizes and resources are found as follows:

Pc(2) =

q1(2) (q1(2)

2q2(2) )

(q3(2) + q6(2) q5 + q6(2) q7(2) ) q1(2) (q1(2)

Rc(a) = 1

2q2(2) )

q2(2) q4(2)

, Pc(1) = 0

(q6(2) q5 + q6(2) q7 ) Pc(2), Rc(p) = (q6(2) q5) Pc(2), Rc(d) = (q6(2) q7 ) Pc(2)

(A.17)

In general, if condition (A.16) is not met, then the stationary solution does not exist. It means that the number of Cс beavers would change with time without tending to a nonzero stationary value. The same would happen if a stationary solution exists but is not stable. Since condition (A.16) is actually satisfied for all Reserves, we use computational experiments to illustrate the conditions under which all the eigenvalues of matrix (A.13) are less than 1 by module. Outcome of the analysis has shown that in, the case when the death rates of both species are not different, the fecundity of Cc is greater than that of Cf, and the condition (A.16) is satisfied, then the stationary population sizes of the both species can be represented by means of (A.17). For the asymptotical stability analysis, the eigenvalues, 1, 2 , 3 , 4 , of matrix (A.13) have been calculated. These values characterize the steady states of system (1) and they are presented in Table A1 for each of the Reserves. Table A1 shows that the eigenvalues 1 have real positive values for all the Reserves, and the eigenvalues that are expressed by complexly conjugated numbers are found for the LNR, DNR, CFNR, and PNTR. For the other Reserves, ONR and KNR, all the eigenvalues are less than 1, so the steady solution is stable and the periodic component is absent. The solution of the equations of population dynamics for the LNR, DNR, CFNR, and PNTR has a periodic component, 13

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

Table A1 Eigenvalues of matrix (A.13) to characterize the asymptotic stability of stationary solutions. Reserve name, Fecundity ratio, Frat

Matrix (A.13)

LNR, 1.17

0.965483 −0.175251 −0.000166 0 0.982072 −0.303165 −8.37E-06 0 0.978051 −0.764077 −8.7E-06 0 0.965623 −0.5838045 −0.000338 0 0.977592 −1.541944 −1.86E-05 0 0.962837 −1.517859 −1.98E-05 0

DNR, 1.57 CFNR, 1.16 PTNR, 1.15 ONR, 1.31 KNR, 1.22

1

Eigenvalues1 0 0.824749 −0.000166 0 0 0.696835 −0.000166 0 0 0.235923 −0.000166 0 0 0.416196 −0.000166 0 0 −0.5419442 −0.000166 0 0 −0.517859 −0.000166 0

0 213.02 1 −0.3427 0 368.5 1 −0.010938 0 928.73 1 −0.004352 0 709.61 1 −0.036723 0 1874.23 1 −0.018352 0 1844.95 1 −0.011489

0 0 0.004985 0.652315 0 0 0.001313 0.98775 0 0 0.007512 0.988136 0 0 0.0239212 0.939356 0 0 0.017014 0.964634 0 0 0.014714 0.973797

0.965483 0.910831 + 0.169119i 0.910831 – 0.169119i 0.655402 0.987678 0.982072 0.848454 + 0.195468i 0.848454 – 0.195468i 0.987967 0.978051 0.618046 + 0.090621i 0.618046 – 0.090621i 0.965623 0.933916 0.710818 + 0.179586i 0.710818 – 0.179586i 0.977592 0.962793 0.763112 0.303215 0.97285 0.962837 0.761338 –0.278249

In the decreasing order of absolute values.

and steady solutions are also stable: the maximum absolute values of the eigenvalues are less than 1 (Table A1). Outcomes show that for the stationary solution to be asymptotically stable (ASSS) in DNR, ONR and KNR, it is necessary that the values of Frat be equal or over 1.57, 1.31 and 1.22, respectively (Table A1). For the three remaining Reserves, LNR, CFNR, and PTNR, the ASSS is realized for the following values of Frat 1.17, 1.16, and 1.15, respectively. It means that for LNR, CFNR, and PTR, which feature the smallest forage resources, the minimal difference in Frat is enough for ASSS; in contrast, for DNR, ONR and KNR Reserves with maximum feed resources, the greatest value of Frat difference is needed for ASSS.

Danilov, P.I., Kan’shiev, V.Y., 1983. The state of populations and ecological characteristics of European (Castor fiber L.) and Canadian (Castor canadensis Kuhl) beavers in the northwestern USSR. Acta Zoologica Fennica. 174, 95–97. Danilov, P.I., Kanshiev, V.Y., Fedorov, F.V., 2007. Eurasian beavers of the European North of Russia. Science, Moscow (in Russian). Danilov, P., Kanshiev, V., Fyodorov, F., 2011. History of beavers in eastern Fennoscandia from the neolithic to the 21st century. In: Sjöberg, G., Ball, J.P. (Eds.), Restoring the European Beaver: 50 years of experience. Pensoft Publishers, Sofia, Bulgaria, pp. 27–38. Dewas, M., Herr, J., Schley, L., Angst, C., Manet, B., Landry, P., Catusse, M., 2012. Recovery and status of native and introduced beavers Castor fiber and Castor canadensis in France and neighbouring countries. Mammal Review. 42, 144–165. https:// doi.org/10.1111/j.1365–2907.2011.00196.x. D’yakov, Y.V., 1975. Beavers of the European part of the Soviet Union. Moscow Worker (in Russian). Ermala, A., Helminen, M., Lathi, S., 1989. Some aspects of the occurrence, abundance and future of the Finnish beaver population. Suomen Riista 35, 108–118. European beaver, 2012. In: Dgebudze, Yu.Yu., Zavyalov, N.A., Petrosyan, V.G. (Eds.), European beaver (Castor fiber L.) as a key species of a small river ecosystem (Prioksko-Terrasnyi Nature Biosphere Reserve). Association of scientific publications KMK, Moscow (in Russian). European Commission, 2011. Communication from the Commission to the Council, European Parliament, the European Economic and Social Committee and the Committee of the Regions. Our life insurance, our natural capital: an EU biodiversity strategy to 2020. SEC 2011https://doi.org/10.2779/24101. 540 final and SEC (2011) 541 final. Fedushin, A.V., 1935. The river beaver, its history, life and breeding experiments. Publishing house of Glavpushnina, Moscow (in Russian). Frosch, C., Kraus, R.H.S., Angst, C., Allgo¨wer, R., Michaux, J., Nowak, C., 2014. The Genetic Legacy of Multiple Beaver Reintroductions in Central Europe. PLoS ONE. 9 (5), e97619. https://doi.org/10.1371/journal.pone.0097619. Gurney, W.S., Lawton, J.H., 1996. The population dynamics of ecosystem engineers. OIKOS, vol. 76. pp. 273–283. Gause, G.F., 1934. The Struggle for Existence. Hafner, New York, USA. Goryainova, Z.I., Katsman, E.A., Zavyalov, N.A., Khlyap, L.A., Petrosyan, V.G., 2014. Evaluation of tree and shrub resources of the Eurasian beaver (Castor fiber L.) and changes in beaver foraging strategy after resources depletion. Russian journal of biological Invasions. 5 (4), 242–254. https://doi.org/10.1134/S207511171404002X. Grave, G.L., 1931. The river beaver within the USSR and its economic importance.

References Albov,, S.A., Khlyap, L.A., 2015. The current state of the mammalian fauna of the Prioksko-Terrasny Reserve. Proceedings of the Prioksko-Terrasny Reserve, 5 109–144 (in Russian). Anonymous, 1972. Finnish beavers shipped to Bavaria. Metsästys ja Kalastus 4, 52–53 (in Finnish). Barabash–Nikiforov, I.I., Dezhkin, V.V., D’yakov, Yu.V., 1961. Beavers of the Don basin, in: Ecology and practical issues (monographic essay. Proceedings of the Khoperskiy Nature Reserve, 5 5–115 (in Russian). Beavers in the reserves of the European part of Russia, 2018. N.A. Zavyalov L.A. Khlyap Velikiye Luki: Velikiye Luki printing houseProceedings of State Nature Reserve «Rdeysky». V. 42018. Zavyalov, N.A., Khlyap, L.A. (Eds.), Proceedings of State Nature Reserve «Rdeysky». V. 4 (in Russian with English abstracts). Beaver as a renewable resource, 2019. In: Sjöberg, G., Belova, Olgirda (Eds.), Beaver as a renewable resource. Swedish Forest Agency, Uppsala. Bobrov, V.V., Warshavsky, A.A., Khlyap, L.A., 2008. Alien species of mammals in the ecosystems of Russia. The Association of Scientific Publications KMK, Moscow (in Russian). Bonesi, L., Tom, M., 2012. Neovison vison Schreber (American mink). A Handbook of Global Freshwater Invasive Species. Earthscan, London, New York, pp. 378 390. Borodina, M.N., 1966. Materials for the study of the dynamics of the Moksha beaver population. Proceedings of the Mordovskiy Nature Reserve, 3 5–38 (in Russian). Brommer, J.E., Alakoski, R., Selonen, V., Kauhala, K., 2017. Population dynamics of two beaver species in Finland inferred from citizen-science census data. Ecosphere. 8 (9), e01947. https://doi.org/10.1002/ecs2.1947. Convention on Biological Diversity, 2001. Decision VI/23. Alien species that threaten ecosystems, habitats or species. Guiding principles for the prevention, introduction and mitigation of impact of alien species that threaten ecosystems, habitats or species. Visited 28.05.2018. http://www.cbd.int/sp/targets/. Danilov, P.I., 1995. Canadian and European beavers in Russian northwest: distribution, number, comparative ecology. Ermala, A., Lahti, S. (Eds.), The 3rd Nordic Beaver Symposium. 10–16. Danilov, P.I., 2009. New mammals in the Russian European North. Karelian Scientific Center of the RAS, Petrozavodsk (in Russian). Danilov, P.I., Fedorov, F.V., 2015. Comparative characterization of the building activity of Canadian and European beavers in the north of European Russia. Russian Journal of Ecology 46 (3), 272–278. https://doi.org/10.1134/S1067413615030029.

14

Ecological Modelling 409 (2019) 108763

V.G. Petrosyan, et al.

Wild furbearer management and conservation in North America. Ashton-Potter Limited, Concord, Ontario, Canada, pp. 283–312. Parker, H., Nummi, P., Hartman, G., Rosell, F., 2012. Invasive North American beaver Castor canadensis in Eurasia: a review of potential consequences and a strategy for eradication. Wildlife Biology. 18, 354–365. https://doi.org/10.2981/12-007. Payne, N.F., 1989. Population dynamics and harvest response of beaver. Craven, S.R. (Ed.), Proceedings of the Fourth Eastern Wildlife Damage Control Conference. USDA Animal and Plant Health Inspection Service, Animal Damage Control 33, 127–134. Petrosyan, V.G., Golubkov, V.V., Goryainova, Z.I., Zav’yalov, N.A., Al’bov, S.A., Khlyap, L.A., Dgebuadze, Yu.Yu., 2013. Modeling of the Eurasian Beaver (Castor fiber L.) population dynamics in the basin of a small Oka river tributary, the Tadenka River (Prioksko Terrasnyi Nature Reserve). Russian journal of biological Invasions. 4, 45–53. https://doi.org/10.1134/S2075111713010086. Petrosyan, V.G., Golubkov, V.V., Zavyalov, N.A., Goryainova, Z.I., Dergunova, N.N., Omelchenkо, A.V., Bessonоv, S.A., Albov, S.A., Marchenko, N.F., Khlyap, L.A., 2016. Patterns of population dynamics of Eurasian beaver (Castor fiber L.) after reintroduction into nature reserves of European part of Russia. Russian journal of biological Invasions. 7, 355–373. https://doi.org/10.1134/S2075111716040068. Ogureeva, G.N., Miklyaeva, I.M., Safronova, I.N., Yurkovskaya, T.K., 1999. The Map “Plant Zones and Their Types in Russia and Adjacent Territories,” 1: 8000000, Moscow. EKOR-M (in Russian). Portenko, L.A., 1926. Beavers (Castor fiber L.) on r. Teterev in Kiev Province. Bulletin of Mosсow Society of Naturalists. Biological series 35, 18–20, (in Russian). Romanko, V.K., 2012. Course of difference equations. Fizmatlit, Moscow (In Russian). Rosell, F., Parker, H., 1995. Forvaltning av bever: dagens tilstand og fremtidig behov. Beaver management: present practice and Norway’s future needs. Telemark University College, Boi Telemark, Norway, pp. 52–55 (in Norwegian). Ruusila, V., 1997. Kanadanmajava on ahkerampi rakentaja. In: Nummi, P. (Ed.), Eds.), Suomen luonto. Nisäkkäät. Weilin & Göös. Porvoo, Finland, pp. 95 (in Finnish). Saveljev, A.P., 2002. Hermann Pohles Biber aus Westsibirien. Mammal Biology. 67, 33 (SH). Saveliev, A.P., Shtubbe, M., Shtubbe, A., Putintsev, N.I., Oleynikov, A.Yu., Saveliev, A.A., 2010. Movements of beavers in a natural setting and in places of introduction. Herald of hunting 7, 340–344 (in Russian). Segal, A.N., Orlova, S.A., 1961. The appearance of beavers in Karelia. Zool. Zhurnal. 40, 1580–1583 (in Russian). Serebrennikov, M., 1929. Review of the beavers of the Palaearctic Region (Castor, Rodentia). Proceedings of the Academy of Sciences USSR 271–276 (in Russian). Sieber, J., Bauer, K., 2001. Europäischer und Kanadischer Biber Castor fiber Linnaeus 1758, C. canadensis Kühl 1820. In: Spitzenberger, F. (Ed.), Die Säugetierfauna Ӧsterreichs, Grüne Reihe des BMLFUW 13. Austria medien service GmbH, pp. 366–374 (in German). Stubbe, M., Dawaa, N., 1986. Die autochthone zentralasiatische Biberpopulation. Zoologische Abhandlungen Staatliches Museum fur Tierkunde in Dresden. 41, 91–103. Stubbe, M., Dawaa, N., Heidecke, D., 1991. The autochtonous central Asiatic beaver population in the Dzungarian Gobi. In: McNeely, A., Neronov, V.M. (Eds.), Mammals in the Palaearctic desert: status and trends in the Sahara–Gobian Region. Man and the Biosphere UNESCO Programme, Moscow, pp. 258–268. Stubbe, M., Romashov, V.A., 1992. Zum Gedenken an den russischen Biberforscher Leonid Sergeevic Lavrov (1911–1992). Semiaquatische Säugetiere. Wiss. Beitr. Univ., Halle, pp. 465–467. Tumanov, I.L., 2009. Rare predatory mammals of Russia (small and medium species). OOO Branko, St. Petersburg (in Russian). Wauters, L.A., Gurnell, J., Martinoli, A., Tosi, G., 2002. Interspecific competition between native Eurasian red squirrels and alien grey squirrels: does resource partitioning occur? Behavioral Ecology and Sociobiology. 52, 332–341. https://doi.org/10.1007/ s00265-002-0516-9. Wright, J.F., Gurney, W.S.C., Jones, C.G., 2004. Patch dynamics in a landscape modified by ecosystem engineers. Oikos. 105, 336–348. https://doi.org/10.1111/j.0030-1299. 2004.12654.x. Zaikin, A.G., 1959. Beavers in the Leningrad Oblast Vol. 6. Hunting and hunting economy, pp. 23 (in Russian). Zavyalov, N.A., Krylov, A.V., Bobrov, A.A., Ivanov, V.K., Dgebuadze, Yu.Yu., 2005. Impact of the European beaver on small river ecosystems. Nauka, Moscow (in Russian). Zavyalov, N.A., Albov, S.A., Petrosyan, V.G., Khlyap, L.A., Goryaynova, Z.I., 2010. Invasion of ecosystem engineer – European beaver (Castor fiber L.) in the Tadenka River basin (Prioksko-Terrasnyi Nature Reserve). Russian Journal of Biological Invasions. 1, 267–281. https://doi.org/10.1134/S2075111710040041. Zavyalov, N.A., Artaev, O.N., Potapov, S.K., Petrosyan, V.G., 2015. Beavers (Castor fiber L.) of the Mordovskii Nature Reserve: Population Development History, Current State, and Prospects. Russian Journal of Biological Invasions. 6, 148–164. https://doi. org/10.1134/S2075111715030091. Zavyalov, N.A., 2015. Ecosystem engineering of the beaver (Castor fiber L.) in the forest zone of European part of Russia. Proceedings of State Nature Reserve «Rdeysky». V. 3 2015 (in Russian with English abstracts). Zavyalov, N.A., Albov, S.A., Khlyap, L.A., 2016. Mobility of settlements and elements of the biological signaling field of Beavers (Castor fiber) in the basin of the Tadenka River (Prioksko-Terrasny Nature Reserve). Biology Bulletin. 43, 1099–1109. https:// doi.org/10.1134/S1062359016090107. Zurowski, W., 1980. Der Europäischer Biber in Poland. Przeglad hodowlany Vol. 48. pp. 18–24 (in German).

Proceedings on forest experimental work 75–140 (in Russian). Halley, D.J., Rosell, F., 2002. The beaver’s reconquest of Eurasia: status, population development and management of a conservation success. Mammal Review. 32, 153–178. https://doi.org/10.1046/j.1365-2907.2002.00106.x. Hartman, G., 2003. Irruptive population development of European beaver (Castor fiber L.) in southwest Sweden. Lutra. 46, 103–108. Heidecke, D., Hörig, H., 1986. Bestands-und Schutzsituation des Elbebiber s. Naturschutzarbeit Halle/Magdeburg 23, 2–14 (in German). Henry, D.B., Bookhout, T.A., 1969. Productivity of beavers in northeastern Ohio. Wildl. Manag. 33, 927–932. Herr, J., Schley, L., 2009. Barbed wire hair traps as a tool for remotely collecting hair samples from beavers (Castor sp). Lutra. 52, 123–127. Ivanov, P.D., 1975. The North American beaver on the Karelian Isthmus of Leningrad oblast. Proceedings of the Voronezh State Nature Reserve 21 (1), 114–120 (in Russian). Kanshiev, V.Ya., 1983. Materials on reproduction of Canadian beavers in Karelia. Fauna and ecology of birds and mammals of the North-West of the USSR. pp. 122–126 Petrozavodsk (in Russian). Kataev, G.D., 2011. Beaver’s Castor fiber in the northern periphery of the range (Kola Peninsula). Bulletin of the Moscow Society of Naturalists, vol. 130. Biological department, pp. 3–11 (in Russian). Kataev, G.D., 2015. Long-term observations over reintroduced beavers Castor fiber orientoeuropaeus on Kola Penninsula, NW Russia. Beavers – from genetic variation to landscape – level effects in ecosystems. 7-th international Beaver Symposium. Book of Abstracts. Biomik Active, Voronezh. 35. Kauhala, K., Turkia, T., 2013. Majavien elinympäristönkäyttö: alkuperäislajin ja vieraslajin alustavaa vertailua. Suomen Riista. 59, 20–33 (in Finnish). Khlyap, L.A., Warshavsky, A.A., Bobrov, V.V., 2011. Diversity of Alien Mammal Species in Different Regions of Russia. Russian journal of Biological Invasions. 2, 293–299. https://doi.org/10.1134/S2075111711040059. Korablev, N., Puzachenko, Y., Zavyalov, N., Zheltukhin, A., 2011. Long-term Dynamics and Morphological Peculiarities of Reintroduced Beaver Population in the Upper Volga Basin. Baltic Forestry. 17, 136–147. Kudryashov, V.S., 1975. Аbout the factors regulating of the number of the European beaver in the Oka Reserve. Mammals. Number, its dynamics and factors determining them. Proceedings of the Okskiy Nature Reserve 11, 5–124 (in Russian). Lathi, S., Helminen, M., 1974. The beaver Castor fiber (L.) and C. canadensis (Kuhl) in Finland. Acta Theriologica. 19, 177–189. Lavrov, L.S., 1956. Revisiting characteristics of population of beavers from the Voronezh Reserve. Proceedings of the Voronezhsky State Nature Reserve 5–11 (in Russian). Lavrov, L.S., 1960. Revisiting biological and morphological differences between the European and the Canadian beavers. Proceedings of the Voronezh state nature reserve 103–120 (in Russian). Lavrov, L.S., 1969. A new subspecies of the Eurasian beaver (Castor fiber) from the Jenisej upper flow. Zoologicheskii Zhurnal. 48, 456–457 (in Russian). Lavrov, L.S., 1981. Beavers of Palearctic. Voronezh University Press, Voronezh (in Russian). Lavrov, L.S., 1983. Evolutionary developments of the genus Castor and taxonomy of the contemporary beavers of Eurasia. Acta Zoologica Fennica. 174, 87–90. Lavrov, L.S., Lavrov, V.L., 1986. Verteilung und Anzahl urspriinglicher und aboriginer Biberpopulationen in der USSR (Distribution and number of aboriginal origin beaver populations in the USSR). Zoologische Abhandlungen Staatliches Museum fur Tierkunde in Dresden Vol. 41. pp. 105–109 (in German). Lavrov, L.S., Lu, H.T., 1961. Present conditions and ecological peculiarities of beavers (Castor fiber L.) in natural colonies in Asia Vol. 9. Vestnik Leningradskogo Universiteta, pp. 72–83 (in Russian). Lavrov, L.S., Orlov, V.N., 1973. Karyotypes and taxonomy of modern beavers (Castor, Castoridae, Mammalia). Zoologicheskii Zhurnal. 52, 734–742 (in Russian). Liukko, U.M., Henttonen, H., Hanski, I.K., Kauhala, K., Kojola, I., Kyheröinen, E.M., Pitkänen, J., 2016. The 2015 Red List of Finnish Mammal Species. Ministry of the Environment and Finnish Environment Institute, Helsinki. Ruusila, V., Ermala, A., Hyvaèrinen, H., 2000. Costs of reproduction in introduced female Canadian beavers (Castor canadensis). Journal of Zoology, Lond. 252, 79–82. Logofet, D.O., Evstigneev, O.I., Aleinikov, A.A., Morozova, A.O., 2015. Succession caused by beaver (Castor fiber L.) life activity: I. What is learnt from the calibration of a simple Markov model. Biology Bulletin Reviews. 5, 28–35. https://doi.org/10.1134/ S2079086415010053. Maran, T., Henttonen, H., 1995. Why is the European mink (Mustela lutreola) disappearing – a review of the process and hypotheses. Annales Zoologici Fennici. 32, 47–54. Martinoli, A., Bertolino, S., Preatoni, D.G., Balduzzi, A., Marsan, A., Genovesi, P., Tosi, G., Wauters, L.A., 2010. Headcount 2010: The multiplication of the grey squirrel introduced in Italy. Hystrix–Italian Journal of Mammalogy. 21, 127–136. Monitoring and conservation of biodiversity of taiga ecosystems of the European North of Russia, 2010. In: Danilov, P.I. (Ed.), Monitoring and conservation of biodiversity of taiga ecosystems of the European North of Russia. Institute of Biology of the Karelian Research Center of the Russian Academy of Sciences, Petrozavodsk (in Russian). Müller-Schwarze, D., Sun, L., 2003. The beaver. Natural History of a wetlands engineer. Comstock Publishing Assotiates, a division of Cornell University Press, Ithaca and London. Nolet, B.A., Rosell, F., 1998. Comeback of the beaver Castor fiber: An overview of old and new conservation problems. Biological Conservation. 83, 165–173. https://doi.org/ 10.1016/S0006-3207(97)00066-9. Novak, M., 1987. Beaver. In: Novak, M., Baker, J.A., Obbard, M.E., Malloch, B. (Eds.),

15