Экология, 2019, № 2, стр. 149-160

К чему привел промысел северного морского котика: результаты калибровки математических моделей по данным наблюдений (на примере популяции о-ва Тюлений)

Е. Я. Фрисман a, О. Л. Жданова ab*, А. Е. Кузин c

a Институт комплексного анализа региональных проблем ДВО РАН
679016 Биробиджан, ул. Шолом Алейхема, 4, Россия

b Институт автоматики и процессов управления ДВО РАН
690041 Владивосток, ул. Радио, 5, Россия

c Тихоокеанский научно-исследовательский рыбохозяйственный центр
690950 Владивосток, пер. Шевченко, 4, Россия

* E-mail: axanka@iacp.dvo.ru

Поступила в редакцию 18.07.2018
После доработки 28.08.2018
Принята к публикации 16.08.2018

Полный текст (PDF)

Аннотация

Приводится подробный анализ тенденций динамики промысловой популяции северного морского котика (Callorhinus ursinus) о-ва Тюлений. Результаты исследования показывают, что замедление процессов воспроизводства (уменьшение численности щенков) и падение выживаемости молодых самок (от рождения до 3 лет) привели в конце 1980-х годов к снижению общей численности самок и существенному росту доли самок старшего возраста. Рост выживаемости молодых самок в позднем периоде наблюдений (после 1988–1989 гг.) стал причиной изменения тенденции – численность самок постепенно восстанавливается, доля молодых самок возрастает, а старых (более 10 лет) уменьшается. Изменился и возрастной состав самцов: доля молодых уменьшилась, а секачей выросла. Более того, количество секачей продолжает расти и уже превысило ту их численность, которая когда-то обеспечивала процветание популяции. Возникает парадоксальная ситуация: численность самок восстанавливается, число секачей растет, а численность приплода остается весьма небольшой. Наша работа посвящена объяснению этой ситуации и содержит обоснование гипотезы о том, что именно промысел мог привести к описанным кардинальным изменениям в популяциях северного морского котика. В результате промысла из популяции были изъяты наиболее продуктивные производители, что привело к существенной перестройке возрастной структуры популяции и резкому замедлению роста ее численности.

Ключевые слова: Callorhinus ursinus, возрастная структура, промысел, оценка выживаемости, динамика популяции, математическое моделирование

Деятельность человека в той или иной степени затронула большинство природных популяций и вызвала угрозу существования многих из них [например 1–3]. Особое место занимают популяции промысловых видов. Достаточно давно разработана теоретическая стратегия оптимального промысла [4 и др.], которая, сохраняя популяцию, позволяет получать максимально возможный урожай. Эта стратегия продолжает развиваться [57] и применяется на практике [810]. Однако даже в условиях оптимальной эксплуатации численность многих популяций значительно снижается; при этом даже полная отмена промысла зачастую не приводит к восстановлению исходного состояния популяции [11, 12]. Ответ на вопрос, какие негативные внутрипопуляционные изменения провоцирует промысел (даже оптимальный), далеко не очевиден. В данной работе на примере локальной популяции северного морского котика, которая эксплуатировалась и наблюдалась в течение длительного периода времени, мы попробуем дать ответ на этот вопрос.

В середине прошлого века северный морской котик (Callorhinus ursinus) оказался уникальным объектом популяционных исследований. Большая ценность меха этого животного и “удачный” ареал обитания привели к созданию четырехсторонней международной конвенции (между Канадой, Японией, СССР и США) по сохранению котиков северной части Тихого океана (Interim Convention of Conservation of North Pacific Fur Seals, 1957). В рамках конвенции были организованы масштабные исследования биологии котиков с ежегодной оценкой численности половых и возрастных групп на каждом из крупных лежбищ этих животных, а также определялись условия и размер допустимой добычи. За почти 30-летний период существования конвенции был накоплен уникальный набор данных по популяционной динамике вида. Естественно, эти данные служили хорошей основой для оценки основных популяционных параметров и базой для построения различных математических моделей динамики численности.

В частности, в 1982 г. была предложена достаточно подробная модель динамики численности морского котика и обоснованы методики вычисления ее параметров [13]. К настоящему времени накоплен значительный объем дополнительной биологической информации, собранной сотрудниками ТИНРО-центра на протяжении 56 лет наблюдений за популяцией о-ва Тюлений. Это позволяет продолжить проверку целесообразности использования построенной модели, а также методик оценки ее параметров – адаптировать их в связи с произошедшими изменениями в популяционных процессах.

На фоне многолетнего регулируемого промысла (несмотря на его предполагаемую оптимальность) в популяции морского котика о-ва Тюлений появились признаки депрессии. К концу 1980-х годов численность новорожденных щенков уменьшилась вдвое относительно средних значений для периода 1960-х–начала 1970-х годов, а затем практически стабилизировалась на этом низком уровне. Депрессия настигла и другие популяции морского котика Северной Пацифики: командорскую и прибыловскую [1416]. Чтобы изменить эту негативную тенденцию, промысел на о-ве Тюлений к 1990-м годам был значительно ограничен, а после 2008 г. прекращен. В результате резко выросла численность секачей, но ожидаемого восстановления популяции не произошло: хотя количество новорожденных щенков растет, уровень середины 1960-х годов так и не достигнут (рис. 1).

Рис. 1.

Численность секачей и объем промысла (1, 2) и численность новорожденных щенков (3) на о-ве Тюлений.

Для выявления и, возможно, коррекции изменений в популяции северного морского котика о-ва Тюлений, которые привели к ее замедленному росту даже в условиях отсутствия эксплуатации, необходимо оценить и проанализировать динамику выживаемости различных половых и возрастных групп. На следующем этапе полученные оценки параметров жизненного цикла морского котика можно использовать для математического моделирования и построения прогноза динамики популяции, а также корректировки стратегии эксплуатации.

ХАРАКТЕРИСТИКА ДАННЫХ НАБЛЮДЕНИЙ

Надежные визуальные подсчеты котиков всех половозрастных групп практически не осуществимы – плотность залегания животных на лежбище высока, они находятся в постоянном движении, стада их многочисленны и занимают большую и часто труднодоступную территорию. Непосредственный визуальный подсчет животных, дающий достаточно надежную оценку, оказался возможен только для секачей и щенков.

Секачи – это наиболее крупные животные, число их сравнительно невелико, перемещения по территории лежбища в гаремный период минимальны; они легко идентифицируются в гаремах, где располагаются на заметном расстоянии друг от друга. Существуют также методы оценки численности щенков, так как до определенного возраста они все время проводят на лежбище, не спускаясь в море. Оценки численности секачей и щенков на лежбище были доступны за весь период наблюдений (1958–2013 гг.).

Имеются практически точные данные о численности и возрастном составе животных, добытых в ходе берегового промысла. Коммерческий промысел ограничивался добычей двух-, пятилетних самцов “холостяков”, в основном добывались трех-, четырехлетние. Ряд данных о промысле заканчивается 2008 г.

Данные о наблюдаемой численности гаремных самок на лежбище доступны за период 1980–2103 гг. Однако нужно учитывать, что это заниженная оценка численности половозрелых самок, так как они не все время проводят на лежбище, подолгу кормятся в море, а часть трехлеток вообще не выходит на берег. Нижней оценкой числа половозрелых самок является количество новорожденных щенков. Кроме того, в период с 1958 по 1988 г. осуществлялась добыча самок котиков в море в научно-исследовательских целях, поэтому имеются значительные выборки самок котиков, для каждой из которых известны ее возраст и наличие или отсутствие беременности.

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ДИНАМИКИ ЧИСЛЕННОСТИ ЛОКАЛЬНОЙ ПОПУЛЯЦИИ СЕВЕРНОГО МОРСКОГО КОТИКА

Обозначим через P(n) – общее число новорожденных щенков в n-м году; M0(n) и F0(n) – численности новорожденных самцов и самок соответственно; Mi(n) – число самцов холостяков в возрасте i лет (i = 2,…,6); M(n) – число секачей (половозрелых особей 7 лет и старше); Fi(n) – число самок в возрасте i лет (i = 3,…,10); F11(n) – число самок в возрасте 11 лет и старше;$F\left( n \right) = \sum\nolimits_{i = 3}^{11} {{{F}_{i}}\left( n \right)} $ – общее число половозрелых самок. Все численности меняются год от года, т.е. играют роль переменных модели. За временной шаг моделирования естественно принять 1 год, так как за это время численности особей успевают измениться вследствие появления нового поколения щенков, перехода особей из младших возрастов в старшие (с учетом выживания и естественной смертности особей) и промыслового изъятия.

Выпишем соотношения, связывающие численности выделенных компонентов популяции котиков в разных поколениях:

(1)
$P(n) = 2{{F}_{0}}(n) = 2{{M}_{0}}(n) = B\left( n \right)F(n);$
(2)
${{M}_{2}}(n) = {{w}_{{02}}}(n - 2){{M}_{0}}(n - 2);$
(3)
$\begin{gathered} {{M}_{{i + 1}}}(n) = {{w}_{{i,i + 1}}}(n - 1)({{M}_{i}}(n - 1) - {{R}_{i}}(n - 1)), \\ i = \overline {2,5} ; \\ \end{gathered} $
(4)
$\begin{gathered} M(n) = {{w}_{{6,7}}}({{M}_{6}}(n - 1) - {{R}_{6}}(n - 1)) + \\ + \,\,w(n - 1)(M(n - 1) - R(n - 1)); \\ \end{gathered} $
(5)
${{F}_{3}}(n) = {{v}_{{03}}}(n - 3){{F}_{0}}(n - 3);$
(6)
${{F}_{{i + 1}}}(n) = {{v}_{{i,i + 1}}}(n - 1){{F}_{i}}(n - 1),\,\,\,\,i = \overline {3,9} ;$
(7)
$\begin{gathered} {{F}_{{11}}}(n) = {{v}_{{10,11}}}(n - 1){{F}_{{10}}}(n - 1) + \\ + \,\,v(n - 1){{F}_{{11}}}(n - 1). \\ \end{gathered} $

Здесь B(n) – коэффициент рождаемости (или доля рожавших в данном году самок среди всех половозрелых самок); w02(n) – коэффициент выживаемости самцов котиков, родившихся в n-м году, на первых двух годах жизни; v03(n) – коэффициент выживаемости самок котиков, родившихся в n-м году, на первых трех годах жизни; wi, i +1(n) и vi, i +1(n) – коэффициенты выживаемости холостяков и самок, доживших в n-м году до i лет, на i +1-м году жизни (т.е. от i-го до i + 1-го года); w(n) и v(n) – выживаемости секачей и самок 11 лет и старше между n-м и n + 1-м годом; Ri(n) – величина промыслового изъятия из i-й группы холостяков ($i = \overline {2,6} $) в n-м году; R(n) – величина промыслового изъятия секачей в n-м году.

Степень детализации представленной модели обусловлена имеющимися в наличии данными. Поскольку практически нет надежных сведений о возрастной структуре секачей и самок старших возрастов (после 10 лет), а также о численности годовалых самцов (годовики никогда не добывались, визуальные их подсчеты невозможны, поскольку они практически отсутствуют на лежбище) и одно-, двухлетних самок (они почти не попадали в пробы), используются агрегированные переменные: возрастные группы секачей всех возрастов (7 лет и старше) и самок старше 10 лет, а также коэффициенты выживаемости самцов от рождения до двух лет и самок до трехлетнего возраста.

Правомерность уравнения (1) подтверждается многочисленными исследованиями физиологии морского котика: так, половое соотношение эмбрионов практически 1 : 1, и хотя к моменту рождения встречаются отклонения от этого соотношения для территориально различных популяций как в одну, так и в другую сторону [17, 18 и др.], в количественном выражении они незначительны. Для тех периодов, для которых может быть оценена возрастная структура группы половозрелых самок, выражение $B\left( n \right)F(n)$ в уравнении (1) можно заменить эквивалентным, но более подробным: $\sum\nolimits_{i = 3}^{11} {{{B}_{i}}\left( n \right){{F}_{i}}\left( n \right)} ,$ где Bi(n) – коэффициент рождаемости самок i-го возраста (или доля рожавших в данном году среди всех i-летних самок ).

Модель (1‒7) не учитывает миграции котиков из одних популяций в другие, так как исследование миграционных процессов на основе данных о возврате меток показало высокую степень изоляции стада о-ва Тюлений от остальных репродуктивных группировок северных котиков [19].

ОЦЕНКА ВЫЖИВАЕМОСТИ РАЗЛИЧНЫХ ПОЛОВОЗРАСТНЫХ ГРУПП КОТИКОВ

Выживаемость самцов

Оценки ювенильной выживаемости самцов северного морского котика широко использовались в математических моделях для изучения причин снижения численности прибыловских популяций [14, 20, 21], для моделирования динамических эффектов изменения ювенильной выживаемости на рождаемость и численность самцов популяции о-ва Тюлений [13], оценки влияния промысла самок на рождаемость [22]; они также использовались для определения максимального устойчивого объема изъятия [19, 20].

Оценка ювенильной выживаемости (самцов от рождения до 2 лет) w02 [21, 23, 24] основывается на определении верхней (w02U) и нижней (w02L) границы в предположении, что среднегодовая выживаемость холостяков рассматриваемого поколения (самцов от 2 до 7 лет, w27) постоянна: w27(n) = w2, 3(n) = w3, 4(n + 1) = = w4, 5(n+2) = w5, 6(n + 3) = w6, 7(n + 4).Здесь необходимо подчеркнуть, что коэффициент w02 обозначает выживание за два первых года жизни, а коэффициент w27 за один год – каждый год от 2 до 7 лет. Такое обозначение носит уже традиционный характер, и мы вынуждены его придерживаться.

Изменение характера промысла привело к необходимости модификации используемых ранее методик оценки параметра w02 [23, 24]. Исключение предположения об объеме изъятия, а также включение информации о пополнении секачей самцами данного поколения при достижении ими возраста семи лет позволило получить удовлетворительные оценки выживаемости ювенильных особей и холостяков [25, 26]. Методы итераций и перебора [26] дают качественно похожие результаты: очень высокая ювенильная выживаемость в начале периода наблюдений (1958–1963 гг.), резкое падение в 1966 г. и т.д. (рис. 2 а). Второй метод увеличивает выживаемость холостяков и стремится уменьшить ювенильную выживаемость относительно первого метода. Корректность каждого набора оценок была проанализирована методом восстановления динамики по полученным оценкам.

Рис. 2.

Диаграммы размаха выживаемости самцов (a) на первом году жизни w01 (слева) и w27 (справа): 1958–1971, 1972–1977, 1978–1987 и 1988–2006 гг., рассчитанные методами итераций (сверху) и перебора (внизу); б – наблюдаемая и расчетная (1 – по методу итераций, 2 – перебора) динамика численности секачей (периоды 1965–1988 и 1989–2011) без поправки в позднем периоде. Средняя ошибка аппроксимации 4.8% и 5.5% соответственно.

Средняя выживаемость секачей (w) была получена из уравнения (4), которое представляет собой модель парной линейной регрессии. Методом наименьших квадратов (МНК) оценивали уравнение вида

(8)
$M(n) - {{M}_{7}}(n) = w(M(n - 1) - R(n - 1)),$
в котором используются оценки M7(n), полученные на основе рассчитанных ранее ювенильных выживаемостей и выживаемостей холостяков.

Оказалось, что данные неоднородны; причем изменение тенденции, произошедшее где-то в конце 1980-х годов, демонстрируют оценки, полученные обоими методами: тест отношения правдоподобия (likelihood ratio test) дает вероятность изменения тенденции для 1988 г., p < 0.01. Исходя из этого вся выборка была разбита на два периода: 1965–1988 гг. и 1989–2013 гг. Ниже приведены уравнения каждой из регрессий со значениями коэффициентов, полученными для оценок M7(n) на основе методов итераций (А) и перебора (Б):

$\begin{gathered} 1965 - 1988\,\,{\text{г г }}.\,:\,\,({\text{А }})\,\,\,\,M(n) - {{M}_{7}}(n) = 0.75(M(n - 1) - R(n - 1))--124,\,\,\,\,{{R}^{2}} = 0.96, \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,({\text{Б }})\,\,\,\,M(n) - {{M}_{7}}(n) = 0.74(M(n - 1) - R(n - 1))--115,\,\,\,\,{{R}^{2}} = 0.97; \hfill \\ \end{gathered} $
$\begin{gathered} 1989 - 2013\,\,{\text{г г }}{\text{.}}\,:\quad\,\,({\text{А }})\,\,\,\,M(n) - {{M}_{7}}(n) = 0.67(M(n - 1) - R(n - 1)) + 109,\,\,\,\,{{R}^{2}} = 0.94, \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,({\text{Б }})\,\,\,\,M(n) - {{M}_{7}}(n) = 0.65(M(n - 1) - R(n - 1)) + 145,\,\,\,\,{{R}^{2}} = 0.95. \hfill \\ \end{gathered} $

Благодаря разделению всего периода наблюдений на ранний (до 1988 г.) и поздний (после 1988 г.) удалось избавиться от гетероскедастичности – дисперсия ошибки стала однородной и это верно для каждой из полученных регрессий. Коэффициенты детерминации остались достаточно высокими (94–97%).

Отметим основные изменения, которые произошли в конце 1980-х годов. Несколько снизилась средняя выживаемость секачей, при этом смещение оценок численности семилеток M7 потеряло статистическую значимость, т.е. оценка M7 стала несмещенной! В раннем периоде (1965–1988 гг.) количество семилеток в среднем было переоценено на 115–124 особи. На рис. 2б представлена модельная (с фиксированной выживаемостью секачей из приведенных выше регрессионных оценок (8) для данных двух методов) и наблюдаемая динамика численности секачей. Средняя ошибка аппроксимации составляет 4.8 и 5.5% для метода итераций и перебора соответственно.

Выживаемость самок

По методике, предложенной ранее [19], были обработаны имеющиеся данные о возрастном составе и физиологическом состоянии самок морских котиков, добытых в море за 31-летний период (1958–1988 гг.). На рис. 3 представлена динамика долей беременных особей для этого периода (а), оценки выживаемостей самок от рождения до 3 лет отдельно по поколениям (б), среднегодовые показатели выживаемости для каждого возраста с 3 до 9 лет (в), а также представлен возрастной состав самок (г), рассчитанный на основе приведенных оценок их выживаемости. Можно заключить, что в период 1958–1988 гг. происходило “старение” женской части популяции – увеличивалась доля старших возрастов и уменьшалась младших на фоне резкого снижения выживаемости младшей возрастной группы (от рождения до 3 лет).

Рис. 3.

Оценки состояния самок по морским пробам: а – доля беременных самок в морских пробах по годам вместе с 95%-ными доверительными интервалами; б – динамика оценок выживаемости самок до трех лет за весь период морских исследований; в – динамика средней выживаемости самок по возрастам; г – возрастной состав самок, % от общей численности (столбец численности формируется снизу вверх от меньших возрастов к большим – 3-летки, потом 4-летки и т.д.); д – динамика наблюдаемой численности самок на лежбище (1) и модельной (2); коэффициенты модели: v03 = 0.625, v3, 4 = v4, 5 = …v9, 10 ≈ 0.87, v10+ = 0.86.

Полученные коэффициенты были использованы для моделирования динамики численности половозрелых самок вплоть до настоящего времени [27]. Оценка адекватности модельной численности самок производилась с учетом данных о наблюдаемой численности новорожденных щенков на лежбище (1958–2013 гг.). Оказалось, что после 1988 г. модельная динамика абсолютно не согласуется с данными учета щенков на лежбище, что позволило предположить изменение коэффициентов выживаемости самок в более позднем периоде в сторону их увеличения, иначе они не смогли бы произвести наблюдаемое количество новорожденных особей. В ходе численных экспериментов были подобраны новые коэффициенты выживаемости самок, на основе которых получена достаточно адекватная модельная динамика их численности, согласующаяся с реальной динамикой численности щенков и визуально оцененной динамикой численности самок (см. рис. 3д). При этом лучшими оказались результаты при разделении всего периода наблюдения на две части и соответственно при использовании различных коэффициентов выживаемости самок в раннем и позднем периодах. Это подтверждает предположение об изменении выживаемости самок после 1988 г., причем уверенно можно констатировать рост среди них выживаемости младшей возрастной группы – от рождения до 3 лет.

МОДЕЛЬ ДИНАМИКИ ЧИСЛЕННОСТИ И ВОЗРАСТНОЙ СТРУКТУРЫ ПОПУЛЯЦИИ СЕВЕРНОГО МОРСКОГО КОТИКА О-ВА ТЮЛЕНИЙ

Динамика численности и возрастной структуры самок популяции северного морского котика о-ва Тюлений

Далее приводится подробный модельный анализ динамики численности и возрастной структуры самок в популяции. Представим их жизненный цикл в виде графа (рис. 4а). В соответствии с графом жизненного цикла (ГЖЦ) динамику численности и возрастного состава самок описывает модель Лесли–Лефковича:

(9)
$F(n + 1) = {\text{L}}F(n).$
Рис. 4.

Граф жизненного цикла самок(а) и самцов (б) северного морского котика.

Вектор-столбец F(n) = (F0(n), F1(n), F2(n), …, F10+(n))T представлен численностями самок соответствующего возраста в n-м году; проекционная матрица L модели имеет следующий вид:

(10)
${\text{L}} = \left( {\begin{array}{*{20}{c}} 0& \cdots &0&{{{b}_{3}}}& \cdots &{{{b}_{9}}}&{{{b}_{{10}}}}&{{{b}_{{10 + }}}} \\ {{{v}_{{01}}}}&0& \cdots & \cdots & \cdots & \cdots & \cdots &0 \\ 0&{{{v}_{{12}}}}&0& \cdots & \cdots & \cdots & \cdots &0 \\ \vdots &0& \ddots & \ddots & \vdots & \vdots & \vdots & \vdots \\ \vdots & \vdots &0& \ddots & \ddots & \vdots & \vdots & \vdots \\ \vdots & \vdots & \vdots & \ddots & \ddots &0& \vdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \ddots &{{{v}_{{910}}}}&0&0 \\ 0& \cdots & \cdots & \cdots & \cdots &0&{{{v}_{{910 + }}}}&{{{v}_{{910 + }}}} \end{array}} \right).$

Численность новорожденных самок F0(n) равна половине численности всех новорожденных щенков P(n), а репродуктивная группа представлена половозрелыми самками всех возрастов от 3 лет и старше; рождаемость bi – это половина от доли беременных в каждом возрастном классе.

Так как матрица вида (10) неразложима и примитивна, здесь применима классическая теорема Перрона–Фробениуса [28]: в спектре неразложимой матрицы L имеется положительное число λ1 = ρ(L) (доминантное собственное число), которому соответствует положительный собственный вектор х* > 0: Lх* = λ1х*, называемый перроновым вектором, при условии нормировки: сумма его компонент равна 1. В математической популяционной биологии такой вектор называется относительной возрастной структурой популяции.

Асимптотически x(t) в матричном уравнении динамики представляет собой геометрический рост со скоростью λ1, и потому λ1 считается аналогом мальтузианского параметра в моделях роста (неструктурированных) популяций. Собственные числа матрицы L определяются из характеристического уравнения: |L – λI| = 0.

Рассмотрим отдельно ранний период 1958– 1988 гг., для которого все коэффициенты модели доступны из данных морских наблюдений, и поздний – с оценками параметров, полученными в результате численного подбора. Средние доли беременных самок по возрастам в морских пробах составили: B3 = 0.01, B4 = 0.44, B5 = 0.74, B6 = 0.80, B7 = 0.83, B8 = 0.85, B9 = 0.86, B10 = 0.91, B10+ = 0.77. Значения выживаемостей приведены на рис. 3. Для каждого периода численно рассчитывалось максимальное по модулю собственное число проекционной матрицы (λ1) и ее перронов вектор (х*), определяющий естественную стационарную относительную возрастную структуру рассматриваемой популяции (рис. 5а).

Рис. 5.

Стационарная возрастная структура и скорость роста самок (а) и самцов (б) для выделенных периодов.

Отметим, что в раннем периоде скорость роста поголовья самок была недостаточной даже для самовоспроизводства, мальтузианский параметр оказался меньше единицы (λ1 = 0.97), и численность снижалась. Рост выживаемости самок в позднем периоде привел к существенному изменению тенденции – численность самок постепенно восстанавливается (λ1 = 1.03). Заметно изменился и стационарный возрастной состав самок: доля молодых самок выросла, а старых (более 10 лет) уменьшилась.

Если для раннего периода все коэффициенты модели (9) получены по данным морских исследований и результаты расчетов вполне обоснованы, то для позднего периода надежных данных наблюдений (доли беременных, выживаемостей самок) у нас нет и полученный результат – следствие адекватных численных экспериментов. Чтобы дополнить вырисовывающуюся картину динамики, вернемся к наблюдаемым размерам гаремов, численности секачей и новорожденных щенков.

Динамика среднего размера гарема и количества щенков на одного секача (как показатель “репродуктивного успеха” производителя)

Для о-ва Тюлений наблюдения по максимальной численности гаремных самок на лежбище доступны только за период 1980–2013 гг., однако и этих данных достаточно, чтобы визуально увидеть изменение соотношения полов (рис. 6а), произошедшее около 1988 г. Средний размер гарема после 1988 г. уменьшился более чем в 2 раза (рис. 6б, медианы периодов: 11.3 и 5.5 соответственно). При этом значительно снизился межквартильный размах (с 4.9 до 1.6), т.е. средний наблюдаемый размер гарема держится на уровне 5.5 в течение длительного времени и значения в разные годы мало отличаются друг от друга.

Рис. 6.

Cреднее количество самок на одного секача на лежбище (или средний размер гарема) и количество самок на лежбище (а); б – диаграммы размаха среднего размера гарема: 1980–1987 гг. и 1988–2013 гг.; в – среднее количество щенков на одного секача (P(n + 1)/M(n)) и количество секачей на лежбище (M(n)); г – диаграммы размаха среднего “репродуктивного успеха” секача: 1958–1971, 1972–1977, 1978–1987 и 1988–2012 гг.

При этом наблюдаемое количество самок на лежбище может давать значительно заниженную оценку общей численности половозрелых самок, поскольку они подолгу кормятся в море и часть времени проводят вне лежбища. По-видимому, для дальнейшего анализа имеет смысл перейти от оценок размера гарема к оценкам количества щенков, приходящихся на одного секача, или среднего “репродуктивного успеха” секача (рис. 6в).

Очевидно и изменение среднего репродуктивного успеха секача: сначала он достигал 20 щенков на секача и выше (1958–1971 гг.: Me = 19.3, IQR = 2.5), затем после непродолжительного всплеска (1972–1977 гг.: Me = 44.3, IQR = 18.6) произошло некоторое снижение до 10–20 (1978–1987 гг.: Me = 16.6, IQR = 4.9) и с конца 1980-х (вместе с уверенным ростом количества секачей) падение до следующего уровня около 10 (1988–2012 гг.: Me = 7.1, IQR = 2), который держится уже более 20 лет. Такое изменение динамики “плодовитости” секача на о-ве Тюлений может быть результатом селективного промысла. Сначала, пока еще были “сильные” секачи, промысел, резко сократив конкуренцию, позволил им увеличить производительность. Затем, после изъятия основных будущих производителей, остались секачи худшего качества. Аналогичный тренд наблюдается и в командорских популяциях [15, с. 111]: “соотношение секачей к родившим самкам было высоким в 1950-е гг. и составляло более 1 : 40, в 1960-е упало до 1 : 10, а к 1970-м составляло около 1 : 35, далее снижалось до 1 : 8 к началу 1980-х гг. и держится на таком уровне до настоящего времени”.

Динамика численности и возрастной структуры самцов популяции северного морского котика о-ва Тюлений

Имея оценки среднего репродуктивного успеха секача и не рассматривая явно динамику самок, мы можем провести подробный модельный анализ динамики численности и возрастной структуры самцовой части популяции. В соответствии с ГЖЦ (см. рис. 4 б) динамику численности и возрастного состава рассматриваемой части популяции описывает модель Лесли–Лефковича:

(11)

Вектор-столбец M(n) = (M0(n), M1(n), M2(n), M3(n), M4(n), M5(n) M6(n), M7+(n))T представлен численностями самцов соответствующего возраста в n-м году; проекционная матрица А модели имеет следующий вид:

(12)
${\text{A}} = \left( {\begin{array}{*{20}{c}} 0& \cdots & \cdots &0&b \\ {{{w}_{{01}}}}&0& \cdots & \cdots &0 \\ 0&{{{w}_{{27}}}}& \ddots & \vdots & \vdots \\ \vdots &0& \ddots &0&0 \\ 0& \cdots &0&{{{w}_{{27}}}}&w \end{array}} \right).$

Численность новорожденных самцов M0(n) равна половине численности всех новорожденных щенков P(n), а репродуктивная группа секачей (M7+) представлена половозрелыми самцами всех возрастов от 7 лет и старше; рождаемость b – это половина от коэффициента репродуктивного успеха секача.

Чтобы оценить скорость популяционного роста и стационарную возрастную структуру (для каждого из выделенных периодов наблюдений, параметры модели на рис. 2а и 6г), были вычислены максимальное по модулю собственное число проекционной матрицы (λ1) и перронов вектор (х*). Обе оценки дают похожие результаты (рис. 5б, по методу итераций): в позднем периоде скорость роста популяции значительно замедлилась – от 17–27% в год в раннем периоде до 6–7% в позднем. Значительно изменился и ее стационарный возрастной состав: доля щенков уменьшилась с 43–54% в раннем периоде до 34–37% в позднем, а доля секачей выросла от 3–6% до 10–11% соответственно.

ОБСУЖДЕНИЕ

Полученные результаты фактически достаточно убедительно свидетельствуют о кардинальных изменениях в популяциях северного морского котика: изменения затрагивают возрастную структуру и самок, и самцов. Замедление процессов воспроизводства (уменьшение численности щенков) и падение выживаемости молодых самок (от рождения до 3 лет) привели к тому, что к концу 1980-х годов общая численность самок стала снижаться (λ1 = 0.97) и в группе самок существенно выросла доля самок старшего возраста.

Рост выживаемости молодых самок в позднем периоде наблюдений (после 1988–1989 гг.) привел к существенному изменению тенденции – численность самок постепенно восстанавливается (λ1 = 1.03) и заметно изменяется стационарный возрастной состав самок: доля молодых возрастает, а старых (более 10 лет) уменьшается.

Заметим, что мы не связываем напрямую снижение выживаемости самок до 3 лет (выявляемое по морским пробам в раннем периоде наблюдений) со старением репродуктивной группы самок, а рост выживаемости молоди самок позднего периода с ее омолаживанием, поскольку известно, что низкая выживаемость приплода свойственна именно молодым самкам. Выживаемость щенков в первые годы жизни зависит от многих факторов; помимо физиологического состояния новорожденных щенков (на что может влиять возраст матери), на нее могли влиять и такие факторы, как доступность и качество пищи, случайные воздействия, плотностная зависимость и т.д. Трудно предположить, что именно физиологическая неполноценность новорожденных щенков стала причиной низкой ювенильной выживаемости самок того времени, поскольку в этом случае аналогичное снижение должно было быть и у ювенильной выживаемости самцов того же периода, однако у самцов этого не наблюдается.

Возрастной состав репродуктивного ядра самок в данной работе детально не анализировался, хотя это было бы, безусловно, интересно и полезно, но для такого исследования нужны более детальные данные о физиологии самок. Из рис. 5а видно, что действительно увеличилась доля самок дорепродуктивного возраста и заметно уменьшилась доля старых самок (возраста 10+). Именно об этом идет речь в статье. Если оценивать физиологическое состояние новорожденных щенков, то, наверное, его в большей степени могла бы охарактеризовать ранняя смертность (или смертность щенков на лежбище). На рис. 7 приведена динамика доли щенков, павших на лежбище по годам. Явных изменений ранняя смертность щенков не демонстрирует (пик в 1965 г. связан с сильным штормом, из-за которого погибла почти половина щенков). Поэтому у нас нет оснований утверждать, что увеличение доли молодых самок в популяции привело к вступлению в размножение физиологически незрелых самок. Наоборот, наши оценки показывают, что среднегодовая выживаемость самок возросла, и их численность восстанавливается.

Рис. 7.

Смертность щенков на лежбище (1) и численность секачей (2).

Вместе с тем уменьшение численности новорожденных самцов полностью не компенсируется: их среднегодовая выживаемость если и увеличивается, то не очень существенно. В итоге в позднем периоде скорость роста самцовой части популяции замедлилась, значительно изменился и ее стационарный возрастной состав: доля щенков уменьшилась, а доля секачей выросла от 3–6% в раннем периоде до 10–11% в позднем. Более того, количество секачей продолжает расти и давно превысило ту их численность, которая когда-то обеспечивала процветание популяции. Возникает парадоксальная ситуация: численность самок восстанавливается, число секачей растет, а численность приплода остается весьма небольшой.

Проведенное нами исследование фактически объясняет эту ситуацию и дает обоснование нашей гипотезы о том, что именно промысел мог привести к описанным кардинальным изменениям в популяциях северного морского котика. Хотя промысел не должен был быть селективным, результаты расчетов показывают, что он все-таки таковым оказался: невольно под промысел чаще попадали активные самцы с большей тягой к лежбищу и соответственно более успешные будущие производители. Средний репродуктивный потенциал самцов снизился почти катастрофически. Это привело к существенной перестройке возрастной структуры популяции и резкому замедлению роста популяции даже на фоне небольшого увеличения естественной выживаемости практически всех возрастных групп. Наблюдаемые явления отражают экологические закономерности эволюции, описанные Станиславом Семеновичем Шварцем [29]: “Колебания качества популяции в зависимости от численности и возрастной структуры – явления взаимосвязанные. Динамика качества столь же характерное их свойство, как и динамика численности”.

Дальнейший рост популяции северного морского котика, если и произойдет, то будет весьма небыстрым. По-видимому, для такого роста (и особенно для его ускорения) необходима новая перестройка качественного состава популяции: увеличение доли самцов с высоким репродуктивным успехом. Такая перестройка наверняка уже идет в популяции под действием естественного отбора после прекращения промысла. Но, поскольку сейчас доля самцов-производителей с высоким репродуктивным успехом невелика (если судить по текущим низким средним значениям), позитивные процессы естественного отбора займут много времени.

Отметим, что рассмотренная здесь локальная популяция северного морского котика лишь один из примеров многочисленных эксплуатируемых видов с выраженной сезонностью жизненного цикла и сложной возрастной структурой. Следовательно, выявленные нами изменения могут быть общими для многих эксплуатируемых популяций и представляют собой возможное объяснение трудностей восстановления многих популяций после промысла.

Работа выполнена при частичной поддержке РФФИ (проект № 18-04-00073). Авторы подтверждают отсутствие у них конфликта интересов, а также соблюдение этических норм при работе с животными.

Список литературы

  1. Caputo F.P., Canestrelli D., Boitani L. Conserving the terecay (Podocnemis unifilis, Testudines: Pelomedusidae) through a community-based sustainable harvest of its eggs // Biol. Conserv. 2005. V. 126. № 1. P. 84–92.

  2. Buddle C.M., Langor D.W., Pohl G.R., Spence J.R. Arthropod responses to harvesting and wildfire: Implications for emulation of natural disturbance in forest management // Biol. Conserv. 2006. V. 128. № 3. P. 346–357.

  3. Humbert J.-Y., Ghazoul J., Richner N., Walter T. Uncut grass refuges mitigate the impact of mechanical meadow harvesting on orthopterans // Biological Conservation. 2012. V. 152. P. 96–101.

  4. Russell E.S. Some theoretical considerations on the overfishing problem // ICES J. Mar. Sci. 1931. V. 6. P. 3–20.

  5. Фрисман Е.Я., Жданова О.Л., Колбина Е.А. Влияние промысла на генетическое разнообразие и характер динамического поведения менделевской лимитированной популяции // Генетика. 2010. Т. 46. № 2. С. 272–281. [Frisman E.Ya., Zhdanova O.L., Kolbina E.A. Effect of harvesting on the genetic diversity and dynamic behavior of a limited mendelian population // Rus. J. Genet. 2010. V. 46. № 2. P. 239–248.]

  6. Skonhoft A., Vestergaard N., Quaas M. Optimal harvest in an age structured model with different fishing selectivity // Environ. Resour. Econ. 2012. V. 51. P. 525–544.

  7. Неверова Г.П., Абакумов А.И., Фрисман Е.Я. Влияние промыслового изъятия на режимы динамики лимитированной популяции: результаты моделирования и численного исследования // Математическая биология и биоинформатика. 2016. Т. 11. № 1. С. 1–13.

  8. Tahvonen O. Harvesting an age-structured population as biomass: does it work? // Nat. Resour. Model. 2008. V. 21. № 4. P. 525–550.

  9. Martin J., O’Connell A.F., Kendall W.L. et al. Optimal control of native predators // Biol. Conserv. 2010. V. 143. № 7. P. 1751–1758.

  10. Abakumov A.I., Il’in O.I., Ivanko N.S. Game problems of harvesting in a biological community // Autom. Remote Control. 2016. V. 77. Is. 4. P. 697–707.

  11. Myers R.A., Hutchings J.A., Barrowman N.J. Why do fish stocks collapse? The example of cod in Atlantic Canada // Ecol. Appl. 1997. V. 7. P. 91–106.

  12. Allendorf F.W., England P.R., Luikart G. et al. Genetic effects of harvest on wild animal populations // Trends Ecol. Evol. 2008. V. 23. № 6. P. 327–337.

  13. Frisman E.Ya., Skaletskaya E.I., Kuzyn A.E. A mathematical model of the population dynamics of a local northern fur seal with seal herd // Ecol. Model. 1982. V. 16. P. 151–172.

  14. Trites A.W., Larkin P.A. The decline and fall of the Pribilof fur-seal (Callorhinus ursinus) – a simulation study // Can. J. Fish. Aquat. Sci. 1989. V. 46. P. 1437–1445.

  15. Корнев С.И., Блохин И.А., Генералов А.А., Семеринов А.П. Исторический тренд командорской популяции северного морского котика за 50 лет (1958–2007 гг.) // Исследования водных биологических ресурсов Камчатки и северо-западной части Тихого океана. 2008. № 11. С.105–120.

  16. Lee O.A., Burkanov V., Neill W.H. Population trends of northern fur seals (Callorhinus ursinus) from a metapopulation perspective // J. Exp. Mar. Biol. Ecol. 2014. V. 451. P. 25–34.

  17. Фаулер Ч. Северные морские котики на островах Прибылова // Северный морской котик: систематика, морфология, экология, поведение. В 2-х ч. М., 1998. С. 450–498.

  18. Болтнев А.И. Северный морской котик Командорских островов. М.: Изд-во ВНИРО, 2011. 264 с.

  19. Фрисман Е.Я., Скалецкая Е.И., Кузин А.Е. Математическое моделирование динамики численности северного морского котика и оптимальное управление котиковым хозяйством. Владивосток: ДВНЦ АН СССР, 1985. 156 с.

  20. Eberhardt L.L. Population dynamics of the Pribilov fur seals // Dynamics of large mammal populations / Eds. Fowler C.W., Smith T.D. N.Y.: J.Wiley&Sons, 1981. P. 197–220.

  21. Smith T.D., Polacheck T. The population dynamics of the Alaska fur seal: what do we really know? // NWAFC Processed Report. 1984. V. 84. № 15. P. 1–122.

  22. York A.E., Hartley J.R. Pup production following harvest of female northern fur seals // Can. J. Fish. Aquat. Sci. 1981. V. 38. №. 1. P. 84–90.

  23. Lander R.H. 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. 1975. V. 32. № 12. P. 2447‒2452.

  24. Trites A.W. Estimating the juvenile survival rate of male northern fur seals (Callorhinus ursinus) // Can. J. Fish. Aquat. Sci. 1989. V. 46. P. 1428–1436.

  25. Жданова О.Л., Кузин А.Е., Фрисман Е.Я. Динамика выживаемости самцов северного морского котика (Callorhinus ursinus) острова Тюлений (Охотское море) по данным многолетних наблюдений // Зоол. журн. 2017. Т. 96. Вып. 6. С. 720–739. [Zhdanova O.L., Kuzin A.E., Frisman E.Y. The survival dynamics of male northern fur seals (Callorhinus ursinus) on Tyulenii Island, Sea of Okhotsk, based on long-term observations // Biol. Bull. 2017. V. 44. № 9. P. 1174–1191. doi. 10.1134/S1062359017090175]

  26. Жданова О.Л., Кузин А.Е., Фрисман Е.Я. Оценка ювенильной выживаемости самцов северного морского котика (Сallorhinus ursinus): математическое моделирование и анализ данных // Математическая биология и биоинформатика. 2018. Т. 13. № 2. С. 360–375. doi: 10.17537/2018.13.360

  27. Жданова О.Л., Кузин А.Е., Фрисман Е.Я. Математическое моделирование динамики выживаемости самок северного морского котика Сallorhinus ursinus (Linnaeus, 1758) стада острова Тюлений // Биология моря. 2017. Т. 43. № 5. С. 310–320. [Zhdanova O.L., Kuzin A.E., Frisman E.Ya. Mathematical Modeling of the variation in the survival of female northern fur seals, Callorhinus ursinus (Linnaeus, 1758), on Tyuleniy Island // Rus. J. Mar. Biol. 2017. V. 43. № 5. P. 348–358.]

  28. Логофет Д.О., Белова И.Н. Неотрицательные матрицы как инструмент моделирования динамики популяций: классические модели и современные обобщения // Фундаментальная и прикладная математика. 2007. Т. 13. № 4. С. 145‒164.

  29. Шварц С.С. Экологические закономерности эволюции. М.: Наука, 1980. 287 с.

Дополнительные материалы отсутствуют.