Why the population of the northern fur seals (Callorhinus ursinus) of Tyuleniy Island does not recover following the harvest ban: Analysis of 56 years of observation data

Ecological Modelling 363 (2017) 57–67

Contents lists available at ScienceDirect

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

Research Paper

Why the population of the northern fur seals (Callorhinus ursinus) of Tyuleniy Island does not recover following the harvest ban: Analysis of 56 years of observation data Oksana L. Zhdanova a,b,∗ , Alexey E. Kuzin c , Elena I. Skaletskaya d , Efim Ya. Frisman d a

Institute for Automation and Control Processes, Far Eastern Branch, Russian Academy of Sciences, Vladivostok 690041, Russia Far Eastern Federal University, Vladivostok 690950, Russia c Pacific Research Fisheries Center, Vladivostok 690950, Russia d Institute for Complex Analysis of Regional Problem, Far Eastern Branch, Russian Academy of Sciences, Birobidzhan 679016, Russia b

a r t i c l e

i n f o

Article history: Received 21 March 2017 Received in revised form 15 August 2017 Accepted 25 August 2017 Keywords: Callorhinus ursinus Mathematical modeling Harvest Conservation Survival rates estimation

a b s t r a c t Following many years of managed harvest, the population size of Northern fur seal on Tyuleniy Island appeared to be depressed. In order to preserve the population, the harvest was first significantly limited, and then banned altogether. However, the birth rates have never recovered. We perform the estimation and analysis of survival rates of different groups of males in the population with the goal of elucidating the underlying changes in intra-population dynamics. We use the estimated numbers of adult males and pups in the rookery along with the exact numbers and ages of animals, harvested over the years. Our investigation shows that the unexpectedly slow growth of the population is most likely not caused by an abrupt decline of survival rates of any age groups. There is no substantial change in juvenile survival, and the survival rates of the older animals are actually increasing. However, the analysis of densitydependent factors of juvenile survival demonstrates significant increase in intra-species competition. This may be caused by change in environmental conditions such as quality and availability of food. Also, the dynamics of pups’ viability is telling: it was the highest at the beginning of the period of observations, and then it sharply decreased following intensification of harvest, and finally started growing again in response to strong conservation measures, but never achieved its original level. These facts may be consequence of changes in the population gene pool, which shows the importance of research of the evolutionary consequences of harvest. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Human activity affects to varying degrees a great number of natural plant and animal populations, many of which are now endangered (Buddle et al., 2006; Caputo et al., 2005 etc.). This especially refers to populations of valuable commercial species. The problem of conservation and optimal exploitation of the commercial species is not new. The theory of optimal harvest strategy was developed long ago (Chapman and Robson, 1960; Clark, 1990; Graham, 1935; Hjort et al., 1933; Ricker, 1958; Russell 1931), and aims to provide the maximum sustainable yield with a guaranteed conservation of the population. This theory has been further

∗ Corresponding author at: Institute for Automation and Control Processes, Far Eastern Branch, Russian Academy of Sciences, Vladivostok 690041, Russia. E-mail addresses: [email protected] (O.L. Zhdanova), [email protected] (A.E. Kuzin), [email protected] (E.Ya. Frisman). http://dx.doi.org/10.1016/j.ecolmodel.2017.08.027 0304-3800/© 2017 Elsevier B.V. All rights reserved.

developed in recent years (e.g., Frisman et al., 2010, 2006; Hauser et al., 2006; Horan and Bulte, 2004; Neverova et al., 2016; Skonhoft et al., 2012), and has been abundantly applied in practice (e.g., Abakumov et al., 2016, 2009; Abakumov and Kucher, 1996; Hauser et al., 2007; Martin et al., 2010; Punt et al., 2014; Tahvonen, 2008). However, many populations began to degrade down to extinction even with presumably optimal harvest. Unfortunately, very often simple cessation of exploitation is not sufficient to restore the original population size. There arises the question of which particular intra-population changes this inability to recover can be attributed to. The answer is far from obvious. In this study, we analyze this problem using the example of the local northern fur seal population, which has been exploited and closely observed for a long period of time. In the middle of the last century, the northern fur seal became a unique object of population studies. The high value of its fur and the specifics of its habitat led to the creation of a four-way international convention (Canada, Japan, USSR, USA) on the conservation of

58

O.L. Zhdanova et al. / Ecological Modelling 363 (2017) 57–67

In order to elucidate and hopefully correct the changes in the population of the northern fur seal of Tyuleniy Island that led to slow growth even in absence of harvesting, we need to solve the following problems: 1. Estimate the model parameters: analyze the dynamics of the survival rates of different age and sex groups of the population. One possible obstacle to restoring the population size may be substantial decline of the juvenile (or some other age group) survival rate. 2. Improve the general model of local population dynamics: a) Estimate the effect of density-dependent factors on the juvenile survival rate. It is possible that change in extrinsic conditions (e.g. food quality) led to the increase of intraspecific competition and the resulting decline of the survival rate and numbers of young animals. b) Analyze the effect of sex ratio on pregnancy rates. The important question is whether the harvest had selective consequences by removing the most reproductively active bulls from the population. 3. Get more accurate predictions of population dynamics and improve the harvesting strategy.

Fig. 1. Field counts of bulls (a) and pups born (b) on Tyuleniy Island.

Northern Pacific fur seals (Interim Convention of the Conservation of North Pacific Fur Seals, 1957). Under this Convention, large-scale studies of northern fur seal biology were organized. They included annual estimation of the sizes of age and sex groups of these animals on every large breeding rookery, as well as determining the conditions and sizes of acceptable harvest. During almost 30 years of this convention, a unique data set on population dynamics of this species has been created. These data served as a good basis for estimating major population parameters and for developing various mathematical models of population dynamics. In particular, Frisman et al. (1982) suggested a detailed model of northern fur seal population dynamics and developed methods for model parameters’ estimation. The model coefficients were calculated using the data from more than twenty years of observations, and the resulting model allowed sufficiently good prediction of population size and estimation of the maximum sustainable yield. A large body of additional information was accumulated in subsequent years by the Pacific Research Fisheries Center, resulting by now in the data on the population of the northern fur seal of Tyuleniy Island over 56 years (Kuzin, 2015, 2014, 2010, 2003). This new data enable the evaluation of the utility of the existing model and the methods of estimating its parameters, and allows us to adapt these methods to recent changes in population processes. After many years of managed harvest (and despite all attempts to optimize it), the population of the northern fur seal of Tyuleniy Island began dropping (Kuzin, 2010, 1999; etc.). By the late 1980s, the number of newborn pups dropped by half compared to the average numbers during the 1960s and early 1970s, and then stabilized at this low level. Other fur seal populations of the North Pacific also suffered reduced populations, particularly Commander and Pribilof Islands (Kornev et al., 2008; Trites and Larkin, 1989; York, 1987; and others). In order to reverse this negative trend, the harvest on Tyuleniy Island was substantially reduced by the 1990s, and then completely banned in 2008. As expected, the number of adult males dramatically increased (Fig. 1a), but the anticipated increase in the number of births in the population did not follow: although the numbers of newborn pups are growing, the level of the mid-1960s was never achieved (Fig. 1b).

In this work, we focus on the first problem, limited to the male part of the population. We use the data on the numbers of bulls and pups in the rookery together with the age composition of on-land harvest (for the years preceding the harvest ban). During most of the considered period of time, only subadult males were harvested. It made the estimation and prediction of the dynamics of these groups of animals essential for rational management of the fur seal populations. Direct visual observations allow us to estimate the numbers of only two parts of the population, namely pups and bulls. It makes the methods of estimating the survival rates for different age groups essential for such prediction. There have been many attempts to determine these rates, and a large body of literature has ensued (Andreev et al., 1978; Chapman, 1964; Corey et al., 2003; Ichihara, 1972; Lander, 1975; Trites, 1989; etc.). Numbers of newborn pups, bulls, and harvest information were used as the input in most of these works. However, there are substantial differences between the estimates by different authors. Our goal here is to get the estimates that accord best with our dataset. Solving this problem should help us answer the question of whether the population does not rebound because the survival rates of juveniles or other age groups decreased. It may also clarify the influence of population density on juvenile survival rates and possible selective effects of past harvests.

2. Material and methods 2.1. Description of observations results Reliable direct counting of fur seals of all age and sex groups is not possible: the animals in the rookeries are very densely grouped and are constantly moving; their numbers are high; and they are spread over a large area of often forbidding terrain. Only adult males (bulls) can be counted directly. They are much bigger than others, their total number is relatively small, and during the reproductive period their movements are minimal. It is easy to identify them in harems, and usually they are at considerable distances from each other. The numbers of bulls based on such visual estimates are available for the entire period of observations (1958–2013). The second group of animals allowing direct count are the pups (their numbers are also available for almost all years between 1958

O.L. Zhdanova et al. / Ecological Modelling 363 (2017) 57–67

and 2013). Counting them is possible because they do not leave land during their first weeks of life. Finally, the exact numbers and ages of the animals harvested on land are also recorded for this period of time. The commercial harvest was limited to subadult males two to five years of age, with the majority being 3–4 years old (while northern fur seals can live to be 20–25 years old). Respectively, the age composition of the commercial harvest does not allow us to directly determine the age distribution of all males in the population. Besides, there is no harvest data since the harvesting ban of 2008. To summarize, we have directly estimated numbers of bulls and pups, and the numbers and ages of harvested animals (Tables 1 and 2). 2.2. Mathematical model of the numbers of males in the local population As the most complete data are available for the male part of the population, we will write the equations describing male numbers only. Let P denote the number of newborn pups, M0 – the number of newborn males, Mi – the number of subadult males of age i (i = 2, 6); M – number of bulls (adult males age 7 and up). All these numbers, changing from year to year, are the variables in our model. As the natural cycle of the population is one year, so is the step in our model. The numbers of all these population groups change annually because a new generation of pups is born, and all animals move from one age group to the next (taking into account the natural death rate and harvest). We get the following equations: M 0 (n) = P(n)/2

(1)

M 2 (n) = w02 (n − 2)M 0 (n − 2)

(2)

M i+1 (n) = wii+1 (n − 1)(M i (n − 1)–Ri (n−1)),i = 2, 5;

(3)

M(n) = w67 (n − 1)(M 6 (n − 1)–R6 (n − 1)) + w(n − 1)(M(n − 1) –R(n − 1))

59

2.3. Age-specific survival rates of the northern fur seal The survival rates of fur seal males from birth to two years of age were estimated using a modified Lander’s method (Lander, 1975), and using a new method based on Trites’ approach (Trites, 1989). For all other age groups, the survival rates were calculated using the procedure suggested by Frisman et al. (1985). The quality of the resulting parameters was estimated when calculating bull numbers. For the populations of large mammals (like pinnipeds), juvenile survival is a key parameter defining population dynamics and an important indicator of the population condition (Eberhardt, 1981). There are few methods for estimating this parameter. In particular, besides the most commonly used Lander’s method (Lander, 1975), there are also Chapman’s (Chapman, 1964) and Smith and Polacheck (Smith and Polacheck, 1984) methods, which only determine the lower bound of the male survival rate from birth to three years of age, and are based on certain assumptions on harvest structure. The reliability of all three methods is studied in details by Trites (1989), with the most critique directed to Lander’s method due to correlation of his estimates with the harvest size. We do not completely abandon Lander’s method, because it is biologically grounded and produces adequate results if the requirements for the harvesting strategy are met. We do not have any reasons to believe that this assumption was frequently wrong for the Tyuleniy population, which was intensively harvested up to 2008. We supplement these estimates with the new ones, based exclusively on biological arguments without any assumptions about harvest, and then compare the results. All methods we use produce both upper and lower bounds of the juvenile survival rates: w02L and w02U , and we accept their average w02 = (w02L + w02U )/2 as the final estimate. 2.3.1. Estimation of the survival rate from birth to 2 years of age using basic and modified lander’s method 2.3.1.1. Basic method. Let us consider the basics of Lander’s method (1975), using Eqs. (2) and (3) of our model. Using the Eq. (3), we get:

(4)

where w02 (n) is the survival rate of fur seal males born in year n during two first years of life, wii+1 (n) – survival rate of males being i years old in the nth year between i and i + 1 year of age, w(n) – survival rate of bulls between nth and (n + 1)th years, Ri (n) – number of harvested males of age i (i = 2, 6) in year n, R(n) – number of harvested bulls in year n. The choice of the model’s variables is defined by data availability. Reliable information about the ages of bulls is absent, as is any information about the numbers of yearlings (they have never been harvested and cannot be counted visually because they spend most of their time in the water). Therefore, aggregated values are used: the group including all bulls (ages 7 and up), and the survival rate of males from birth to two years of age. The validity of the Eq. (1) is confirmed by numerous studies: sex ratio of fur seal embryos is very close to 1:1 (York, 1987), and although by birth there may be some deviation from 1:1 ratio both ways in different populations (Antonelis et al., 1994; Boltnev, 2011; Fauler, 1998; etc.), these deviations are insignificant. Model (1–4) does not take into account the migration of fur seals between populations. It does not affect the validity of the model for the Tyuleniy population, because a previous detailed study of migration flows (based on tag returns) showed a very high level of isolation of this population from all others (Frisman et al., 1985). The inflow of Commander and Pribilof fur seals is negligible, and local fur seals in turn are rarely spotted near other islands. Let us now consider the methods of calculation of the model coefficients based on the available data.

M 5 (n) = w45 (n − 1)(M 4 (n − 1) − R4 (n − 1)) = w45 (n − 1)((w34 (n − 2)(M 3 (n − 2) − R3 (n − 2)) − R4 (n − 1)) = w45 (n − 1)(w34 (n − 2)(w23 (n − 3)(M 2 (n − 3) − R2 (n − 3)) − R3 (n − 2)) − R4 (n − 1)) Substituting M2 in the last expression, we get: (∗)M 5 (n) = w45 (n − 1)(w34 (n − 2)(w23 (n − 3)(w02 (n − 5)M 0 (n − 5) − R2 (n − 3)) − R3 (n − 2)) − R4 (n − 1)) According to Lander, the yearly survival rate of males from 2 to 5 years of age is constant and equal to w25 : w45 (n − 1) = w34 (n − 2) = w23 (n − 3) = w25 (n − 3) Substituting it into (Equation *) after some algebra we get (∗∗)M 5 (n) = (w25 (n − 3))3 w02 (n − 5)M 0 (n − 5) − (w25 (n − 3))3 R2 (n − 3) − (w25 (n − 3))2 R3 (n − 2) − (w25 (n − 3))R4 (n − 1)) Based on this formula, with some additional assumptions about the size and structure of harvest, we can get interval estimates of the juvenile survival rate. It is obvious that the number of fiveyear-old males before the harvest was no less than the number of

60

O.L. Zhdanova et al. / Ecological Modelling 363 (2017) 57–67

Table 1 Field counts of bulls and pups born on Tyuleniy Island. Year

Number of bulls

Number of pups born

Year

Number of bulls

Number of pups born

1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985

1016 1360 2279 2528 2628 2465 2452 2213 2934 2439 2242 2004 2331 1864 1147 564 540 553 738 861 1285 1603 2000 1833 1405 1622 978 875

32 200 35 000 38 000 41 200 44 700 49 000 51 400 48 400 44 900 56 500 45 800 43 500 31 500 41 100 44 050 35 300 33 170 27 000 30 828 28 580 29 900 24 185 22 000 21 000 22 300 24 000 23 200 20 000

1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013

990 1007 1604 2124 2503 2721 2868 2739 3035 2812 2588 1990 3153 3169 4335 2896 4014 3980 4200 3137 3785 3460 5827 4781 5346 4958 5287 5615

18 000 16 000 18 800 18 500 18 200 18 000 15 000 17 000 18 700 17 600 16 500 18 500 19 200 20 650 23 630 24 000 26 400 27 000 32 400 33 500 33 500 36 000 38 000 41 400 39 300 35 500 35 100 34 700

harvested five-year-old males, so the lower bound of their number is M5L = R5 . Substituting it into (Equation **) and solving for w02 , we get the lower bound for w02 : w02L (n − 5) = +

R2 (n − 3) R4 (n − 1) −2 R3 (n − 2) −1 + w + w M0 (n − 5) M0 (n − 5) 25 M0 (n − 5) 25

R5 (n) w−3 . M0 (n − 5) 25

(5)

Lander assumed that the upper bound was M5U (n) = R5 (n) + R4 (n − 1). This is apparently true for the generations from which at least half of 4 year old males was harvested: R4 (n) ≥ M4 (n)/2. Substituting M5U (n) into (Equation **) and solving for w02 , we get the upper bound of the survival rate from 0 to two years (w02U ): w02U (n − 5) = +

R2 (n − 3) R4 (n − 1) −2 R3 (n − 2) −1 + w + w M0 (n − 5) M0 (n − 5) 25 M0 (n − 5) 25

R5 (n) + R4 (n − 1) −3 w25 . M0 (n − 5)

calculated and inequality 02L ≥ 25 is verified. The smallest w25 for which this inequality is true is accepted as the lower bound w25L . When w25L is calculated, and assuming w25U = 1, we can get the following set of estimates: M5L = R5 ,

M5U = R4 + R5 ,

M4L = R5 + R4 ,

M4U = M5U /w25L + R4 ,

M3L = R5 + R4 + R3 ,

M3U = M4U /w25L + R3 .

Here, two upper equalities mirror Lander’s assumptions about the harvest, and the others are their direct consequence. ˆ 25 = These formulae allow to estimate w25 . Let us assume w ˆ 4 /(M ˆ 3 − R3 ), M ˆ 4 = 0.5(M4U + M4L ), M ˆ 3 = 0.5(M3U + M3L ). SubstiM tuting M4U , M4L , M3U , M3L from (Eq. (7)) we get: w25 =

(6)

Now we need to estimate w25 . This requires additional assumptions and calculations. We can assume with confidence that the natural death rate of fur seals between birth and 2 years of age is higher than that of males between two and five years of age. The quantitative measure of the death rate is the “instant death rate” (␮), related to the survival rate (w) by the formula  = (−ln w)/t, where t is the analyzed time interval. R. Lander suggests using one month as the basis, and respectively to estimate male death rate as 02 = (−lnw02 )/24. Further it is assumed that the natural death rate during the harvest is zero. This assumption is based on the rarity of finding dead young males in the rookery. Thus, the denominator for calculation of yearly survival rate is 11 rather than 12 months: 25 = (−lnw25 )/11. Taking into account that 02L = (− ln w02U )/24 ≥ 25 = (−lnw25 )/11 and (Eq. (6)) we can fit the lower bound w25L . Trial values of w25 decrease (starting with w25 = 1) with small step (equal to 1E-6), and for each value respective w02U is

(7)

ˆ4 M ˆ 3 − R3 M

= ... =

−1 ) R4 + (R4 + R5 )(1 + w25L −2 −1 ) + R4 w25L (R4 + R5 )(1 + w25L

.

(8)

We accept w25 from (Eq. (8)) as its estimate. Substituting it into (Eq. (5)) and (Eq. (6)), we can get w02L and w02U . 2.3.1.2. Modification 1. For the periods of time when mostly 3 year old males were harvested, Frisman et al. (1982) developed the following modification of Eqs. (5)–(6) using assumptions (R3 (n) ≥ M3 (n)/2) and M4U = R4 + R3 : w02L = w02U = where −1 R3 w24L ),

R2 R3 −1 R4 −2 + w + w M0 M0 24 M0 24

R2 R3 −1 (R4 + R3 ) −2 + w + w24 M0 M0 24 M0

,

(9)

−1 −2 w24 = (R3 + (R4 + R3 )(1 + w24L ))/((R4 + R3 )(1 + w24L )+

and the lower bound w24L is obtained similar to w25

before.

2.3.1.3. Modification 2. The generation of fur seals born in 1992 was not harvested at the age of three years or at the age of four (R4 = R3 = 0), which means that the two previous methods are not

O.L. Zhdanova et al. / Ecological Modelling 363 (2017) 57–67

61

Table 2 Catchs of mail fur seals on Tyuleniy Island. Year

Estimated composition by age class 2

3

4

5

6

7

8

9

10

10 +

1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008

89 88 124 340 480 459 1032 1336 2238 1014 1486 1833 874 613 1657 935 556 508 469 572 794 530 659 721 982 592 453 356 653 0 516 206 424 486 455 355 375 0 0 121 212 340 286 387 285 220 104 273 380 0 593

1386 1976 2116 3221 3080 3387 4320 4678 4199 3582 1341 3360 4580 3525 2975 2609 1750 1757 1737 2911 1911 1719 1755 2396 1516 1206 1521 103 1287 0 1132 842 789 788 950 862 862 0 0 769 612 665 562 731 849 770 811 774 1057 0 905

858 1795 1670 1584 1711 1734 1852 1869 1598 1839 811 575 1026 1662 786 544 185 239 341 507 448 560 567 289 336 244 292 0 99 0 120 282 292 204 301 324 222 0 0 526 353 370 248 337 493 378 279 417 514 0 344

566 554 645 460 496 684 379 450 557 1119 573 242 275 517 362 30 11 6 22 10 47 91 19 94 31 22 28 0 4 0 6 14 43 33 37 127 31 0 0 73 167 111 41 117 128 89 23 80 125 0 30

89 137 391 227 214 205 52 77 139 594 420 93 127 200 107 10 5 0 1 0 0 0 0 0 5 1 3 0 3 0 0 5 1 7 7 23 6 0 0 6 32 8 8 18 13 8 1 8 17 0 14

41 87 463 354 273 226 73 13 110 255 175 34 96 70 13 3 4 0 0 0 0 0 0 0 1 0 4 0 1 0 2 7 1 3 13 6 0 0 0 0 74 1 2 3 11 1 0 9 8 0 0

42 22 356 353 306 249 166 3 157 141 105 28 42 53 2 0 4 0 0 0 0 0 0 0 0 2 2 0 0 0 0 4 2 1 27 6 0 0 0 0 80 1 3 1 6 1 3 7 1 0 0

0 0 166 237 231 156 162 3 113 121 83 18 9 69 1 0 1 0 1 0 0 0 0 0 0 2 2 0 0 0 0 1 2 2 35 9 0 0 0 0 43 1 2 0 3 1 1 5 0 0 0

0 0 94 67 74 60 50 2 46 74 45 7 4 87 1 0 10 0 0 0 0 0 0 0 3 4 0 0 0 0 1 1 2 0 57 6 0 0 0 0 34 1 0 0 1 1 0 3 0 0 0

56 192 184 104 169 105 123 1 26 19 25 31 9 207 6 1 0 0 2 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 68 11 0 0 0 0 60 1 0 1 4 2 0 3 1 0 0

applicable for this year. The most significant harvest from this generation consists of two- and five-year-old animals. It means that the equation M5U = R2 + R5 should be rewritten as M5U = R4 + R5 . The formulae (5)–(6) now have to be rewritten as w02L =

R2 R5 −1 R2 R2 −1 R5 −1 + w and w02U = + w + w M0 M0 25 M0 M0 25 M0 25

(10)

−1 −2 )/(R5 + (R2 + R5 )w25L ). where w25 = (R5 + (R2 + R5 )w25L

2.3.2. Estimation of the survival rate from birth to two years of age without any assumptions about harvest structure In order to get the estimate independent of harvest structure, we modify the Trites’ method (Trites, 1989) for calculation of the upper and lower bounds of juvenile survival (from birth to two years): (a)w02L = (R2 + R3 w27 −1 + R4 w27 −2 + R5 w27 −3 + R6 w27 −4 )/M0 (b)w02U = (R2 + R3 w27 −1 + 2R4 w27 −2 )/M0

(11)

For the upper bound calculation, the survival rate of subadult males (w27 ) is replaced by maximum possible value (w27U = 1 in Eq. (11)a) and minimum possible (w27U = (w02 )0.5 , in Eq. (11)b). Note that the lower bound estimate (w02L ) in the Eq. ((11)a) is based exclusively on the information about the total harvest from this generation (until the group reaches seven years of age). Including the harvest of later ages (six and seven years) certainly improves the estimate quality, because it uses more information than Lander’s method does. However, if the total harvest from the generation is low, then the lower bound w02L also decreases. It is undesirable, because this decrease has nothing to do with the real value of w02 . It is also risky to use the assumption R4 > M4 /2 in (Eq. (11)b) while getting upper bound (w02U ), because if it is not true, then this estimate also unduly decreases. Let us try to correct the estimates (Eq. (11)) considering not only the harvest from each age group of the generation, but also the

62

O.L. Zhdanova et al. / Ecological Modelling 363 (2017) 57–67

information about the number of bulls provided by this generation upon reaching seven years of age (M7 ):

Finally, w27 = (w27U + w27L )/2.

(a)w02L = (R2 + R3 w27 −1 + R4 w27 −2 + R5 w27 −3 + R6 w27 −4 + M7L w27 −5 )/M0

(b)w02U = (R2 + R3 w27

−1

+ R4 w27

−2

+ R5 w27

−3

+ R6 w27

−4

+ M7U (n)w27

−5

)/M0

(12)

M7L (n) = M(n) − w27 (M(n − 1) − R(n − 1)) M7U (n) = M(n)

Note that these estimates are based only on the known harvest from each age group of this generation, rather than on assumptions about the fraction of the age group harvested. This is possible because the estimates of the number of new bulls are available. In order to determine the minimum possible survival rate from birth to two years of age (w02L), let us assume that the survival rates of older ages are 1: wU = w27U = 1 w02L = (R2 + R3 + R4 + R5 + R6 + M7L )/M0 , M7L (n) = M(n) −(M(n − 1) − R(n − 1)).

2.3.4. Average survival rate of the fur seal bulls We use Eq. (14) to estimate the survival rate for the bulls. This equation can be considered as a linear model, and we can estimate at once two unknown parameters (the average survival rate of the bulls w and the average number of the seven years old males M7 ) using the least squares method. We already have other estimates of the numbers of seven-year-old males M7 (n), so at this time we can verify them. In particular, we can rewrite (Eq. (11)) as M(n) − M 7 (n) = w(M(n − 1) − R(n − 1)),

(15)

and expect the intercept of the linear regression (Eq. (15)) to be close to zero. If it is not substantially different from zero, we can assume that our previous estimates were unbiased. 2.4. Some remarks about choice of methods for calculation of age-specific survival rates

Now, let us compare the estimated value of М7L with the total harvest from the generation in subsequent years R6+ (if М7L < R6+ , then we increase М7L to R6+ , because there could not be fewer seven-year-old males than were subsequently harvested). The maximum possible juvenile survival rate (w02U ) corresponds to the minimum possible survival rate of older animals (w27L = (w02U )0.5 , w = 0): w02U = (R2 + R3 (−0.5 + R4 −1 + R5 −1.5 + R6 −2 + M7U (n)−2.5)/M0 , M7U (n) = M(n).

Let’s consider possibility to run forward in time the population model and estimate the initial condition and model parameters by fitting to the observational data. Indeed, the model under consideration can be reduced to a single equation, where the number of bulls (M(n)) depends on the previous year’s bulls (M(n − 1)), pups born a few years ago (M1(n − 7)) and harvested animals (Ri): M(n) = (M(n − 1) − R(n − 1))w27 + (((((w02 M 1 (n − 7) − R2 (n − 5))w27 − R3 (n − 4))w27 − R4 (n − 3))w27 − R5 (n − 2))w27

2.3.3. Calculation of male fur seal survival rate from two to seven years of age We estimated the survival rates from two to seven years of age (w27 ) using the method from (Frisman et al., 1982). Let us recap its basic assumptions. Assume there is a generation of males for which observations until the age of seven are available. The males from this generation join the group of bulls when they reach the age of seven (M7 ). The losses from two to seven years include harvest (R2 –R6 ) and natural death: M 7 = (((((M 2 –R2 )w27 –R3 )w27 –R4 )w27 –R5 )w27 –R6 )w27 .

(13)

If the males from this generation turn seven in the nth year, then the total number of bulls in that year is M(n) = M 7 (n) + w(M(n − 1) − R(n − 1)).

(14)

With the maximum possible survival rate w27 = w27U , the number of new bulls in the nth year is not greater than the total number of bulls M(n) ≥ M7 (n): M(n) = (((((M 2 (n − 5) − R2 (n − 5))w27U − R3 (n − 4))w27U − R4 (n − 3))w27U − R5 (n − 2))w27U − R6 (n − 1))w27U . We can solve this equation for w27U , and if it appears greater than one, then we replace it with w27U = 1. The estimate of minimum possible survival rate w27L can be found as a solution of the system (13)–(14) with the maximum possible bull survival rate (w = wU ). Multiple observations show that the yearly survival rate of the fur seal bulls is lower than the yearly survival rate of subadult males (from two to seven years), i.e. wU = w27L , and M(n) − (M(n − 1) − R(n − 1))w27L = = (((((M 2 (n − 5) − R2 (n − 5))w27L − R3 (n − 4))w27L − R4 (n − 3))w27L − R5 (n − 2))w27L − R6 (n − 1))w27L .

− R6 (n − 1))w27 All the three parameters of this model (w02 , w27 and w) can be estimated statistically (by the maximum likelihood method, for example), the data set is quite large (56-7–3 = 46 ◦ of freedom). Here we assume the survival rates did not change during the whole period of observations, i.e. the model parameters are constant. However, the population dynamics observed indicates some changes occurred and the purpose of the current study is to identify and analyze these changes. So we have to divide the entire observation period into some parts and then make fitting within each series, but it is not clear how should we divide it (e.g., how many periods should be and how long they are). With such fragmentation of entire observation period the number of degrees of freedom decreases considerably and, consequently, the confidence intervals of estimates becomes wider. We use the methods based on those previously developed for fur seals populations. They allows to obtain feasible parameter intervals using only the life-cycle model of this species and some biologically justified assumptions (such as w02 < w27 , w < w27 , w23 = w34 = . . . = w67 = w27 ). So we consider the survival of the selected age groups as varying in time (in fact, for each year of observations) and its dynamics obtained allows drawing conclusions whether the observed changes in the parameters are random or regular. 3. Calculation One of the basic assumptions of Lander’s method is that at least half of four-year-old males (or three-year-old males) are harvested from each generation, i.e. R4 (n) ≥ M4 (n)/2 (or R3 (n) ≥ M3 (n)/2). The available data about the age composition of harvested animals (Ri (n)) do not provide direct estimates of the numbers of animals of a certain age remaining after harvest. We can get some

O.L. Zhdanova et al. / Ecological Modelling 363 (2017) 57–67

63

Fig. 2. Dynamics of the estimates of male survival rate from birth to two years (w02 (n), left) and number of two years old males (M2 (n + 2), right). Values of survival rates w02 l are calculated using Lander’s method and its modifications (5–10). Corresponding estimates of the numbers of two-year-old males are M2 l . Lower and upper bounds of the survival rates (w02L , w02U ) are calculated using Eq. (12), corresponding estimates are w02 2 = (w02L + w02U )/2, M2 L , M2 U , and M2 2 .

information about the number of remaining four-year-old males by summarizing the number of the older animals harvested from the same generation (unfortunately, any information about the age of harvested animals older than 10 years is absent): R4+ (n) =

6 

R4+i (n + i), where R4+ (n) is the number of males being 4 years

i=1

old in the year n, harvested in subsequent years after reaching ages from 5 to 10 years. In this way, it is possible to calculate the number of animals form the given age groups harvested in subsequent years (R4+ (n)) and compare it with the number of animals harvested at four years of age (R4 (n)). If the subsequent harvest is greater (R4+ (n) > R4 (n)), then the main Lander’s assumption is apparently wrong. For the Tyuleniy population, these are the years 1966, 1985, 1987, 1995, 1996, 2007, and 2008, and in some years, there was no harvest at all (1987, 1995, 1996, 2007, and after 2008). It means that for the males born in 1962, 1981, 1983, 1991, 1992, 2003, and 2004, the basic Lander’s method should be replaced with modified method. The applicability of the modified method was verified in a similar manner, by comparing the number of harvested three-year-old animals with the number of harvested males from the same generation in subsequent years: R3+ (n) =

7 

R3+i (n + i), where R3+ (n) is

i=1

the number of males which were three years of age in the year n and were harvested later, until reaching the age of 10. This verification shows that the modified method is applicable for all considered years except 1992. There was no harvest in two consecutive years, 1995 and 1996, and the males born in 1992 were not harvested at the ages of three and four. In order to estimate w02 for this generation, we used modification (Eq. (10)), developed specifically for this case. The new method of w02 estimation based on the Eq. (12) does not require any knowledge of the numbers of animals remaining after the harvest, while the additional estimates of seven-year-old males (М7L and М7U ) are straightforward.

4. Results 4.1. Estimation of model coefficients The comparison of the estimates of juvenile survival rates w02 (n), obtained by the modified Lander’s method and the new method, based on Eq. (12) over the entire observation period, is presented on Fig. 2 left. The corresponding estimated numbers of two-year-old males (M2 (n + 2) = M0 (n)w02 (n)) are presented on

Fig. 2 right. Lower and upper bounds of possible values w02 , w02U and w02L , are calculated using formulae (12). We see that the Lander’s estimate of w02 gravitates to the lower bound (Eq. (12)a), and for some years crosses it. Next figure (Fig. 3) presents the dynamics of the estimates of the coefficient w27 (n) and corresponding estimates of the number of seven-year-old males. M7 (n + 7), calculated using the Eq. (13). Prior to the estimation of regression (Eq. (15)) coefficients, we investigated whether the data are homogenous enough to use the observation over the entire period, and whether there were changes of major trends. Substituting calculated numbers of M7 (n) and observed numbers of M(n), M(n − 1), and R(n − 1), we evaluated the survival rate of bulls for each year: w(n) = (M(n) − M7 (n))/(M(n − 1). The results are presented on Fig. 4. It is obvious that the data are not homogenous, and both sets of estimates show the change of trends that happened in the late 1980s. A likelihood ratio test confirms this change of trend in 1988, p < 0.01. Note that the dynamics of the coefficients calculated based on Lander’s method changed the most. Not only do both mean and variance of w(n) 1 sharply increase in the late period (1989–2011), but the estimated bull survival rates for some years (1998, 2002, and 2008) are greater than one, which means that the calculated number of seven-year-old males for these years is significantly underestimated. Most likely, the illustrated by Fig. 4 dynamics not only shows changes in the population processes, but also demonstrates that after 1988, Lander’s method stopped producing adequate estimates of juvenile survival. The latter can be explained by changes in harvesting. Keeping all these in mind, we further used the estimates based on the Eq. (12). All periods of observation we separated into two intervals: 1965–1988 and 1989–2013. For each of these intervals, we calculated regression (Eq. (15)) coefficients by the least squares method (OLS): 1965–1988гг. :M(n) − M 7 (n) = 0.446·(M(n − 1) − R(n − 1)) − 99.2 (16)

1989–2013гг. :M(n) − M 7 (n) = 0.48·(M(n − 1) − R(n − 1)) + 5.35 (17)

After calculating OLS estimates we tested the errors assumption of the Gauss–Markov theorem (i.e., the errors should be uncorrelated and homoscedastic) to be aware that our estimator is efficient in the class of linear unbiased estimators.

64

O.L. Zhdanova et al. / Ecological Modelling 363 (2017) 57–67

Fig. 3. Dynamics of the estimates of male survival rates from two to seven years (left) and the numbers of seven-year-old males M7 (n + 7) (right). Values w27 calculated using the estimates of w02 1 (basic and modified Lander’s method); w27 2 and M7 2 –using estimates of w02 2 (new method).

Testing for homoscedasticity and autocorrelation (with White Heteroskedasticity and Breusch-Godfrey Serial Correlation LM Tests) shows that residuals of late period (17) are uncorrelated and homoscedastic, but residuals of early period (16) are correlated and heteroscedastic. Therefor we used Newey-West HAC Standard Errors & Covariance for early period (16). Using Shapiro-Wilk and Jarque-Bera tests we found that errors normality cannot be rejected for both periods (16) and (17), so we proceed with using t- and F-statistics. Both equations are statistically significant, with a very high coefficient of determination: R-squared = 0.97 for the early period (Eq. (16)) and R-squared = 0.91 for the late period (Eq. (17)). Both coefficients in the Eq. (16) are statistically significant, which means that the number of seven-year-olds is on average overestimated by 99 animals. The intercept in (Eq. (17)) is not statistically different from zero (p = 0.96), i.e. the estimate of the number of seven-year-olds is unbiased. Fig. 5 presents the regression model (Eq. (15)) for both periods (Eq. (16)), (Eq. (17)). Two major changes happened in the late 1980s and necessitated building separate regression models for two periods of time (Eqs. (16), (17)). First, the survival rate of bulls slightly increased (from 0.446 to 0.48). Second, the observed values became much more dispersed around the theoretical line, while the bias in estimation of the number of seven-year-old males lost its statistical significance; i.e., an estimate of M7 became unbiased! The number of seven-yearolds was overestimated on average by 99 during the early period (1965–1988). 4.2. Modelling of the dynamics of bull numbers Let us evaluate the calculated survival rates of different age groups of fur seal males. The question is whether using these

1

and M7 1 are

Fig. 5. Values of the variables of the regression model (Eq. (15)) and theoretical regression lines separately for each period (Eqs. (16), (17)). Blue circles denote numbers of bulls in the late period (1989–2011); green diamonds denote numbers of bulls in the early period (1965–1988). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

rates produces adequate model dynamics. Using the available data about the numbers of newborn pups (M0 ) and harvest (Ri ) over the entire period of observations, together with calculated values of w02 , w27 , and w, we can calculate the numbers of two years old males: M2 (n) = w02 (n − 2)M0 (n − 2), seven years old males: M7 (n + 5) = (((((M2 (n) − R2 (n))w27 (n) − R3 (n + 1))w27 (n) − R4 (n + 2))w27 (n) − R5 (n + 3))w27 (n) − R6 (n + 4))w27 (n), and bulls M(n) (formulae (13)–(14) with the correction M7 (n) – 99.2 for the early period, because the bias is significant, and without correction for the late period, because the bias is not significant). Fig. 6 presents comparison of modelled (with fixed bull survival rates from Eqs. (16), (17) and observed dynamics of bull numbers. The average error of approximation is 3.2%. It is calculated using the formula

2013 

|M(n) − Mobs (n)|/Mobs (n)/49 · 100%, where M(n)

n=1965

and Mobs (n) are modelled and observed numbers of bulls in the nth year, respectively. 4.3. About density-dependent juvenile survival We now consider the density effects on the juvenile survival rate. Although there is no visible decline in the survival rates of young animals, it is plausible that adverse changes in the environment (like diminished quantity and availability of food) could lead to an increase of intraspecific competition. The simplest formalization of the density dependence of the juvenile survival rate is a linear equation (analogous to Verhulst model (1838)): Fig. 4. Dynamics of estimated bull survival rates. Values w(n) 1 are calculated using estimates w02 1 (basic and modified Lander’s method); w(n) 2 –using estimates w02 2 (new method).

w02 (P) = a − bP.

O.L. Zhdanova et al. / Ecological Modelling 363 (2017) 57–67

65

Fig. 6. Observed and calculated numbers of bulls. Average error is 3.2%.

Fig. 7. The dependence of the juvenile survival rate w02 (n) on the number of newborn pups (P(n)) and theoretical regression line (left). The dynamics of residuals (right).

Estimating the coefficients by the least squares method, we get the following: w02 (P) = 0.51–1.99 ∗ 10−6 P. The coefficient of determination is very small (R-squared = 0.1), which means that only a small part of the w02 variation can be explained by density-dependent factors (Fig. 7, left). However, the regression is significant, and coefficient b is significant also (p < 0.03). Looking at the graph of residuals (Fig. 7, right), one may notice that the data are not homogenous and may be split into at least three distinct periods: 1958–1965, 1966–1985, and 1986–2006. This choice of periods is confirmed by a likelihood ratio test (change of population trends in 1965 and 1985 years, p < 0.001). Below are the estimates of the coefficients, the coefficients of determination (R-squared), and the significance of the slope b (p, Student test) for the three periods of time:

Fig. 8. Dependencies of juvenile survival rates (w02 (n)) on the number of newborn pups (P(n)) and theoretical regression lines for three periods of time.

1958–1964гг. :w02 (P) = 0.75–5.1 ∗ 10−6 P, (R2 = 0.53, p = 0.06) (18)

1964–1985гг. :w02 (P) = 0.46–2 ∗ 10−6 P, (R2 = 0.28, p = 0.01) (19)

1985–2006гг. :w02 (P) = 0.6–5.7 ∗ 10−6 P, (R2 = 0.4, p < 0.01) (20) Fig. 8 presents the same data as Fig. 7 but is split into the three periods. Periods are denoted by the color and shape of the marker,

and the period-specific regression line is shown. Although the three data series are not long (especially the first period, 1958–1964), and the statistical methods are not very reliable, one cannot help but notice the changes in the population dynamics. While the juvenile survival rates increased in the late period, so did the slope b (almost three times compared to 1964–1985). This is a sign of increased intraspecific competition. Note that, before the year 1964, this coefficient was almost as big, so its change can be cyclical, related to cycles in the food supply or other extrinsic factors. The value of the intercept a, which characterizes the juvenile survival without competition (viability), also varied over time. It was the highest in the early period, a little lower in the late period, and the lowest in the middle (1965–1985). Overall, the latest period is characterized by increased intraspecific competition and the general decrease of the viability of pups.

66

O.L. Zhdanova et al. / Ecological Modelling 363 (2017) 57–67

Coefficients of determination in all three cases are modest, which means that there are other factors influencing juvenile survival, which must be included in the analysis.

5. Discussion and conclusions We obtained satisfactory estimates of the coefficients characterizing the major stages of the life cycle of the male fur seals: survival rate between birth and two years of age (w02 ), subadult (w27 ) and adult (w) males. The estimation of the juvenile survival rate is the most important step, and the survival rates of the next age groups are chosen to fit the observed numbers of adult males. If w02 is underestimated, it causes automatic overestimation of the other rates (w27 and w), and vice versa. Lander’s method (including the modifications of the basic method for some years) for w02 estimate and subsequent modelling leads to a quite dramatic scenario of the population dynamics. The survival rates from birth to two years of age notably decrease (Fig. 2, left), and the numbers of two-year-old males sharply decline (Fig. 2, right). The survival rate of adult males after the late 1980s increases significantly, while the variation of this parameter also increases (Fig. 4). Lander’s method depends on certain assumptions about harvest structure (relative share of harvested animals form each age group), which are not verifiable with the available data. In order to avoid this problem, we developed the new method, independent of such assumptions. The formulae (12) allow to calculate for each generation intervals of possible values of juvenile survival rates [w02L , w02U ], with the real w02 values being in between. The survival rates of older groups are calculated using these w02 values. Calculated in this way, the numbers of M7 contain only one outlier (in 1973). The calculations also show the change of population dynamics in the late 1980s, with higher survival rates of bulls. Fig. 2 shows that the survival rate of males between birth and two years of age was changing during the entire period of observation. It significantly decreased by the mid-1960s, remained low during the 1970s, and then oscillated with overall slight growth, which stopped in the early 2000s. (Fig. 2, left). The number of two-year-old males (Fig. 2, right) decreased between 1970 and the 1990s, but not as dramatically as Lander’s method shows. New estimates (Eq. (12)) look more realistic, as this method provides good model approximation of the numbers of adult males (with an average error of 3.2%). The obtained dynamics of the survival rates of different age groups agrees with the understanding of biology of this species. Low survival rates (high death rates) of males during their first two years of life are explained by their high death rate as they adapt to life at sea (up to 70%, Lander, 1975). For the surviving and already adapted animals, the death rate sharply decreases. It then goes up again for the adult animals because of their competition during reproductive age and because of natural aging. Our study shows that the slow growth of the Tyuleniy populations is most likely not caused by the decrease of the survival rates of juvenile or any other age group. Indeed, the juveniles do not exhibit visible decline, and older groups even seem to be slowly growing. However, the study of density-dependent effects on the early survival rate showed a significant increase of intraspecific competition, which is most likely related to some extrinsic factors, such as the quality and availability of food. These negative changes are happening on the seemingly normal background. There is no considerable decline of juvenile survival because the numbers of newborn pups are low. Above all, the most daunting question is: Why does the number of pups not grow while the number of bulls is at an all-time high? Also, the dynamics of pup viability is telling (coefficient a in Eqs. (18)–(20)) – it was the highest in the early period, then declined with intensifying harvest, and then resumed

Fig. 9. Average number of pups per bull (P(n)/M(n-1), on the left axis, and the number of bulls in the rookery (M(n-1)), on the right axis).

growing after strong conservation measures, but never reached the original level. All these facts may be the consequence of negative genetic changes in the population, and raise questions about the genetic consequences of the harvest. We believe that the most important question for future studies is the relationship between the sex ratio in the population and the pregnancy rate. This is the way to answer the question about the possible selectiveness of harvest. However, even now, without information about the physiological status of the females, it is possible to estimate the average reproductive efficiency of one bull (which is a proxy for an average harem size), using the estimated number of the bulls in previous year and the number of newborn pups (Fig. 9). The change in the reproductive efficiency is readily seen. At the start of observations, it was 20 pups per bull and higher, with the sharp peak in 1970s; then it dropped to 10–20; and since the late 1980s (parallel with the steady growth of the number of bulls), it dropped to the next level around 10, which has remains stable for over 20 years. Such dynamics may be a consequence of selective harvest. At first, when the bulls were still “strong,” the harvest allowed them to increase efficiency by decreasing the competition. Later, when the harvest removed the best potential sires, only the inferior bulls remained. Certainly, this is only one possible version, and more studies are necessary (including analysis of the population and physiological parameters of females) in order to confirm or disprove it. Note that this local northern fur seal herd is only one example of many commercial species with a complex age structure and strictly seasonal life cycle. This means that the intra-population changes revealed in this population may be characteristic of many other exploited populations, and these changes appear to be responsible for difficulties with the restoration of the populations after their exploitation is prohibited. Acknowledgement This work was partially supported by the Russian Foundation for Basic Researches [grant number 15-29-02658]. References Abakumov, A.I., Kucher, A.I., 1996. Optimal exploitation of a fish population. Hydrobiol. J. 32 (7), 38–46. Abakumov, A.I., Bocharov, L.N., Reshetnyak, T.M., 2009. The optimal quota distribution in multispecies fishing on an example of Karagandinsky region. Voprosy rybolovstva 10, 352–363 (in Russian).

O.L. Zhdanova et al. / Ecological Modelling 363 (2017) 57–67 Abakumov, A.I., Il’in, O.I., Ivanko, N.S., 2016. Game problems of harvesting in a biological community. Autom. Remote Control 77 (4), 697–707. Andreev, V.L., Bulgakova, T.I., Chelnokov, F.G., 1978. An estimation method of some parameters of fur seal population using mark materials. Trudy VNIRO 128, 23–32, in Russian. Antonelis, G.A., Ragen, T.J., Rooks, N.J., 1994. Male-biased secondary sex ratios of northern fur seals on the Pribilof Islands, Alaska, 1989 and 1992. In: Sinclair, E.H. (Ed.), Fur Seal Investigations. NOAA Tech. Mem., pp. 84–89. Boltnev, A.I., 2011. Nothern Fur Seal of Comandor Islands. Izdatel’stvo VNIRO, Мoskow, in Russian. Buddle, C.M., Langor, D.W., Pohl, G.R., Spence, J.R., 2006. Arthropod responses to harvesting and wildfire: implications for emulation of natural disturbance in forest management. Biol. Conserv. 128 (3), 346–357. Caputo, F.P., Canestrelli, D., Boitani, L., 2005. Conserving the terecay (Podocnemis unifilis, Testudines: Pelomedusidae) through a community-based sustainable harvest of its eggs. Biol. Conserv. 126 (1), 84–92. Chapman, D.G., Robson, D.S., 1960. The analysis of a catch curve. Biometrics 16 (3), 354–368. Chapman, D.G., 1964. Critical study of Pribilov fur seal population estimates. Fish. B-NOAA 63 (3), 657–669. Clark, C.W., 1990. Mathematical Bio-Economics: The Optimal Management of Renewable Resources, second ed. Wiley, New York. Corey, J.A., Bradshaw, R.J., Barker, R.G., Harcourt, A., Lloyd, S.D., 2003. Estimating survival and capture probability of fur seal pups using multistate mark–recapture models. J. Mammal. 84 (1), 65–80. Eberhardt, L.L., 1981. Population dynamics of the Pribilov fur seals. In: Fowler, C.W., Smith, T.D. (Eds.), Dynamics of Large Mammal Populations. J. Wiley&Sons, New York, pp. 197–220. Fauler, Ch., 1998. Northern fur seals of Pribilof Islands, Severnyj morskoj kotik: sistematika, morfologija, ecologija, povedenie. Мoskow, pp. 450–498 (in Russian). Frisman, E. Ya., Skaletskaya, E.I., Kuzyn, A.E., 1982. A mathematical model of the population dynamics of a local northern fur seal with seal herd. Ecol. Model. 16, 151172. Frisman, E. Ya., Skaletskaya, E.I., Kuzin, A.E., 1985. Mathematical Modelling the Number Dynamics of the Northern Fur Seal and Optimal Control of Fur Seal Farm. FESC AS USSR, Vladivostok, in Russian. Frisman, E. Ya., Last, E.V., Skaletskaya, E.I., 2006. Population dynamics of harvested species with complex age structure (for Pacific salmons fish stocks as an example). Ecol. Modell. 198, 463–472. Frisman, E. Ya., Zhdanova, O.L., Kolbina, E.A., 2010. Effect of harvesting on the genetic diversity and dynamic behavior of a limited mendelian population. Russ. J. Genet. 46 (2), 239–248. Graham, M., 1935. Modern theory of exploiting a fishery, and application to North Sea trawling. J. Conseil 10 (3), 264–274. Hauser, C.E., Cooch, E.G., Lebreton, J.-D., 2006. Control of structured populations by harvest. Ecol. Modell. 196, 462–470. Hauser, C.E., Runge, M.C., Cooch, E.G., Johnson, F.A., Harvey, W.F., 2007. Optimal control of Atlantic population Canada geese. Ecol. Modell. 201, 27–36. Hjort, J., Jahn, G., Ottestad, P., 1933. The optimum catch. Hvalrådets Skrifter 7, 92–127. Horan, R.D., Bulte, E.H., 2004. Optimal and open access harvesting and multi-use species in a second best world. Environ. Resour. Econ. 28 (3), 251–272. Ichihara, T., 1972. Maximum sustainable yield from the Robben fur seal herd//Far Seas Fish. Res. Lab. Bull. 6, 77–94.

67

Interim Convention on Conservation of North Pacific Fur Seal Washington, 1957s. Washington, 1957. http://sedac.ciesin.columbia.edu/entri/texts/acrc/1957FS. txt.html (Accessed 17 February 2007). Kornev, S.I., Blokhin, I.A., Generalov, A.A., Semerinov, A.P., 2008. Historical trend of the population of Northern fur seals in the Commander islands for the 50-year-period (1958–2007). Issledovaniya vodnyh biologicheskih resursov Kamchatki i Severo-Zapadnoj chasti Tihogo okeana. 11, 105–120. (in Russian). Kuzin, A.E., 1999. Northern fur seal. Sovet po morskim mlekopitayuschim, Мoskow, in Russian. Kuzin, A.E., 2003. Resources of northern fur seal from Sakhalin region and their usage. Izv. TINRO 133, 138–144, in Russian. Kuzin, A.E., 2010. The intrapopulation structure of the northern fur seal (Callorhinus Ursinus) on Tyuleniy Island during the post-depression years (1993–2009). Rus. J. Mar. Biol. 36 (7), 507–517. Kuzin, A.E., 2014. New data on the abundance of the northern fur seal (Callorhinus Ursinus), steller sea lion (Eumetopias Jubatus), and spotted seal (Phoca Largha) on Tyuleniy Island, Sea of Okhotsk. Russ. J. Mar. Biol. 40 (7), 532–538. Kuzin, A.E., 2015. Analisys of northern fur seal (Callorhinus ursinus) harvest on Tyuleniy Island. Izv. TINRO 183, 71–80, in Russian. Lander, R.H., 1975. Method of determining natural mortality in the northern fur seal (Callorhinus ursinus) from known pups and kill by age and sex. J. Fish. Res. Board Can. 32 (12), 2447–2452. Martin, J., O’Connell, A.F., Kendall, W.L., Runge, M.C., Simons, T.R., Waldstein, A.H., Schulte, S.A., Converse, S.J., Smith, G.W., Pinion, T., Rikard, M., Zipkin, E.F., 2010. Optimal control of native predators. Biol. Conserv. 143 (7), 1751–1758. Neverova, G.P., Abakumov, A.I., Frisman, E.Y., 2016. Dynamic modes of exploited limited population: results of modeling and numerical study. Math. Biol. Bioinf. 11 (1), 1–13, http://dx.doi.org/10.17537/2016.11.1 (in Russian). Punt, A.E., Smith, A.D.M., Smith, D.C., Tuck, G.N., Klaer, N.L., 2014. Selecting relative abundance proxies for BMSY and BMEY. ICES J. Mar. Sci. 71 (3), 469–483, http://dx.doi.org/10.1093/icesjms/fst162. Ricker, W.E., 1958. Maximum sustained yields from fluctuating environments and mixed stocks. J. Fish. Res. Board Can. 15, 991–1006. Russell, E.S., 1931. Some theoretical considerations on the overfishing problem. ICES J. Mar. Sci. 6, 3–20. Skonhoft, A., ·Vestergaard, N., Quaas, M., 2012. Optimal harvest in an age structured model with different fishing selectivity. Environ. Resour. Econ. 51, 525–544. Smith, T.D., Polacheck, T., 1984. The population dynamics of the Alaska fur seal: what do we really know? NWAFC Proc. Rep. 84 (15), 1–122. Tahvonen, O., 2008. Harvesting an age-structured population as biomass: does it work? Nat. Resour. Modell. 21 (4), 525–550. Trites, A.W., Larkin, P.A., 1989. The decline and fall of the Pribilof fur-seal (Callorhinus Ursinus) – a simulation study. Can. J. Fish. Aquat. Sci. 46, 1437–1445. Trites, A.W., 1989. Estimating the juvenile survival rate of male northern fur seals (Callorhinus ursinus). Can. J. Fish. Aquat. Sci. 46, 1428–1436. Verhulst, P.F., 1838. Notice sur la loi que la population poursuit dans son accroissement. Correspondance mathématique et physique 10, 113–121. York, A.E., 1987. Northern fur seal, Callorhinus ursinus, Eastern Pacific population (Pribilof Islands, Alaska, and San Miguel Island, California). In: Croxall, J.P., Gentry, R.L. (Eds.), Status, Biology, and Ecology of Fur Seals. NOAA Tech. Rep, 51. British Antarctic Survey, Cambridge, pp. 133–140.