Моделирование экологической ниши и пригодности местообитаний бурого медведя (Ursus arctos) в подзоне южной тайги с помощью метода факторного анализа GNESFA Текст научной статьи по специальности «Биологические науки»

МОДЕЛИРОВАНИЕ ЭКОЛОГИЧЕСКОЙ НИШИ И ПРИГОДНОСТИ МЕСТООБИТАНИЙ БУРОГО МЕДВЕДЯ (URSUSARCTOS) В ПОДЗОНЕ ЮЖНОЙ ТАЙГИ С ПОМОЩЬЮ МЕТОДА ФАКТОРНОГО АНАЛИЗА GNESFA

С. С. Огурцов

Центрально-Лесной государственный природный биосферный заповедник, Россия

e-mail: [email protected]

Поступила: 05.04.2020. Исправлена: 19.06.2020. Принята к опубликованию: 12.07.2020.

В статье представлены результаты моделирования экологической ниши и пригодности местообитаний бурого медведя (Ursus arctos), выполненные для территории Центрально-Лесного заповедника и его охранной зоны (Тверская область, Россия) с помощью метода факторного анализа GNESFA. Вкратце рассмотрены основные правила выполнения подобных работ, подходы и техники моделирования. Моделирование экологической ниши проводили с помощью анализов FANTER, ENFA и MADIFA, а моделирование пригодности местообитаний - с помощью анализа по расстояниям Махаланобиса. В качестве эколого-географических переменных использовали вегетационные индексы, морфометрические характеристики рельефа, растры близости и типы ландшафтного покрова. Точки встреч вида регистрировали на постоянных маршрутах с помощью GPS-навигатора в течение периода 2008-2018 гг. В качестве итоговых выбраны 6 сценариев с различным сочетанием входных данных (точек встреч и параметров среды). В качестве набора независимых тестовых данных использовали точки встреч бурого медведя, зарегистрированные сотрудниками отделов охраны и науки. Отдельно смоделированы сценарии влияния весовых коэффициентов на точки встреч и мультиколлинеарности на переменные. Маргинальность лучшей модели составила 1.17, а толерантность - 0.69. Расстояние до антропогенных источников пищи, индекс NDVI и проективное покрытие лугов внесли наибольший вклад по оси маргинальности. По первым трем осям специализации наибольший вклад внесли проективное покрытие хвойных лесов, сфагновых сосняков и расстояние до рек. Значения AUC варьировали от 0.66 до 0.73. Значения непрерывного индекса Бойса - от 0.85 до 0.93. Территория частично нарушенной охранной зоны оказалась более пригодной для бурых медведей, чем малонарушенная территория Центрально-Лесного заповедника. Весовые коэффициенты и число эколого-географических переменных оказали существенное влияние на результаты всех анализов.

Ключевые слова: Carnivora, ENFA, HSM, MADIFA, Ursidae, биотоп, ГИС, пространственное моделирование, расстояния Махаланобиса

Введение

Оценка отношений между живыми организмами в популяциях и их окружающей средой требуется в большинстве экологических исследований, как с точки зрения теории, так и с точки зрения практики. Среди доступных решений этого вопроса особенно интересно моделирование пригодности местообитаний (Habitat Suitability Modelling (HSM): Guisan & Zimmermann, 2000; Manly et al., 2002; Brotons et al., 2004; Elith et al., 2006; Peterson et al., 2011; Guisan el al., 2017). Большинство методов HSM основаны на хорошо известной концепции экологической ниши Хатчинсона (Guisan & Zimmermann, 2000; Hirzel & Le Lay, 2008). Эта концепция (Hutchinson, 1957) предоставляет удобный способ анализа данных о присутствии вида наряду с информацией об изменении окружающей среды («фона»; далее presence/background). Согласно ей, экологическая ниша вида определяется как n-мерный гиперобъем, в котором каждая точка соответствует состоянию окружающей среды, позво-

ляющей организмам поддерживать жизнеспособные популяции. Тогда каждая переменная среды соответствует измерению в так называемом экологическом пространстве, которое определяет доступную среду обитания для изучаемого вида. Наблюдаемые присутствия применяются для оценки использования пространства, то есть одного из аспектов экологической ниши (Basille et al., 2008; Calenge & Basille, 2008). Далее под понятием «ниша» мы будем подразумевать именно проекцию пространственного аспекта экологической ниши, не касаясь трофического и временного.

Несколько статистических анализов были разработаны для оценки пригодности местообитаний с использованием presence/background данных (Elith et al., 2006; Peterson et al., 2011). Наиболее популярными среди них являются Maximum Entropy (MaxEnt: Phillips et al., 2006) и Ecological Niche Factor Analysis (ENFA: Hirzel et al., 2002). Согласно структуре Austin (2007) и Peterson et al. (2011), ENFA относится к теории экологической ниши, опираясь на presence/

background данные, и занимает свое место в хорошо изученном семействе многомерных анализов. В нем ENFA принадлежит к факторным анализам, которые обладают многочисленными положительными качествами для исследования данных в многомерных пространствах, особенно для визуализации их структур (например, Thioulouse & Chessel, 1992; Dray et al., 2003). По этой причине факторный анализ может играть ключевую роль в исследованиях взаимосвязей вида с местообитаниями. Он может быть использован для сортировки факторных осей, которые объясняют большую часть различий между нишей вида и его средой. Паттерны системы «ниша-среда», выявленные в результате анализа, могут затем быть связаны с предпочтениями (пригодностью местообитаний) или требованиями (характеристиками ниши) вида в отношении его среды обитания (Calenge & Basille, 2008).

Факторный анализ экологической ниши (ENFA) и факторный анализ по расстояниям Махаланобиса (Mahalanobis Distance Factor Analysis (MADIFA): Calenge et al., 2008) являются двумя такими методами. С одной стороны, ENFA различает два вида информации, измеряемой в системе «ниша-среда»: маргинальность и специализацию. Маргинальность является мерой своеобразия ниши относительно распределения доступных ячеек в экологическом пространстве, тогда как специализация является мерой ограничения ниши относительно распределения доступных ячеек. Таким образом, ENFA предоставляет возможность измерения реализованной ниши в доступном пространстве на основе вычисления двух указанных параметров. С другой стороны, MADIFA опирается только на одну характеристику - расстояние Махаланобиса (D2), которое может быть интерпретировано как индекс пригодности местообитаний (Habitat Suitability Index (HSI): Calenge & Basille, 2008; Calenge et al., 2008).

С формальной точки зрения, ENFA и MADIFA тесно связаны (Calenge et al., 2008). Оба этих метода можно обобщить единым подходом, который получил название факторный анализ системы «ниша-среда» (General Niche-Environment System Factor Analysis, GNESFA). Эта структура также включает в себя еще один анализ - факторный анализ ниши со средой в качестве референтного распределения (Factor Analysis of the Niche, Taking the Environment as the Reference, FANTER), который предла-

гает третью точку зрения относительно этой системы. Эти три исследовательских метода органично дополняют друг друга. При одновременном использовании они представляют обширную характеристику экологической ниши вида (Calenge & Basille, 2008).

ENFA успешно применялась в экологических исследованиях в широком диапазоне таксонов организмов. Целый ряд работ был выполнен по растениям (Zaniewski et al., 2002; Valle et al., 2011; Silva et al., 2017), насекомым (Gallego et al., 2004; Chefaoui et al., 2005), водным беспозвоночным (Galparsoro et al., 2009; Hemery et al., 2011), амфибиям (Préau et al., 2018), рептилиям (Santos et al., 2006; Halstead et al., 2010; Carretero & Sillero, 2016), птицам (Brotons et al., 2004; Hirzel et al., 2004a; Braunisch et al., 2008; Strubbe & Matthysen, 2008; Zimaroeva et al., 2015), млекопитающим (Hirzel et al., 2002; Dettki et al., 2003; Reutter et al., 2003; Mestre et al., 2007; Sattler et al., 2007; Traill & Bigalke, 2007). Отдельно стоит отметить работы по крупным хищникам (Zimmermann, 2004; Basille et al., 2008; Pettorelli et al., 2009; De Angelo et al., 2011; Caruso et al., 2015), в частности по различным видам медведей (Unger et al., 2008; Xuezhi et al., 2008; Qi et al., 2009; Nazeri et al., 2014), в том числе по бурому медведю (Ursus arctos Linnaeus, 1758) (Zajec et al., 2005; Mertzanis et al., 2008; Falcucci et al., 2009; Martin et al., 2012; Nawaz et al., 2014). В последнее время ENFA был адаптирован для оценки чувствительности ниш к изменениям климата (CENFA: Rinnan & Lawler, 2019). До возникновения и начала активного использования MaxEnt именно ENFA долгое время был наиболее популярным методом моделирования для presence-only данных, о чем говорит большое разнообразие выполненных исследований. Работ с использованием MADIFA встречается значительно меньше (например, Thiebot et al., 2011), однако этот анализ успешно сочетается вместе с ENFA (Halstead et al., 2010; Hemery et al., 2011; Martin et al., 2012) или в рамках GNESFA (Calenge & Basille, 2008; Milanesi, 2014; Caruso et al., 2015; Milanesi et al., 2015; Neupane et al., 2019).

Многие авторы единодушно признают, что для более уверенных выводов о предпочтении видом местообитаний необходимо использовать разные подходы на одном и том же наборе данных (Peterson et al., 2011; Guisan et al., 2017; Préau et al., 2018). В рамках наших

исследований бурого медведя с помощью методов экологического моделирования (Огурцов, 2016, 2019), применение GNESFA является следующим этапом многостороннего и комплексного изучения экологической ниши вида. Данная работа является логическим продолжением нашей первой статьи, посвященной моделированию в MaxEnt (Огурцов, 2019), где уже были изложены основные этапы подготовки данных, общие для HSM-исследований. Насколько нам известно, это первая работа в России по моделированию экологической ниши и пригодности местообитаний, выполненная с помощью GNESFA. Именно поэтому большое внимание нами было уделено теоретическим предпосылкам моделирования с помощью данного метода. На русском языке подобные работы опубликованы только для Украины (Демидов и др., 2013; Андрющенко, Жуков, 2016). Поэтому мы будем следовать русскоязычной терминологии, принятой в указанных публикациях. Подробно о принципах ENFA, MADIFA и GNESFA изложено в работах Hirzel et al. (2002), Basille et al. (2008), Calenge & Basille (2008), Calenge et al. (2008), а на русском языке - у Демидова и др. (2013).

Нами были поставлены следующие задачи: 1) продемонстрировать основные возможности моделирования экологической ниши и пригодности местообитаний при помощи GNESFA; 2) построить проекции экологической ниши бурого медведя и выявить ведущие переменные среды, формирующие пространство ниши; 3) построить модели пригодности местообитаний бурого медведя; 4) выявить влияние весовых коэффициентов точек встреч и мультиколлинеарности переменных на качество моделей.

Материал и методы

Территория исследований

Материал собирали на территории Центрально-Лесного государственного природного биосферного заповедника и его охранной зоны (56.4333°-56.5167° N, 32.4833°-33.0167° E), которая представляет собой южнотаежный ландшафт, в разной степени трансформированный человеческой деятельностью. Центрально-Лесной заповедник расположен на юго-западе Валдайской возвышенности (Тверская область, Россия). Это грядово-холмистая возвышенная равнина с абсолютными высотами от 200 м

до 240 м. Его площадь составляет 244.15 км2, а площадь охранной зоны вокруг него - 460.61 м2. Здесь представлено три основных типа местностей: моренные гряды, сформированные рисским и вюрмским оледенением, плоские озерные равнины на месте приледниковых озер и флювиогляциальные равнины в южной части территории исследования (Пузаченко, Козлов, 2007). В схеме ботанико-географиче-ского районирования Центрально-Лесной заповедник размещается в самой южной части Валдайско-Онежской подпровинции в полосе южнотаежных лесов (Исаченко, 1980). Для дренированных частей моренных гряд характерны южнотаежные широколиственно-еловые леса и их производные, для озерных и флювиогля-циальных равнин - бореальные еловые (Picea abies L.) леса. На пологих склонах моренных гряд и в отдельных понижениях озерных равнин развиты верховые болота (Желтухин и др., 2016). По болотам и их окраинам распространены сообщества с участием Pinus sylvestris L., а местами встречаются участки сомкнутых сосняков возрастом 200-300 лет (рис. 1).

В охранной зоне во второй половине XX в. произведены промышленные рубки, уничтожившие почти все еловые леса. В конце XX в. и первом десятилетии XXI в. ранее уцелевшие леса продолжали вырубаться. Небольшие по площади сплошные вырубки были также проведены в южной части Центрально-Лесного заповедника в период его закрытия (1951-1960 гг.). Рубки конца 1950-х гг. и ветровалы обуславливают в составе лесов территории исследования большую долю мелколиственных лесов в разной стадии восстановительной сукцессии. Естественными факторами динамики растительного покрова являются ветровалы, охватывающие большие площади (Желтухин и др., 2016). До 1980-х гг. большая часть охранной зоны Центрально-Лесного заповедника была заселена. Здесь функционировали колхозы, и велась различная сельскохозяйственная деятельность. На сегодняшний день все деревни, кроме двух, заброшены, а общая численность населения не превышает 120 человек. На местах бывшей хозяйственной деятельности в настоящее время развиваются луговые сообщества, распространение которых связано с моренными грядами. Многие из них находятся на различных стадиях зарастания лесной растительностью (Cherednichenko et al., 2016).

Рис. 1. Территория Центрально-Лесного государственного природного биосферного заповедника и его охранной зоны. Классификация типов ландшафтного покрова, выполненная на основе дешифрирования сцены спутника Landsat 8, наложена на трехмерную цифровую модель рельефа Aster GDEM в разрешении 30 м.

Fig. 1. Map of the Central Forest Nature Reserve and its buffer zone. The classification of landscape cover types, based on the Landsat 8 satellite image interpretation, has been superimposed on a three-dimensional Aster GDEM at a resolution of 30 m.

Переменные окружающей среды

Более подробно о выборе переменных среды и их подготовке было оговорено ранее Огурцо-вым (2019). В настоящей работе мы использовали тот же самый набор переменных, которые принято называть в рамках GNESFA эколого-географи-ческими переменными (Ecogeographical Variables (EGV): Hirzel et al., 2002). В этот набор вошли вегетационные индексы (NDVI, GNDVI, NDWI, GRVI, ARVI), морфометрические характеристики рельефа (высота н.у.м., уклон и экспозиция склонов, индекс пересеченности), топографические характеристики местности (расстояния до рек и антропогенных источников пищи (яблоневые сады и овсяные поля)), типы ландшафтного покрова (темнохвойные леса (ельники бореальные, неморальные и папортниковые), сосняки сфагновые, лиственные леса (спелые и молодые), молод-няки и зарастающие вырубки, верховые болота, луга (пойменные и материковые), сельскохозяйственные земли и свежие вырубки). Всего мы использовали 22 переменные (Электронное прило-

жение, табл. 1П). Они отражают орографические, экологические, ландшафтные и антропогенные свойства среды, которые неоднократно были признаны экологически значимыми для моделирования ниши бурого медведя в ENFA (Mertzanis et al., 2008; Falcucci et al., 2009; Martin et al., 2012).

Все переменные среды были приведены к единому географическому экстенту и системе координат WGS84 в проекции UTM 36N с разрешением 30 м. Такой уровень был выбран по минимальному разрешению данных (снимок Landsat и продукт Aster GDEM). В ArcMap был создан грид (полигональная сеть с ячейками 30 х 30 м), охватывающий всю территорию исследования. Несмотря на критику подобного крупномасштабного уровня, он достаточно распространен в HSM-исследованиях (Chefaoui et al., 2005; Xuezhi et al., 2008; Pettorelli et al., 2009; Qi et al., 2009; Nazeri et al., 2014; Neupane et al., 2019). Ячейку грида, для которой измеряются значения всех переменных окружающей среды, в GNESFA принято называть ресурсной единицей («Resource

unit» (RU) или «Environment unit» (EU): Calenge & Basille, 2008; Calenge, 2011).

Переменные, тип распределения которых отличался от нормального, были трансформированы с помощью преобразования Йео-Джонсона (аналога преобразования Бокса-Кокса) и квадратичного преобразования, чтобы их распределение имело симметричный унимодальный вид (Электронное приложение, рис. 1П). Хотя считается, что ENFA достаточно устойчив к отклонению от нормальности (Hemery et al., 2011), такие трансформации повышают качество моделей (Basille et al., 2008; Halstead et al., 2010).

Одной из ключевых особенностей ENFA является невосприимчивость к мультиколлинеар-ности (Hirzel et al., 2002; Dormann et al., 2013). Таким образом, во-первых, снимается требование к независимости переменных. Во-вторых, совершенно не обязательно располагать предварительными сведениями о связях вида со средой обитания и важности для него определенных переменных (Hirzel et al., 2002; Mestre et al., 2007; Nawaz et al., 2014). Для проверки гипотезы о влиянии мультиколлинеарности и числа переменных на результаты моделирования нами использовано три различных набора EGV. Два первых набора были такими же, как и в предыдущем исследовании (Огурцов, 2019), и включали в себя 14 EGV и 18 EGV соответственно.

Третий набор состоял из всех 22 переменных без учета мультиколлинеарности. Подобный подход для выяснения уровня значимости переменных также встречается в моделировании в ENFA (Strubbe & Matthysen, 2008).

Поскольку в нашем случае уже известны предварительные связи вида с местообитаниями, основные сценарии мы рассматриваем на наборе из 14 переменных, прошедших все тесты на мультиколлинеарность (табл. 1). Это также упрощает дальнейшую биологическую интерпретацию результатов (Galparsoro et al., 2009).

Данные о присутствии вида

Как и в случае с MaxEnt, во многих исследованиях крупных хищников с применением GNESFA в качестве точек присутствия вида выбирают различные следы жизнедеятельности (отпечатки лап, экскременты, поеди, визуальные встречи и т.п. - Zajec et al., 2005; Xuezhi et al., 2008; Qi et al., 2009; Martin et al., 2012), только визуальные встречи (Traill & Bigalke, 2007; Pettorelli et al., 2009), регистрации с фотоловушек (Nazeri et al., 2014; Caruso et al., 2015), сочетание следов жизнедеятельности с данными телеметрии (Mertzanis et al., 2008; Falcucci et al., 2009) или только данные телеметрии (Unger et al., 2008). Некоторые работы основываются на неинвазивных генетических пробах (Milanesi et al., 2015).

Таблица 1. Переменные окружающей среды, используемые в качестве EGV в моделировании экологической ниши и пригодности местообитаний бурого медведя (Ursus arctos) в условиях Центрально-Лесного заповедника и его охранной зоны за период 2008-2018 гг.

Table 1. Environmental variables used as EGV in ecological niche and habitat suitability modelling of brown bear (Ursus arctos) in the Central Forest State Nature Reserve and its buffer zone during 2008-2018

№ Название эколого-географической переменной (EGV) Обозначение переменной в моделях Источник данных

Вегетационные индексы

1 NDVI ndvi Landsat 8 OLI

2 NDWI ndwi Landsat 8 OLI

Морфометрические характеристики рельефа

3 высота н.у.м. elevation Aster GDEM

4 экспозиция склонов aspect Aster GDEM

5 уклон (крутизна) склонов slope Aster GDEM

Топографические характеристики местности

6 расстояние до рек river топокарта

7 расстояние до яблоневых садов и овсяных полей artifict съемка с БПЛА

Проективное покрытие типов ландшафтного покрова

8 темнохвойные леса coniferous Landsat 8 OLI *

9 сосняки сфагновые pine bogs Landsat 8 OLI *

10 лиственные леса deciduous Landsat 8 OLI *

11 молодняки и зарастающие вырубки clearcuts Landsat 8 OLI *

12 верховые болота bogs Landsat 8 OLI *

13 луга all fields Landsat 8 OLI *

14 с/х земли, ветровалы, гари и свежие вырубки agriculture Landsat 8 OLI *

Примечание: * - на основе полуавтоматической классификации.

90

В качестве данных о присутствии вида были использованы материалы регистрации точек встреч в ходе регулярных маршрутов по территории исследования. Любые признаки присутствия бурого медведя (отпечатки лап, следы питания, экскременты, лежки, визуальные встречи и т.п.) отмечались с помощью GPS-навигатора Garmin (моделей GPS 60, GPSmap 60CSx, GPSmap 62s, GPSmap 64st) и в дальнейшем рассматривались нами как точки достоверного присутствия вида на территории исследования. Более подробно подготовка точек встреч обсуждена у Огурцова (2019). Поскольку результаты GNESFA следует рассматривать, в том числе, как пригодность местообитаний, был выбран подход, который отражает именно это (HSM-подход, как подробнее показано у Огурцова (2019)). Он включает только те точки, которые имели отношение к использованию видом местообитания. Это любые проявления пищевой деятельности (поеди, жировочные тропы и т.п.), лежки, посещение водопоев и т.п.

В течение 2008-2018 гг. автором собрано 2804 точки присутствия бурого медведя на территории исследования. Для независимых проверок (тестирования) качества моделей из архива Центрально-Лесного заповедника были взяты 1811 точек присутствия бурого медведя, собранных сотрудниками отделов охраны и науки за тот же период времени отдельно от автора.

После выбора из первоначальной («сырой») выборки только тех регистраций, которые несут информацию об использовании видом местообитания осталось 1144 авторских и 417 архивных (тестовых) точек. С ними были проведены все те же процедуры фильтрации (пространственной и временной), в результате чего осталось 566 авторских и 341 тестовых точек. На этапе пространственной фильтрации для каждой ячейки грида оставили только одну точку (центроид), удалив все остальные. Временная фильтрация включала в себя выравнивание сезонных неод-нородностей при сборе данных (подробнее см. Огурцов, 2019). Для проверки гипотезы о влиянии весовых коэффициентов все точки были разбиты на два набора. Первый набор включал в себя все «сырые точки местообитаний» (1144 точки). Второй набор включал в себя эти же точки, но после процедур фильтрации (566 точек). Поскольку одной из задач являлась проверка влияния весовых коэффициентов на качество моделей, мы не использовали тест Морана для снижения пространственной автокорреляции точек (Milanesi et al., 2015).

Весовые коэффициенты

Отличительной особенностью GNESFA является возможность работы с весами, что может быть весьма полезно (Calenge & Basille, 2008), но при этом она не всегда используется в исследованиях (например, Hal stead et al., 2010). Так, для всего грида был создан дополнительный вектор весовых коэффициентов использования (utilisation weight), который отражает число «сырых» точек в каждой ячейке. Все ячейки грида были признаны одинаково доступными для бурых медведей, т.к. даже те из них, которые содержат водные объекты, сельскохозяйственные земли и деревни, включают участки, которые могут быть ими посещены. В связи с этим весовые коэффициенты доступности (availability weight) для всех ячеек были одинаковыми. Суммы коэффициентов весов использования и доступности равнялись 1 (Calenge & Basille, 2008). Совокупности весов использования составляют распределение использованных ресурсов, тогда как совокупности весов доступности - распределение доступных ресурсов.

Желательным требованием к применимости ENFA является предположение о том, что распределения весов использования и весов доступности являются симметричными и унимодальными. В случае MADIFA такое требование распространяется только на веса использования (Calenge & Basille, 2008; Calenge et al., 2008). FANTER лоялен к типу распределений и подходит для анализа мультимодальных ниш (Calenge & Basille, 2008).

Построение моделей

Согласно Calenge & Basille (2008), мы использовали моделирование по методу GNESFA, которое состоит из трех отдельных анализов: FANTER, ENFA и MADIFA. Целью GNESFA является идентификация направления в экологическом пространстве, где два распределения весовых коэффициентов (использования и доступности) максимально различны, что может быть интерпретировано как паттерн ниши. Основной принцип этого анализа заключается в выборе одного из этих двух распределений в качестве референтного, а другого - в качестве фокусного. В зависимости от этого выбора, результат будет принципиально различен (Calenge & Basille, 2008; Calenge, 2011).

Совместное использование трех анализов с графическим отображением ниши в рамках среды позволяет проводить комплексное исследование этой системы. Это придает ценность GNESFA

как общей основы статистического исследования экологической ниши. В большинстве случаев для эколога представляют интерес два различных подхода: система с точки зрения вида и система с точки зрения среды. Поэтому три дополняющих друг друга анализа могут использоваться для построения концептуальной модели исследуемой системы «ниша-среда» (Calenge & Basille, 2008; Caruso et al., 2015; Milanesi et al., 2015; Neupane et al., 2019).

В качестве основных сценариев для моделирования были выбраны шесть, которые отражают рабочие гипотезы о взаимоотношениях вида со средой обитания в различных комбинациях входных данных (табл. 2). Рабочие гипотезы представляют собой проверки значимости весовых коэффициентов на наборе сгруппированных переменных среды, прошедшем проверку на мультиколлинеарность (сценарии 1-2), на наборе всех переменных среды, прошедшем проверку на мультиколлинеарность (сценарии 3-4), на наборе всех переменных среды без учета мультиколли-неарности (сценарии 5-6).

Анализ главных компонент

Поскольку переменные среды имеют разную размерность, первоначально они были стандартизованы, что также является обязательным требованием перед выполнением анализа главных компонент (Principal Component Analysis, PCA). На первом этапе был проведен PCA, куда были включены все переменные. Используя метод сломанного стержня Ма-кАртура и график кумулятивной дисперсии, из них были выделены главные компоненты, которые объясняли наибольшую долю дисперсии (Jackson, 1993; Hirzel et al., 2002, 2004a). В дальнейшем эти компоненты использовали в факторных анализах.

Факторный анализ ниши со средой в качестве референтного распределения (FANTER)

Как видно из самого названия, это особый вид анализа GNESFA который принимает распределение весов доступности в качестве референтного. В этом случае в качестве фокусного распределения выступают веса использования. В связи с этим значение может иметь как первая, так и последняя оси GNESFA. Первые оси соответствуют направлениям, в которых распределение использования в целом является наиболее удаленным от центра тяжести распределения доступности. Эти направления часто соответствуют направлениям, в которых маргинальность является самой сильной. Последние оси соответствуют направлениям, где ширина распределения использования является наименьшей по отношению к ширине распределения доступности. Это те оси, на которых используемые ячейки наиболее сконцентрированы вокруг среднего значения доступности, т.е. где самой сильной является специализация (Calenge, 2011).

Идею FANTER удобно представить, если рассматривать проблему с точки зрения специалиста, хорошо знакомого с изучаемой территорией. Если исследователь обладает предварительными данными о структуре переменных среды в интересуемой области (например, известны взаимосвязи между ними), изучение распределения весов доступности в экологическом пространстве может не представлять первостепенного интереса. В таких случаях целью часто является выявление закономерностей самой ниши и того, чем она отличается от области исследования (Са1е^е & ВавШе, 2008). В связи с этим, FANTER целесообразно использовать в качестве предварительного анализа к Е№А, потому что с его помощью можно провести оценку адекватности рассмотрения экологической ниши в компонентах маргинальности и специализации.

Таблица 2. Входные данные шести рассматриваемых сценариев моделирования экологической ниши и пригодности местообитаний бурого медведя (Ursus arctos) с помощью GNESFA в условиях Центрально-Лесного заповедника и его охранной зоны за период 2008-2018 гг.

Table 2. The input data of the six considered scenarios of brown bear (Ursus arctos) ecological niche and habitat suitability modelling with GNESFA in the Central Forest State Nature Reserve and its buffer zone during 2008-2018

№ сценария Название модели* Набор данных Число переменных Число точек встреч

1 FAHR1 M + G + R 14 1144

2 FAHC1 M + G + C 14 566

3 FAHR2 M + R 18 1144

4 FAHC2 M + C 18 566

5 FAHR3 R 22 1144

6 FAHC3 C 22 566

Примечание: *F - factor analysis (факторный анализ), A - author's (авторские данные), H - habitat suitability (пригодность местообитаний), R - raw («сырые» точки), C - centroid (центроиды точек). M - multicollinearity (проверка на мультиколлинеарность); G - group (групповые переменные).

Факторный анализ экологической ниши (ENFA)

По своей сути, анализ ENFA находится посередине между FANTER и MADIFA. Оба распределения весов используются одновременно как референтные и фокусные. Эта симметричная роль достигается в ущерб экологическому пространству, потому что теряется одно из его измерений через ось маргинальности. ENFA может использоваться для выявления различий положения и узости распределений относительно друг друга. Этот анализ также является оптимальным для моделирования экологической ниши относительно окружающей среды (Calenge & Basille, 2008).

ENFA, как и PCA, объединяет все эколого-гео-графические переменные в несколько некоррелированных факторов, сохраняющих большую часть информации (Hirzel et al., 2002). Основное различие этих двух методов заключается в характере набора данных. PCA вычисляется на одном наборе данных в многомерном пространстве среды. В то же время ENFA использует два набора данных: один - глобальный набор, содержащий значения всех ячеек во всей области исследования, а другой - набор, содержащий только значения для ячеек, где присутствует вид (Hirzel et al., 2004a).

Будучи по своей сути факторным анализом, ENFA выделяет несколько не коррелирующих между собой факторов (осей): одну ось маргиналь-ности и несколько осей специализации (Hirzel et al., 2002; Basille et al., 2008). Ось маргинальности (M) отражает разницу между распределением ресурсных единиц, используемых видом, и распределением ресурсных единиц для всей доступной среды. Таким образом, маргинальность измеряет позицию и величину отклонения ниши относительно экологического пространства. Чем выше маргинальность, тем больше ниша отклоняется от «средних условий» доступной среды обитания, т.е. вид имеет тенденцию существовать в маргинальных (пограничных) условиях (Hirzel et al., 2004a).

Последующие факторы называются факторами специализации (S) и сортируются по убыванию количества объясненной дисперсии. Эти факторы описывают, насколько специализирован вид относительно набора переменных в районе исследований (Hirzel et al., 2002). Только некоторые из первых факторов объясняют основную часть всей информации. Специализация колеблется от 1 до бесконечности, и поэтому ее трудно интерпретировать. По этой причине лучше использовать толерантность (обратную специализацию, 1/S), которая измеряет предпочтения

вида относительно доступного диапазона значений EGV и колеблется от 0 до 1 (Hirzel et al., 2002, 2004a). Ожидается, что случайно выбранный набор ячеек будет иметь толерантность 1, т.е. любое значение ниже 1 указывает на некоторую форму специализации (Pettorelli et al., 2009). Низкие значения (близкие к 0) свидетельствуют о том, что вид является специалистом, а высокие (близкие к 1) - генералистом, т.е. указывают на толерантность вида к условиям среды (Santos et al., 2006; Mestre et al., 2007; Falcucci et al., 2009). Чем больше толерантность - тем шире экологическая ниша вида (Nazeri et al., 2014).

Принцип ENFA заключается в сравнении распределений значений EGV между наборами данных о присутствии вида и всей области исследования. Анализ ищет направления в экологическом пространстве так, чтобы были максимальны два параметра: 1) различие между ресурсами среды, используемыми видом, и ресурсами, доступными в исследуемой области (т.е. маргинальность); 2) отношение между дисперсией доступных ресурсов и дисперсией ресурсов, используемых видом (т.е. специализация, согласно Hirzel et al., 2002; Basille et al., 2008; Calenge et al., 2008; Nazeri et al., 2014).

Биплот ENFA (и в целом GNESFA) - полезная визуализация экологической ниши вида. Он проецирует используемые и доступные ресурсные единицы в экологическом пространстве на плоскость, определяемую осью маргинальности и одной осью специализации. Векторы представляют проекции всех EGV на каждую ось факто-риальной плоскости. Для оси маргинальности их направление и длина выражают влияние на положение экологической ниши. Для оси специализации направление стрелки неинформативно; только ее длина имеет биологический смысл. Самые длинные вектора (вдоль маргинальности или специализации) представляют наиболее влиятельные переменные для соответствующей оси (Nawaz et al., 2014).

Мы можем определить вклад каждой переменной среды в оси маргинальности и специализации, рассчитав коэффициенты корреляции Спирмена (R) и построив их матрицу. Коэффициенты матрицы, относящиеся к фактору мар-гинальности, указывают на корреляцию между каждой EGV и фактором, т.е. определяют величину и направление (предпочтение/избегание) для каждой переменной. Чем больше абсолютное значение коэффициента, тем больший вклад эта EGV вносит в маргинальность. Низкое значение (близкое к 0) указывает на то, что вид

имеет тенденцию жить при средних показателях этой переменной на всей территории исследования. Высокое значение (близкое к 1) указывает на тенденцию к жизни при экстремальных показателях, т.е. условия окружающей среды, используемые видом, сильно отличаются от тех, которые «в среднем» встречаются в районе исследований (Préau et al., 2018). Положительный коэффициент означает, что вид «предпочитает» высокие значения этой EGV, в то время как отрицательный - означает, что вид «предпочитает» низкие значения (Basille et al., 2008). Коэффициент специализации рассчитывается как отношение глобальной дисперсии EGV к дисперсии в части, используемой видом. Это мера ширины ниши в пределах доступной среды обитания. Более высокие абсолютные коэффициенты (знак произвольный) указывают на более специфическое распространение вида для этой EGV (Hirzel et al., 2002; Mestre et al., 2007; Basille et al., 2008; Calenge & Basille, 2008). Таким образом, ENFA по своей сути очень похож на анализ главных компонент. Но в отличие от последнего, его оси имеют непосредственный биологический смысл (Hirzel et al., 2002; Chefaoui et al., 2005).

Полезными показателями являются доли объясненной специализации (Explained Specialisation (ExS), идентичной традиционной объясненной дисперсии) и объясненной информации (Explained Information (ExI) - модифицированной версии ExS, в которой учитывается значение мар-гинальности, согласно Hirzel et al., 2002). Их применяют для того, чтобы показать, как расчетные HS-модели объясняют наблюдаемые результаты.

Несмотря на то, что ENFA использовалась в качестве HSM-метода (Zaniewski et al., 2002; Reutter et al., 2003), важно понимать, что его основная цель - исследование ниши как таковой, а не построение карт пригодности местообитаний (Habitat Suitability Index Map, HSIM) и предсказание значений индексов HSI. Главной причиной этого является различный математический смысл осей маргинальности и специализации, которые некорректно объединять в единое значение HSI (Calenge & Basille, 2008; Calenge et al., 2008). Именно поэтому для построения HSIM в GNESFA лучше использовать анализы на основе расстояний Махаланобиса.

Анализы по расстояниям Махаланобиса (MADIFA и MD)

Анализ по расстояниям (дистанциям) Махаланобиса (Mahalanobis distances, MD) уже до-

статочно давно используется в различных HSM-исследованиях для расчета индексов пригодности местообитаний (Knick & Dyer, 1997; Farber & Kadmon, 2003; Browning et al., 2005; Rotenberry et al., 2006; Tsoar et al., 2007; Halstead et al., 2010), в том числе для крупных хищников (Clark et al., 1993; Corsi et al., 1999; Cayuela, 2004).

MADIFA - это логическое продолжение MD, но в виде факторного анализа. Он опирается только на одну меру пригодности местообитания - расстояние Махаланобиса (D2), рассчитанное в экологическом пространстве от центра тяжести распределения вида на все доступные точки. D2 дает показатель пригодности местообитания в каждой точке, поскольку он указывает на ее отклонение от оптимума вида. Предполагается, что низкое значение расстояния указывает на высокую пригодность (Clark et al., 1993; Knick & Dyer, 1997). MADIFA возвращает направления в экологическом пространстве, где расстояния Махаланобиса от доступных участков самые большие (Calenge & Basille, 2008). Поскольку все они имеют одинаковый математический смысл, эти компоненты могут быть легко объединены для создания карт пригодности местообитаний (Calenge et al., 2008). В остальном, моделирование с помощью MADIFA похоже на таковое с ENFA (Calenge, 2011).

В отличие от FANTER, идею MADIFA проще представить, если рассматривать проблему с точки зрения специалиста, хорошо знакомого с исследуемым видом. Когда известны экологические условия, требуемые для вида, изучение его ниши уже не является первостепенной целью. Вместо этого интерес представляет поиск различий между характеристиками всей среды и той ее части, которую занимает вид. Это и является смыслом моделирования пригодности местообитаний (Calenge & Basille, 2008).

Оценка качества моделей

Для каждого анализа мы проводили тестирование Монте-Карло с помощью процедуры рандомизации, чтобы выяснить значимость первого и последнего собственных значений (у) для FANTER и первого собственного значения для ENFA и MADIFA (Calenge & Basille, 2008; Caruso et al., 2015). Такое же тестирование было выполнено для оценки значимости индексов мар-гинальности и специализации (Fonderflick et al., 2015; Zimaroeva et al., 2015; Préau et al., 2018). Чтобы проверить нулевую гипотезу о случайном предпочтении бурым медведем местообитаний,

случайным образом выбирали соответствующее число точек встреч и повторяли это 999 раз. Мар-гинальность и специализация были рассчитаны для каждого набора этих точек и сравнивались со значениями, вычисленными для фактических точек присутствия бурого медведя (Basille et al., 2008; Hemery et al., 2011).

Большинство работ, посвященных ENFA, проводят оценку качества моделирования с помощью k-мерной кросс-валидации (Reutter et al., 2003; Sattler et al., 2007; Unger et al., 2008; Xuezhi et al., 2008; Qi et al., 2009), что является одним из видов зависимых оценок. Наиболее приемлемым считается совмещение кросс-валидации с независимыми оценками (Guisan & Zimmermann, 2000). Поэтому, как и в предыдущем исследовании, мы провели независимое тестирование с помощью данных, собранных отдельно от тех, на основе которых строили модель. Для этого мы использовали точки встреч бурого медведя, зарегистрированные за тот же период другими сотрудниками Центрально-Лесного заповедника (341 тестовая точка). Такое тестирование является наиболее подходящим для HSM-исследований в целом, и для GNESFA в частности (Hirzel et al., 2006; Traill & Bigalke, 2007; Braunisch et al., 2008; Falcucci et al., 2009; Valle et al., 2011; Martin et al., 2012).

Для начала рассчитывали самые простые оценки, такие как Absolute Validation Index (AVI) и Contrast Validation Index (CVI), которые очень популярны при проверке качества моделей ENFA (Hirzel et al., 2004a, 2006; Sattler et al., 2007; Pettorelli et al., 2009; Qi et al., 2009; Nazeri et al., 2014). AVI представляет собой долю точек присутствия, попадающих в категорию пригодных местообитаний выше заданного порога. Его значения варьируют от 0 до 1. CVI - это AVI, скорректированный с помощью его же теоретически ожидаемого значения, т.е. случайной модели (предсказывающей присутствие вида на всей территории исследования). Значение этого индекса лежит в диапазоне от 0.0 до 0.5 (Hirzel et al., 2006; Guisan et al., 2017). Считается, что AVI и CVI не всегда могут хорошо оценивать дискриминирующую способность моделей (Nazeri et al., 2014). Индекс Бой-са лучше показывает различия и является более точным и надежным критерием прогнозирующей способности (Boyce et al., 2002; Hirzel et al., 2006; Strubbe & Matthysen, 2008).

Индекс Бойса является показателем того, насколько предсказания модели отличаются от случайного ожидания. Он находится в диа-

пазоне от -1 до +1 и представляет собой ранговую корреляцию Спирмена между отношением предсказаний от ожиданий и средним значением каждого класса пригодности местообитаний. Положительные значения, близкие к 1, указывают на высокую прогностическую способность модели. Индекс Бойса часто используется для оценки качества предсказаний моделей ENFA и MADIFA (Hirzel et al., 2006; Sattler et al., 2007; Strubbe & Matthysen, 2008; Falcucci et al., 2009; Pettorelli et al., 2009; Qi et al., 2009; Nazeri et al., 2014; Carretero & Sillero, 2016). Для оценки качества моделей MADIFA и MD мы использовали непрерывную модификацию данного индекса по методу скользящего окна (moving window) с шириной W = 0.1 (Bcont(01), согласно Hirzel et al. (2006)). При его расчете в качестве тестовых точек использовали архивные данные.

Также оценку качества предсказаний проводили с помощью ROC-анализа (Hanley & McNeil, 1982; Brotons et al., 2004; Lobo et al., 2010; Milanesi, 2014; Préau et al., 2018; Robin et al., 2020). Мы использовали показатель площади под ROC-кривой (Area Under the Curve, AUC), рассчитанный на тестовых точках для оценки производительности моделей.

Построение карт пригодности местообитаний (HSIM)

Для получения карт пригодности местообитаний использовали алгоритмы расчета расстояний Махаланобиса, реализованные в MADIFA и MD (Calenge et al., 2008), поскольку ENFA не является математически оптимальным методом для прогнозирования (Tsoar et al., 2007; Calenge & Basille, 2008; Calenge et al., 2008). MADIFA изначально был разработан для описания ширины ниши и построения карт пригодности местообитаний с использованием меньшего числа переменных. Тем не менее, для построения таких карт лучше применять не сам MADIFA, а именно MD, поскольку этот анализ использует всю доступную информацию, а MADIFA опирается лишь на одну или несколько извлеченных осей (Halstead et al., 2010). Рассчитанные расстояния Махалано-биса могут быть переведены в вероятности, что дает пространственное представление значений HSI от 0 до 1, вычисленных для каждого пикселя (Préau et al., 2018).

Используемое программное обеспечение

Расчеты вегетационных индексов, морфо-метрических характеристик рельефа и растров

близости (расстояний до объектов) проводили в ГИС QGIS 2.18.14. Дешифрирование сцены Landsat и создание слоев типов ландшафтного покрова проводили в ГИС ArcGIS 10.4 (Esri Inc.) и QGIS с помощью плагина SCP (Semi-automatic Classification Plugin). Трехмерная модель территории подготовлена в ArcScene 10.4.

ENFA целиком реализован в отдельной программе Biomapper (Hirzel et al., 2004b), однако на сегодняшний день данный анализ проще проводить в R (Calenge & Basille, 2008; Halstead et al., 2010; Nawaz et al., 2014; Caruso et al., 2015; Андрющенко, Жуков, 2016; Préau et al., 2018). Моделирование в GNESFA выполняли с использованием пакета «adehabitatHS» (Calenge, 2011, 2020) в программе RStudio (R Development Core Team, 2020). Расчет AUC проводили при помощи пакета «pROC» (Robin et al., 2020). Цветовую гамму карт подбирали при помощи пакета «colorRamps» (Keitt, 2015). Анализ главных компонент и оценку его результатов проводили с помощью пакета «factoextra» (Kassambara & Mundt, 2020). Расчет непрерывного индекса Бойса выполняли с помощью пакета «ecospat» (Broennimann et al., 2018). Оценку информативности осей факторных анализов проводили при помощи пакета «CENFA» (Rinnan, 2020).

Все прочие статистические тесты и проверки, а также дальнейшие анализы и оценки выполняли в RStudio. Для ГИС-обработки, подготовки данных к моделированию и их визуализации были использованы следующие пакеты: «rgdal», «maptools», «gridExtra», «grid», «lattice», «reshape», «reshape2», «rasterVis», «igraph», «rgeos», «raster». Для тестов и проверок на мультиколлинеарность применяли пакеты «ecospat» и «usdm».

Результаты

Результаты моделирования экологической ниши бурого медведя в FANTER

По результатам исследования ниши бурого медведя с помощью FANTER для всех моделей удалось установить четкие паттерны и достоверность собственного значения первого фактора. Поскольку общие результаты для всех моделей оказались весьма схожи, рассмотрим лишь два первых сценария. Для модели FAHR1 первое собственное значение было достоверно (Yj = 8.30, p < 0.05) в отличие от последнего (y14 = 0.32, p = 0.89), как и для модели FAHC1 (Yj = 6.89, p < 0.05; y14 = 0.37, p = 0.90). Первый фактор модели FAHR1 объяснял 38%

информации и сильнее всего коррелировал с переменными «all_fields» (R = -0.61), «artifict» (R = 0.61) и чуть менее с «ndvi» (R = -0.47), тогда как последний - с переменными «pine_ bogs» (R = -0.77) и «coniferous» (R = -0.47; рис. 2C). Первый фактор модели FAHC1 объяснял 34% информации и сильнее всего коррелировал с переменными «all_fields» (R = -0.61), «artifict» (R = 0.61) и чуть менее с «ndvi» (R = -0.49), тогда как последний - с переменной «coniferous» (R = -0.88), а предпоследний - с «pine_bogs» (R = 0.98) (рис. 2F).

Визуализация ниши бурого медведя в экологическом пространстве наглядно проявила эвритопный характер вида, поскольку для обеих моделей реализованная ниша заняла большую часть доступного пространства (рис. 2B,E). Тем не менее, наблюдалось отчетливое вытягивание ниши влево вдоль оси первого фактора. Это вызвано большим влиянием переменных «all_fields» и «ndvi», что особенно хорошо видно при использовании весовых коэффициентов в модели FAHR1 (рис 2А).

Результаты моделирования экологической ниши бурого медведя в ENFA

По результатам моделирования ниши бурого медведя с помощью ENFA было установлено, что значение маргинальности (М) для всех моделей больше 1 (табл. 3). При этом во всех случаях прослеживается влияние весовых коэффициентов. Так, значение маргинальности всегда больше, если веса учитываются (сценарии 1, 3, 5), и всегда меньше, если точки присутствия представлены центроидами (сценарии 2, 4, 6). Значения толерантности (1/S) находились в диапазоне от 0.6 до 0.8. Только модель FAHC1 показала самое больше значение толерантности (0.79). Для всех моделей значения индексов маргинальности и толерантности были значимы (тест Монте-Карло, p < 0.05). Достоверность первого собственного значения осей специализации подтверждена только для модели FAHR1 (Sj = 3.15, p < 0.05). Для нее это означает, что дисперсия доступной среды была только в 3.15 раз больше использованной среды в измерениях первого фактора экологического пространства. Это говорит о несильной специализации вида и достаточно широкой нише. Модели не показали заметных различий по долям объясненной информации и специализации, кроме последней (FAHC3), у которой обе они были самые большие.

Рис. 2. Результаты моделирования экологической ниши бурого медведя (Ursus arctos) с помощью FANTER для моделей FAHR1 (слева) и FAHC1 (справа) для территории Центрально-Лесного заповедника и его охранной зоны за период 20082018 гг. A, D - диаграмма рассеяния доступных точек экологического пространства (серый цвет) и точек использования бурым медведем (черный цвет); размер точек для А отображает значение весового коэффициента использования. B, E - биплот экологической ниши бурого медведя; справа вверху показана диаграмма собственных значений всех факторов, слева - диаграмма их обратных значений (темным цветом выделены оставленные факторы). Светло-серый цвет - доступное экологическое пространство; темно-серый - реализованная экологическая ниша. Ось абсцисс - первая компонента; ось ординат - последняя компонента. Эколого-географические переменные спроецированы в виде векторов. С, F - корреляции переменных с осями первых двух факторов.

Fig. 2. The results of a brown bear (Ursus arctos) ecological niche modeling using FANTER for models FAHR1 (left) and FAHC1 (right) in the Central Forest State Nature Reserve and its buffer zone during 2008-2018. A - scatterplot of the ecological available space (gray colour) and the presence points of a brown bear (black colour); the size of the dots for 1A displays the value of the utilisation weight. B, E - biplot of a brown bear ecological niche; the upper right shows a diagram of the eigenvalues of all factors, on the left is a diagram of their inverse values (the left factors are highlighted in dark colour). Light gray colour - available ecological space; dark gray colour - realised ecological niche. The horizontal axis represents the first component taken, the vertical axis the last component taken. EGVs are projected as vectors. C, F - correlation between the EGVs and first two components.

Поскольку увеличение числа переменных не сказалось на достоверности собственных значений моделей, в дальнейшем мы будем рассматривать результаты ENFA только для двух первых сценариев. Первая модель объясняла 90% общей информации и 88% специализации, тогда как вторая модель - 89% информации и 86% специализации, что является хорошими показателями. Для наглядности можно рассмотреть только первые три фактора специализации, которые вместе с маргинальностью объясняли в сумме большую часть информации (56% и 52%, табл. 4). Для обеих моделей FAHR1 и FAHC1 вклад оси марги-нальности оказался существенным (22% и 19%). По своей значимости он превосходил таковой первой оси специализации. Для обеих моделей коэффициенты корреляции по оси маргиналь-ности оказались практически идентичны. Смещение ниши бурого медведя связано с выбором местообитаний с большими значениями индекса NDVI (R = 0.45 и 0.46), большим проективным

покрытием лугов (Д = 0.47 и 0.45) и меньшей удаленностью от заброшенных деревень и овсяных полей (Д = -0.66; рис. 3). Однако уже с первой оси специализации наблюдались различия между моделями. Помимо резкого снижения значения антропогенных источников пищи, обе модели демонстрировали наибольшую специализацию бурого медведя, связанную с проективным покрытием хвойных лесов ^ = 0.41 и 0.85), но только модель FAHR1 дополнительно выявила сильное влияние сфагновых сосняков (Д = 0.87). Вторая ось специализации еще больше увеличивала значение хвойных лесов для модели FAHR1 (Д = 0.86), а сфагновые сосняки уже оказывали влияние для модели FAHC1 (Д = -0.98). Третья ось специализации выделяла верховые болота (Д = 0.60 и 0.58) и расстояние до рек (Д = 0.48 и 0.59) как наиболее значимые EGV, которые сужают нишу U. arctos для обеих моделей. Первая модель выявила также некоторое влияние высоты (Д = -0.42).

Таблица 3. Результаты моделирования в ENFA для шести рассматриваемых моделей экологической ниши бурого медведя (Ursus arctos) в условиях Центрально-Лесного заповедника и его охранной зоны за период 2008-2018 гг. Table 3. The results of modelling in ENFA for the six considered models of the ecological niche of brown bear (Ursus arctos) in the conditions of the Central Forest State Nature Reserve and its buffer zone during 2008-2018

№ сценария Название модели Число осей ExpI, % ExS, % M 1/S

Y, P

1 FAHR1 9 90.4 87.7 1.17 0.69 3.15 0.01

2 FAHC1 9 88.75 86.05 1.01 0.79 2.63 0.11

3 FAHR2 11 88.57 85.66 1.35 0.62 3.88 0.15

4 FAHC2 12 90.17 88.0 1.17 0.71 3.65 0.36

5 FAHR3 14 90.45 88.16 1.45 0.62 4.28 0.15

6 FAHC3 15 91.66 89.94 1.25 0.71 3.75 0.36

EGV Модель FAHR1 Модель FAHC1

M (22%) Sj (13%) S2 (12%) S3 (9%) M (19%) Sj (13%) S2 (11%) S3 (8%)

ndvi 0.45 0.03 0.05 0.08 0.46 -0.03 -0.07 0.10

ndwi 0.11 0.10 0.07 -0.02 0.04 0.21 -0.02 0.02

elevation 0.12 -0.08 -0.23 -0.42 0.06 -0.29 0.07 -0.34

aspect -0.02 -0.06 -0.04 -0.13 -0.02 -0.05 0.04 -0.12

slope -0.01 0.02 0.02 0.02 0.04 0.02 -0.03 0.04

river -0.03 -0.01 0.08 0.48 -0.06 0.07 -0.05 0.59

coniferous -0.13 0.41 0.86 0.12 -0.14 0.85 -0.06 0.20

decidious -0.19 0.05 0.06 -0.07 -0.24 0.08 -0.11 -0.33

clearcuts 0.20 0.03 0.00 0.03 0.25 -0.02 -0.07 0.04

pine bogs -0.10 0.87 -0.36 0.21 -0.09 0.11 -0.98 -0.02

bogs -0.05 -0.08 0.09 0.60 -0.04 0.05 -0.01 0.58

all fields 0.47 -0.01 -0.07 -0.05 0.45 -0.09 -0.06 -0.01

agriculture 0.01 -0.03 -0.05 -0.34 0.03 -0.06 0.02 -0.12

artifict -0.66 -0.20 -0.19 -0.15 -0.66 -0.34 0.08 0.05

Кумулятивная информация 21.83 34.92 47.39 56.39 19.30 32.31 43.41 51.90

Примечание: Представлены первые четыре оси. Жирным шрифтом выделены коэффициенты > |0.4|. Для осей специализации (8) важны только абсолютные значения. В скобках указана доля объясненной информации.

Таблица 4. Коэффициенты корреляции (R) для эколого-географических переменных по результатам ENFA для двух моделей ниши бурого медведя (Ursus arctos) в условиях Центрально-Лесного заповедника и его охранной зоны за период 2008-2018 гг. Table 4. Correlation coefficients for ecogeographic variables according to ENFA results for the two considered models of the ecological niche of brown bear (Ursus arctos) under the conditions of the Central Forest State Nature Reserve and its buffer zone during 2008-2018

Рис. 3. Результаты моделирования экологической ниши бурого медведя (Ursus arctos) с помощью ENFA для моделей FAHR1 (слева) и FAHC1 (справа) для Центрально-Лесного заповедника и его охранной зоны за период 2008-2018 гг. A, D - диаграмма рассеяния доступных точек экологического пространства (серый цвет) и точек использования медведем (черный цвет); размер точек для А отображает значение весового коэффициента использования. B, E - биплот экологической ниши бурого медведя; справа вверху показана диаграмма собственных значений всех факторов, темным цветом выделены оставленные факторы. Светло-серый цвет - доступное экологическое пространство; темно-серый цвет - реализованная экологическая ниша. Ось абсцисс - маргинальность; ось ординат - первая ось специализации. Эколого-географические переменные спроецированы в виде векторов. С, F - корреляции переменных с осями первых двух факторов (маргинальности и первой специализации). Fig. 3. The results of a brown bear (Ursus arctos) ecological niche modeling using ENFA for models FAHR1 (left) and FAHC1 (right) in the Central Forest State Nature Reserve and its buffer zone during 2008-2018. A, D - scatterplot of the ecological available space (gray colour) and the presence points of a brown bear (black colour); the size of the dots for A displays the value of the utilisation weight. B, E - biplot of a brown bear ecological niche; the upper right shows a diagram of the eigenvalues of all factors; the left factors are highlighted in dark colour. Light gray colour - available ecological space; dark gray colour - realised ecological niche. The horizontal axis is marginality; the vertical axis is the first axis of specialisation. EGVs are projected as vectors. C, F - correlation between the EGVs and first two components.

Построенный биплот экологической ниши для двух рассматриваемых моделей позволяет наилучшим образом визуализировать выявленные взаимосвязи (рис. 3). Биплоты обеих моделей наглядно демонстрируют смещения экологической ниши бурого медведя к периферии экологического пространства по оси маргинальности (рис. 3B,E). Это происходит за счет сильного влияния переменных «all_fields» и «ndvi», что также демонстрирует график распределения точек пространства с весовыми коэффициентами (рис. 3A). Переменные «coniferous» и «pine_bogs» тесно связаны со специализацией, как мы уже убедились из значений их коэффициентов корреляции. Они растягивают область ниши вверх, вдоль первой оси специализации. Переменная «artifict» занимает промежуточное положение, имея большой вес как по маргинально-сти, так и по первой специализации. Ее вектор направлен в сторону от обеих осей, «вытягивая» левый нижний угол. Точка G по оси маргинальности обозначает центр области ниши бурого медведя, а расстояние от нее до центра осей и есть маргиналь-ность. Еще нагляднее влияния переменных проявляются на изображении их корреляций с осями маргинальности и первой специализации (рис. 3C,F). В целом обе модели проявляют схожие черты.

Основной паттерн ниши, выявленный FANTER для моделей FAHR1 и FAHC1, выражается также через маргинальность ENFA, т.е. результаты обоих анализов хорошо согласуются между собой. Это проявляется в высокой корреляции между первой компонентой FANTER и осью маргинальности ENFA (R = -0.77 для FAHR1 и R = -0.79 для FAHC1).

Результаты моделирования ниши и пригодности местообитаний бурого медведя в MADIFA и MD

Только у первой модели (FAHR1) первое и второе собственные значения оказались достоверны (Yj = 3.78, y2 = 3.14; p < 0.05), их оси вместе объясняли 32% информации. В связи с этим далее мы будем рассматривать только два первых сценария. Первые две модели MADIFA объясняли 85% и 83% информации соответственно. Для модели FAHR1 наиболее влиятельными переменными по первой компоненте оказались «coniferous» (0.74) и «artifict» (-0.63), а по второй - «pine_bogs» (0.91) (рис. 4). Для модели FAHC1 наиболее влиятельными переменными по первой компоненте также оказались «coniferous» (-0.74) и «artifict» (0.63), а по второй - «pine_bogs» (0.94) (рис. 4). В целом, результаты MADIFA хорошо согласовались с резуль-

татами предыдущих анализов. Для модели FAHR1 первая компонента MADIFA умеренно коррелировала с первой компонентой специализации ENFA (R = -0.51), а вторая практически полностью совпадала с ней (R = 0.99) и с последней компонентой FANTER (R = 0.99). Для модели FAHC1 первая компонента MADIFA сильно коррелировала с первой компонентой специализации ENFA (R = -0.85) и с последней компонентой FANTER (R = 0.91).

По результатам моделирования с помощью MD были рассчитаны основные оценочные показатели качества моделей. Так, значения AUC оказались небольшими для всех сценариев, что вполне ожидаемо для эвритопного вида. Наибольшие значения наблюдали для моделей FAHR2 и FAHC3 (0.73 и 0.72 соответственно; табл. 5). Непрерывный индекс Бойса указывал в целом на очень хорошую предсказательную способность всех моделей. Тем не менее, он был всегда больше для сценариев с использованием центроидов. Индексы AVI и CVI оказались в целом невысокими для всех моделей. Только для FAHR1 и FAHR2 CVI был сравнительно большой (0.41 и 0.43).

С помощью MD были вычислены прогнозируемые расстояния Махаланобиса для всех ячеек грида, на основании которых построены карты пригодности местообитаний (рис. 5). Все модели дали примерно один и тот же результат. Наиболее пригодные местообитания сосредоточены в охранной зоне Центрально-Лесного заповедника в восточной части территории исследования, а максимальные значения индекса пригодности растянулись вдоль долготной моренно-камовой гряды. Значение проективного покрытия верховых болот и сфагновых сосняков было различным. Использование весовых коэффициентов (модели FAHR1, FAHR2, FAHR3) в целом занижает значения HSI по территории исследования. Это наглядно проявляется для участков верховых болот и сфагновых сосняков, которые в основном выделялись в моделях с центроидами (FAHC1, FAHC2, FAHC3). Увеличение числа эколого-географических переменных в целом сужает пространство пригодных местообитаний и делает их более концентрированными.

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

Рис. 4. Результаты моделирования экологической ниши бурого медведя (Ursus arctos) с помощью MADIFA для моделей FAHR1 (слева) и FAHC1 (справа) для Центрально-Лесного заповедника и его охранной зоны за период 2008-2018 гг. A, D - диаграмма рассеяния доступных точек экологического пространства (серый цвет) и точек использования медведем (черный цвет); размер точек для А отображает значение весового коэффициента использования. B, E - биплот экологической ниши бурого медведя; справа вверху показана диаграмма собственных значений всех факторов, темным цветом выделены оставленные факторы. Светло-серый цвет - доступное экологическое пространство; темно-серый цвет - реализованная экологическая ниша. Эколого-географические переменные спроецированы в виде векторов. С, F - корреляции переменных с осями первых двух факторов.

Fig. 4. The results of a brown bear (Ursus arctos) ecological niche modeling using MADIFA for models FAHR1 (left) and FAHC1 (right) in the Central Forest State Nature Reserve and its buffer zone during 2008-2018. A, D - scatterplot of the ecological available space (gray colour) and the presence points of a brown bear (black colour); the size of the dots for A displays the value of the utilisation weight. B, E - biplot of a brown bear ecological niche; the upper right shows a diagram of the eigenvalues of all factors; the left factors are highlighted in dark colour. Light gray colour - available ecological space; dark gray colour -realised ecological niche. EGVs are projected as vectors. C, F - correlation between the EGVs and first two components.

Таблица 5. Оценки качества моделирования по расстояниям Махаланобиса для шести рассматриваемых моделей пригодности местообитаний бурого медведя (Ursus arctos) в условиях Центрально-Лесного заповедника и его охранной зоны за период 2008-2018 гг.

Table 5. Quality of modeling using Mahalanobis distances for the six considered habitat suitability models of brown bear (Ursus arctos) in the conditions of the Central Forest State Nature Reserve and its buffer zone during 2008-2018

№ сценария Название модели AVI CVI AUC B cont(0.1)

1 FAHR1 0.59 0.41 0.69 0.87

2 FAHC1 0.61 0.39 0.68 0.92

3 FAHR2 0.57 0.43 0.73 0.85

4 FAHC2 0.61 0.39 0.67 0.93

5 FAHR3 0.63 0.37 0.66 0.89

6 FAHC3 0.68 0.32 0.72 0.92

Рис. 5. Карты пригодности местообитаний бурого медведя (Ursus arctos) для территории Центрально-Лесного заповедника и его охранной зоны, построенные с помощью моделирования по расстояниям Махаланобиса (MD). Представлены шесть рассматриваемых моделей. Цветовая шкала отражает значения индексов пригодности местообитаний (HSI). Fig. 5. Habitat suitability index maps for a brown bear (Ursus arctos) in the Central Forest State Nature Reserve and its buffer zone, constructed using the Mahalanobis distances (MD) modeling. It presents six considered models. The colour scale reflects the values of habitat suitability indices (HSI).

Таблица 6. Обобщение основных EGV для ведущих факторов двух первых сценариев (1 - FAHR1; 2 - FAHC1) моделирования экологической ниши бурого медведя (Ursus arctos) в рамках GNESFA для территории Центрально-Лесного заповедника и его охранной зоны за период 2008-2018 гг.

Table 6. Main EGVs for the leading factors of the first two models (1 - FAHR1; 2 - FAHC1) of a brown bear (Ursus arctos) ecological niche modeling with the GNESFA in the Central Forest State Nature Reserve and its buffer zone during 2008-2018

№ FANTER ENFA MADIFA

1 Fj*: all fields (+); artifict (-); ndvi (+) M**: artifict (-); all fields (+); ndvi (+) F**: coniferous; artifict

F14: pine bogs; coniferous S1*: pine bogs; coniferous F2**: pine bogs;

2 Fj*: all fields (+); artifict (-); ndvi (+) M**: artifict (-); ndvi (+); all fields (+) F1: coniferous; artifict

F : coniferous S1: coniferous F2: pine bogs

Примечание: Знак в скобках указывает предпочтение переменной (+) или ее избегание (-). Еп - ведущий фактор (номер); М - маргинальному Б - специализация (номер). * - значение у для фактора достоверно (р < 0.05). ** - значение у для фактора достоверно (р < 0.01).

Обсуждение

После появления GNESFA со времен его первой публикации (Calenge & Basille, 2008) достаточно мало исследователей стало его использовать (Milanesi, 2014; Caruso et al., 2015; Milanesi et al., 2015; Neupane et al., 2019). Более того, ряд авторов продолжал использовать ENFA не в качестве исследовательского анализа, а в качестве инструмента построения HSIM (Galparsoro et al., 2009; De Angelo et al., 2011; Valle et al., 2011; Zimaroeva et al., 2015; Андрющенко, Жуков, 2016), что противоречило рекомендациям Clément Calenge и Mathieu Basille. Другие же авторы применяли только ENFA и MADIFA для моделирования пригодности местообитаний (Halstead et al., 2010; Hemery et al., 2011; Thiebot et al., 2011) вместо того, чтобы использовать комплексный подход GNESFA для более глубокого и всестороннего анализа.

Важное различие между GNESFA и стандартными методами анализа присутствия/отсутствия (например, GAM или GLM) заключается в том, что GNESFA рассчитывает пригодность местообитаний, а не вероятность присутствия вида. Это было убедительно продемонстрировано на примере моделей для 43 видов папоротников в Новой Зеландии (Zaniewski et al., 2002) и трех видов амфибий во Франции (Préau et al., 2018). В начале своего появления ENFA, в основном, использовался для построения карт пригодности местообитаний (Zaniewski et al., 2002; Reutter et al., 2003; Hirzel et al., 2004a; Chefaoui et al., 2005; Santos et al., 2006; Traill & Bigalke, 2007). Однако в дальнейшем его возможности для этого подверглись серьезной критике (Tsoar et al., 2007; Calenge & Basille 2008; Calenge et al., 2008), и на сегодняшний день это не лучший метод для построения HSIM. Наиболее адекватным применением ENFA является выявление общих закономерностей использования местообитаний на основе двух метрик экологической ниши вида - маргинальности и специализации (Basille et al., 2008).

В нашем исследовании ENFA является не самым предпочтительным методом для моделирования пригодности местообитаний еще и потому, что мы имеем дело с небольшой территорией, целиком заселенной бурым медведем. Известно, что ENFA хорошо работает для стенотопных видов или популяций с невысокой плотностью, ограниченно распределенных по территории исследования (Hirzel et al., 2002; Nawaz et al., 2014). Тем не менее, в нашем случае ENFA помогает лучшим образом визуализировать проекцию экологической ниши бурого медведя и оценить ее положение относительно пространства переменных среды. В то же время MADIFA и MD предоставляют возможности показать это в виде карт, поэтому все эти анализы следует применять совместно для достижения максимальной эффективности (Calenge & Basille, 2008; Martin et al., 2012).

В то время, как методы присутствия/отсутствия дают наиболее надежные прогнозы при наличии достоверных данных об отсутствии (Brotons et al., 2004), предсказания ENFA считаются менее надежными, хотя иногда они показывают схожие результаты (Hirzel et al., 2006). По данным Tsoar et al. (2007), предсказательная способность ENFA находится наряду с BIOCLIM на самом нижнем уровне. Так, она ниже, чем у таких методов, как HABITAT и DOMAIN, GARP и MD. В то же время, среди анализов GNESFA лучшая предсказательная способность отмечалась для MADIFA (Milanesi et al., 2015). Именно поэтому для предсказаний мы использовали MADIFA и MD (Calenge et al., 2008).

Известно, что GNESFA больше подходит для подвижных животных, таких как Canis lupus Linnaeus, 1758 (Milanesi, 2014; Milanesi et al., 2015), Puma concolor (Linnaeus, 1771) (Caruso et al., 2015) или Elephas maximus Linnaeus, 1758 (Neupane et al., 2019), для которых трудно собрать достоверные данные об отсутствии. В данной работе мы впервые продемонстрировали это и на Ursus arctos.

В нашем случае результаты всех трех анализов GNESFA хорошо согласовались между собой и выявили общие особенности экологической ниши медведя на исследуемой территории. По результатам ENFA, лучшей моделью следует признать FAHR1, поскольку она объяснила достаточно много информации и специализации, а, главное, имела достоверные значения индексов маргинальности и специализации. Применение весовых коэффициентов в целом повышало качество моделей в ENFA. Кроме того, при их использовании возрастало значение маргинальности, они способствовали его достоверности в тесте Монте-Карло и повышали долю объясненной информации.

Согласно итогам ENFA на территории исследования бурый медведь хотя и склонен проявлять себя эвритопным видом, но предпочитает существовать при определенных параметрах среды, которые отличаются от наиболее распространенных. Значения маргинальности для всех моделей были больше 1, т.е. экологическая ниша U. arctos смещалась на периферию общего доступного экологического пространства. Это говорит о том, что «средние» местообитания, предпочитаемые бурым медведем, значительно отличались от «средних» доступных местообитаний на территории исследования. В то же время значения толерантности (от 0.6 до 0.8) характеризуют бурого медведя как экологически пластичного вида, способного обитать в широком диапазоне условий. Необходимо учитывать, что поскольку в качестве точек встреч мы использовали данные, отражающие функциональные связи вида со средой (в первую очередь, пищевые), то наши результаты следует рассматривать именно с этой точки зрения.

Лучшая модель FAHR1 показала значение общей маргинальности 1.17, а толерантности 0.69. С помощью индекса маргинальности доказано значительное различие местообитаний, которые предпочитает бурый медведь, от средних местообитаний по территории исследования для разных регионов. Это продемонстрировано для Швейцарии (маргинальность - 0.95, толерантность - 0.55, согласно Zajec et al. (2005)), Апеннин (маргинальность - 0.76, толерантность - 0.52, согласно Falcucci et al. (2009)) и западных Гималаев (маргинальность - 0.80, согласно Nawaz et al. (2014)). Для пиренейской популяции U. arctos индекс маргинальности был еще выше и составил 1.2 (Martin et al., 2012). Средние значения индексов наблюдались для горных районов Греции (маргинальность - 0.55, толерантность - 0.75, со-

гласно Mertzanis et al. (2008)) и штата Кентуки, США (маргинальность - 0.51, согласно Unger et al. (2008)). Однако следует отметить, что эти коэффициенты относятся к области исследования и могут использоваться только при сравнении видов, для которых строятся модели с одним и тем же набором переменных (Hirzel et al., 2004a).

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

На материковых лугах Ursus arctos поедает травянистую растительность и раскапывает почвенные гнезда муравьев. Особое значение имеют периферии лугов и лесные поляны небольшой площади, где произрастают различные виды семейства Umbelliferae (излюбленный корм местных бурых медведей) (Ogurtsov, 2018). На пойменных лугах вдоль рек U. arctos также потребляет различные виды травянистой растительности. Заброшенные деревни предоставляют один из основных нажировочных кормов - яблоки, которые продуцируют Malus domestica L. в старых фруктовых садах. Овсяные поля, засеиваемые недалеко от таких деревень, также являются важными кормовыми стациями в летне-осенний период. В лесах группы ельников-черничников и сфагновых сосняках бурые медведи питаются плодами черники (Vaccinium myrtillus L.), а лиственные леса и зарастающие вырубки обеспечивают их ягодами рябины (Sorbus aucuparia L.) и орехами лещины (Corylus avellana L.). Более подробно значение различных местообитаний для U. arctos обсуждалось ранее (Ogurtsov, 2018; Огурцов, 2019).

Наши результаты согласуются с общими экологическими требованиями вида, которые проявляются по всему его ареалу. Луговые сообщества играют большую роль для Ursus arctos горных районов Италии (Zajec et al., 2005) и Греции, но питаться на них бурые медведи предпочитают в ночное время суток (Mertzanis et al., 2008). То же характерно для U. arctos провинции Альберта, Канада (Nielsen et al., 2004). Большое значение покрытия широколиственных лесов продемонстрировано на моделях для Центральной Италии (Posilico et al., 2004) и Греции (Mertzanis et al., 2008), где U. arctos находят нажировочные корма - плоды дуба (Quercus cerris L.) и бука (Fagus sylvatica L.). Предпочтение U. arctos небольших речных долин с богатой прибрежной раститель-

ностью также доказано для Греции (Mertzanis et al., 2008). Антропогенные ландшафты (humandominated landscapes) повсеместно определяют низкую пригодность местообитаний для бурых медведей (Zajec et al., 2005; Falcucci et al., 2009; Martin et al., 2012). Тем не менее, звери проявляют определенное тяготение к ним. Так, в горных районах Греции с незначительной плотностью сельского населения U. arctos избегают участков, близких к жилым деревням (< 500 м), но предпочитают места на промежуточных расстояниях. Эти местообитания предоставляют им ключевые пищевые ресурсы (фруктовые сады, пашни и возделанные поля с зерновыми культурами) (Mertzanis et al., 2008). Тяга U. arctos к приграничным биотопам и экотонным сообществам была продемонстрирована на других моделях ENFA, где бурые медведи предпочитали участки на границе луговых сообществ с лесом, избегая при этом больших открытых пространств (Mertzanis et al., 2008). В том же исследовании было продемонстрировано тяготение U. arctos к мозаичным ландшафтам.

Поскольку первые оси MADIFA объяснили мало информации, построенные по ним HSIM, не отражали всех отношений вида к местообитаниям. В нашем случае было целесообразно использовать всю информацию и применять MD (Halstead et al., 2010). Среди всех моделей MD лучшей оказалась последняя - FAHC3. Наряду со второй (FAHC1) она продемонстрировала самые большие значения непрерывного индекса Бойса.

Карты пригодности местообитаний оказались в целом похожими и согласуются с теми, что были построены с помощью MaxEnt (Огурцов, 2019). Наиболее пригодные местообитания во всех сценариях сосредоточены в охранной зоне Центрально-Лесного заповедника и приурочены к долготной моренно-камовой гряде, которая определяет распространение как луговых сообществ, так и заброшенных деревень - основного источника кормов для местных особей бурых медведей. В то же время наблюдались и серьезные отличия. Карты, построенные с помощью MD, выделяли больше пригодных местообитаний, чем MaxEnt, и занижали значение верховых болот.

Использование весовых коэффициентов с одной стороны повышало значение маргиналь-ности и увеличивало достоверность собственных значений, с другой - занижало значение непрерывного индекса Бойса. Во всех случаях HSIM выглядели более реалистично без использования весов. По нашему мнению, прогнозные модели и карты лучше строить без учета весовых коэффи-

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

Заключение

Мы провели комплексное моделирование пространственного аспекта экологической ниши бурого медведя. Это было выполнено с помощью трех факторных анализов в составе GNESFA. Модели самой ниши построили и изучили с использованием FANTER, ENFA и MADIFA, а модели пригодности местообитаний - с помощью анализа расстояний Махаланобиса. Несмотря на то, что метод GNESFA не является таким популярным, как MaxEnt, он способен полноценно решать требуемые задачи и давать новую информацию о характеристиках экологической ниши вида. Основные выводы о предпочтении бурым медведем местообитаний, в целом, согласуются с нашими предыдущими результатами (Огурцов, 2019). Весовые коэффициенты и число эколого-географических переменных оказали существенное влияние на результаты всех анализов. Преимущества использования GNESFA заключаются в том, что в нем сочетается как исследование экологической ниши, так и моделирование пригодности местообитаний. Помимо этого он прост в исполнении и может быть легко интерпретирован, поскольку основан на концепции экологической ниши Хатчинсона.

Благодарности

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

Дополнительная информация

Информация об использованных эколого-географических переменных и характере распределения весовых коэффициентов для моделей GNESFA (Электронное приложение. Используемые эколого-географические переменные и характер распределения весовых коэффициентов для моделей GNESFA), может быть найдена здесь.

Литература

Андрющенко А.Ю., Жуков А.В. 2016. Масштабно-зависимые эффекты в структуре экологической ниши лебедя-шипуна Cygnusolor (Gmelin, 1803) в период зимовки в пределах залива Сиваш // Бюлопчний вюник МДПУ iменi Богдана Хмельницького. Т. 6(3). C. 234-247. Демидов А.А., Кобец А.С., Грицан Ю.И., Жуков А.В. 2013. Пространственная агроэкология и рекультивация земель. Днепропетровск: Изд-во «Свидлер А.Л.». 560 с.

Желтухин А.С., Пузаченко Ю.Г., Котлов И.П., Широня И.И., Сандлерский Р.Б. 2016. Пространственно-временная динамика следовой активности лесной куницы (Martes martes L.), зайца-беляка (Lepus timidus L.) и белки (Sciurus vulgaris L.) в Европейской южной тайге // Журнал общей биологии. Т. 77(4). С. 262-283. Исаченко Т.И. 1980. Южнотаежные леса // Растительность европейской части СССР. Л.: Наука. С. 93-96. Огурцов С.С. 2016. Моделирование пригодности местообитаний бурого медведя на основе дистанционной информации // Териофауна России и сопредельных территорий. М.: Товарищество научных изданий КМК. С. 294. Огурцов С.С. 2019. Моделирование пригодности местообитаний и распределения бурого медведя (Ursus arctos) в подзоне южной тайги с помощью метода максимальной энтропии // Nature Conservation Research. Заповедная наука. Т. 4(4). С. 34-64. DOI: 10.24189/ncr.2019.061 Пузаченко Ю.Г., Козлов Д.Н. 2007. Геоморфологическая история развития территории Центрально-лесного заповедника // Труды Центрально-Лесного государственного заповедника. Вып. 4. С. 125-159. Austin М. 2007. Species distribution models and ecological theory: a critical assessment and some possible new approaches // Ecological Modelling. Vol. 200(1-2). P. 1-19. DOI: 10.1016/j.ecolmodel.2006.07.005 Basille M., Calenge C., Marboutin E., Andersen R., Gaillard J.-M. 2008. Assessing habitat selection using multivariate statistics: some refinements of the ecological-niche factor analysis // Ecological Modelling. Vol. 211(1-2). P. 233-240. DOI: 10.1016/j.ecolmodel.2007.09.006 Boyce M.S., Vernier P.R., Nielsen S.E., Schmiegelow F.K.A. 2002. Evaluating resource selection functions // Ecological Modelling. Vol. 157(2-3). P. 281-300. DOI: 10.1016/S0304-3800(02)00200-4 Braunisch V, Bollmann K., Graf R.F., Hirzel A.H. 2008. Living on the edge - Modelling habitat suitability for species at the edge of their fundamental niche // Ecological Modelling. Vol. 214(2-4). P. 153-167. DOI: 10.1016/j.ecolmodel.2008.02.001 Broennimann O. 2018. Package «ecospat», version 3.0.

Spatial Ecology Miscellaneous Methods. 107 p. Brotons L., Thuiller W., Araujo M.B., Hirzel A.H. 2004. Presence-absence versus presence-only modelling methods for predicting bird habitat suitability //

Ecography. Vol. 27(4). P. 437-448. DOI: 10.1111/j.0906-7590.2004.03764.x Browning D., Beaupré S., Duncan L. 2005. Using partitioned Mahalanobis D2(K) to formulate a GIS-based model of timber rattlesnake hibernacula // Journal of Wildlife Management. Vol. 69(1). P. 3344. DOI: 10.2193/0022-541X(2005)069<0033:UPM DTF>2.0.C0;2 Calenge C. 2011. Exploratory analysis of the habitat selection by the wildlife in R: the «adehabitatHS» package. 60 p. Calenge C. 2020. Package «adehabitatHS», version 0.3.15.

Analysis of Habitat Selection by Animals. 67 p. Calenge C., Basille M. 2008. A general framework for the statistical exploration of the ecological niche // Journal of Theoretical Biology. Vol. 252(4). P. 674-685. DOI: 10.1016/j.jtbi.2008.02.036 Calenge C., Darmon G., Basille M., Loison A., Jullien J.-M. 2008. The factorial decomposition of the Mahalanobis distances in habitat selection studies // Ecology. Vol. 89(2). P. 555-566. DOI: 10.1890/06-1750.1 Carretero M.A., Sillero N. 2016. Evaluating how species niche modelling is affected by partial distributions with an empirical case // Acta Oecologica. Vol. 77. P. 207216. DOI: 10.1016/j.actao.2016.08.014 Caruso N., Guerisoli M., Luengos Vidal E.M., Castillo D., Casanave E.B., Lucherini M. 2015. Modelling the ecological niche of an endangered population of Puma concolor: First application of the GNESFA method to an elusive carnivore // Ecological Modelling. Vol. 297. P. 11-19. DOI: 10.1016/j.ecolmodel.2014.11.004 Cayuela L. 2004. Habitat evaluation for the Iberian wolf Canis lupus in Picos de Europa National Park, Spain // Applied Geography. Vol. 24(3). P. 199-215. DOI: 10.1016/j.apgeog.2004.04.003 Chefaoui R.M., Hortal J., Lobo J.M. 2005. Potential distribution modelling, niche characterization and conservation status assessment using GIS tools: a case study of Iberian Copris species // Biological Conservation. Vol. 122(2). P. 327-338. DOI: 10.1016/j. biocon.2004.08.005 Cherednichenko O., Gorik V., Borodulina V. 2016. Herb vegetation diversity in the north of Central Forest Reserve (Tver province, Russia) // 25th Meeting of European Vegetation Survey. Book of Abstracts. Italy, Roma: Sapienza University. P. 26. Clark J., Dunn J., Smith K. 1993. A multivariate model of female black bear habitat use for a geographic information system // Journal of Wildlife Management. Vol. 57(3). P. 519-526. DOI: 10.2307/3809276 Corsi F., Dupré E., Boitani L. 1999. A Large-Scale Model of Wolf Distribution in Italy for Conservation Planning // Conservation Biology. Vol. 13(1). P. 150-159. DOI: 10.1046/j.1523-1739.1999.97269.x De Angelo C., Paviolo A., Di Bitetti M. 2011. Differential impact of landscape transformation on pumas (Puma concolor) and jaguars (Panthera onca) in the Upper Paraná Atlantic Forest // Diversity and Distributions. Vol. 17(3). P. 422-436. DOI: 10.1111/j.1472-4642.2011.00746.x

Dettki H., Löfstrand R., Edenius L. 2003. Modeling habitat suitability for moose in coastal northern Sweden: empirical vs. process-oriented approaches // Ambio. Vol. 32(8). P. 549- 556. DOI: 10.1639/0044-7447(2003)032[0549:MHS FMI]2.0.CO;2

Dormann C.F., Elith J., Bacher S., Buchmann C., Carl G., Carré G., García Marquéz J.R., Gruber B., Lafourcade B., Leitäo P.J., Münkemüller T., McClean C., Osborne P.E., Reineking B., Schröder B., Skidmore A.K., Zurell D., Lautenbach S. 2013. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance // Ecography. Vol. 36(1). P. 27-46. DOI: 10.1111/j.1600-0587.2012.07348.x Dray S., Chessel D., Thioulouse J. 2003. Co-inertia analysis and the linking of ecological data tables // Ecology. Vol. 84(11). P. 3078-3089. DOI: 10.1890/03-0178 Elith J., Graham C.H., Anderson P.R., Dudík M., Ferrier S., Guisan A., Hijmans R.J., Huettmann F., Leathwick J.R., Lehmann A., Li J., Lohmann L.G., Loiselle B.A., Manion G., Moritz C., Nakamura M., Nakazawa Y., Overton J.McC.M., Townsend Peterson A., Phillips S.J., Richardson K., Scachetti-Pereira R.K., Schapire R.E., Soberón J., Williams S., Wisz M.S., Zimmermann N.E. 2006. Novel methods improve prediction of species' distributions from occurrence data // Ecography. Vol. 29(2). P. 129-151. DOI: 10.1111/j.2006.0906-7590.04596.x Falcucci A., Ciucci P., Maiorano L., Gentile L., Boitani L. 2009. Assessing habitat quality for conservation using an integrated occurrence-mortality model // Journal of Applied Ecology. Vol. 46(3). P. 600-609. DOI: 10.1111/j.1365-2664.2009.01634.x Farber O., Kadmon R. 2003. Assessment of alternative approaches for bioclimatic modeling with special emphasis on the Mahalanobis distance // Ecological Modelling. Vol. 160(1-2). P. 115-130. DOI: 10.1016/ S0304-3800(02)00327-7 Fonderflick J., Azam C., Brochier C., Cosson E., Quékenborn D. 2015. Testing the relevance of using spatial modeling to predict foraging habitat suitability around bat maternity: a case study in Mediterranean landscape // Biological Conservation. Vol. 192. P. 120-129. DOI: 10.1016/j.biocon.2015.09.012 Gallego D., Cánovas F., Esteve M.A., Galían J. 2004. Descriptive biogeography of Tomicus (Coleoptera: Scolytidae) species in Spain // Journal of Biogeography. Vol. 31(12). P. 2011-2024. DOI: 10.1111/j.1365-2699.2004.01131.x Galparsoro I., Borja Á., Bald J., Liria P., Chust G. 2009. Predicting suitable habitat for the European lobster (Homarus gammarus), on the Basque continental shelf (Bay of Biscay), using Ecological-Niche Factor Analysis // Ecological Modelling. Vol. 220(4). P. 556567. DOI: 10.1016/j.ecolmodel.2008.11.003 Guisan A., Thuiller W., Zimmermann N.E. 2017. Habitat suitability and distribution models. Cambridge: Cambridge University Press. 462 p. Guisan A., Zimmermann N.E. 2000. Predictive habitat distribution models in ecology // Ecological Modelling.

Vol. 135(2-3). P. 147-186. DOI: 10.1016/S0304-3800(00)00354-9 Halstead B.J., Wylie G.D., Casazza M.L. 2010. Habitat Suitability and Conservation of the Giant Gartersnake (Thamnophis gigas) in the Sacramento Valley of California // Copeia. Vol. 4. P. 591-599. DOI: 10.1643/CE-09-199 Hanley J.A., McNeil B.J. 1982. The meaning and use of the area under a receiver operating characteristic (ROC) curve // Radiology. Vol. 143(1). P. 29-36. DOI: 10.1148/ radiology. 143.1.7063747 Hemery L., Galton-Fenzi B., Ameziane N., Riddle M., Rintoul S., Beaman R., Post A., Eleaume M. 2011. Predicting habitat preferences for Anthometrina adriani (Echinodermata) on the East Antarctic continental shelf // Marine Ecology Progress Series. Vol. 441. P. 105116. DOI: 10.3354/meps09330 Hirzel A.H., Hausser J., Chessel D., Perrin N. 2002. Ecological-niche factor analysis: how to compute habitat-suitability maps without absence data? // Ecology. Vol. 83(7). P. 2027-2036. DOI: 10.1890/0012-9658(2002)083[2027:EN FAHT]2.0.CO;2 Hirzel A.H., Posse B., Oggier P.A., Crettenand Y., Glenz C., Arlettaz R. 2004a. Ecological requirements of reintroduced species and the implications for release policy: the case of the bearded vulture // Journal of Applied Ecology. Vol. 41(6). P. 1103-1116. DOI: 10.1111/j.0021-8901.2004.00980.x Hirzel A.H., Hausser J., Perrin N. 2004b. Biomapper 3.1. Division of Conservation Biology. Bern: University of Bern. Available from http://www.unil.ch/biomapper Hirzel A.H., Le Lay G., Helfer V, Randin C., Guisan A. 2006. Evaluating the ability of habitat suitability models to predict species presences // Ecological Modelling. Vol. 199(2). P. 142-152. DOI: 10.1016/j. ecolmodel.2006.05.017 Hirzel A.H., Le Lay G. 2008. Habitat suitability modelling and niche theory // Journal of Applied Ecology. Vol. 45(5). P. 1372-1381. DOI: 10.1111/j.1365-2664.2008.01524.x Hutchinson G.E. 1957. Concluding remarks // Cold Spring Harbour Symposium on Quantitative Biology. Vol. 22. P. 415-427. DOI: 10.1101/SQB.1957.022.01.039 Jackson D.A. 1993. Stopping rules in principal components analysis: a comparison of heuristical and statistical approaches // Ecology. Vol. 74(8). P. 2204-2214. DOI: 10.2307/1939574 Kassambara A., Mundt F. 2020. Package «factoextra», version 1.0.7. Extract and Visualize the Results of Multivariate Data Analyses. 84 p. Keitt T. 2015. Package «colorRamps», version 2.3. Builds

color tables. 9 p. Knick S., Dyer D. 1997. Distribution of black-tailed jackrabbit habitat determined by GIS in southwestern Idaho // Journal of Wildlife Management. Vol. 61(1). P. 75-85. DOI: 10.2307/3802416 Lobo J.M., Jimenez-Valverde A., Hortal J. 2010. The uncertain nature of absences and their importance in species distribution modelling // Ecography. Vol. 33(1). P. 103-114. DOI: 10.1111/j.1600-0587.2009.06039.x

Manly B.F., McDonald L.L., Thomas D.L., MacDonald T.L., Erickson W.P. 2002. Resource selection by animals. Statistical design and analysis for field studies. London: Kluwer Academic Publisher. 222 p. DOI: 10.1007/0-306-48151-0 Martin J., Revilla E., Quenette P.-Y., Naves J., Allainé D., Swenson J.E. 2012. Brown bear habitat suitability in the Pyrenees: transferability across sites and linking scales to make the most of scarce data // Journal of Applied Ecology. Vol. 49(3). P. 621-631. DOI: 10.1m/j.1365-2664.2012.02139.x Mertzanis G., Kallimanis A.S., Kanellopoulos N., Sgardelis S.P., Tragos A., Aravidis I. 2008. Brown bear (Ursus arctos L.) habitat use patterns in two regions of northern Pindos, Greece - management implications // Journal of Natural History. Vol. 42(5-8). P. 301-315. DOI: 10.1080/00222930701835175 Mestre F.M., Ferreira J.P., Mira A. 2007. Modelling the distribution of the european polecat Mustela putorius in a Mediterranean agricultural landscape // Revue d Écologie (Terre Vie). Vol. 62(1). P. 35-47. Milanesi P. 2014. Landscape ecology and genetics of the wolf in Italy. PhD Thesis. Bologna: Università di Bologna. 126 p. Milanesi P., Caniglia R., Fabbri E., Galaverni M., Meriggi A., Randi E. 2015. Non-invasive genetic sampling to predict wolf distribution and habitat suitability in the Northern Italian Apennines: implications for livestock depredation risk // European Journal of Wildlife Research. Vol. 61(5). P. 681-689. DOI: 10.1007/s 10344-015-0942-4 Nawaz M.A., Martin J., Swenson J.E. 2014. Identifying key habitats to conserve the threatened brown bear in the Himalaya // Biological Conservation. Vol. 170. P. 198206. DOI: 10.1016/j.biocon.2013.12.031 Nazeri M., Kumar L., Jusoff K., Bahaman A.R. 2014. Modeling the potential distribution of sun bear in Krau wildlife reserve, Malaysia // Ecological Informatics. Vol. 20. P. 27-32. DOI: 10.1016/j. ecoinf.2014.01.006 Neupane D., Kwon Y., Risch T.S., Williams A.C., Johnson R.L. 2019. Habitat use by Asian elephants: Context matters // Global Ecology and Conservation. Vol. 17. Article: e00570. DOI: 10.1016/j.gecco.2019.e00570 Nielsen S.E., Boyce M.S., Stenhouse G.B. 2004. Grizzly bears and forestry: I. Selection of clearcuts by grizzly bears in west-central Alberta, Canada // Forest Ecology and Management. Vol. 199(1). P. 51-65. DOI: 10.1016/j. foreco.2004.04.015 Ogurtsov S.S. 2018. The Diet of the Brown Bear (Ursus arctos) in the Central Forest Nature Reserve (West-European Russia), Based on Scat Analysis Data // Biology Bulletin. Vol. 45(9). P. 1039-1054. DOI: 10.1134/S 1062359018090145 Peterson A.T., Soberon J., Pearson R.G., Anderson R.P., Martinez-Meyer E., Nakamura M., Araùjo M.B. 2011. Ecological Niches and Geographic Distributions. Princeton and Oxford: Princeton University Press. 329 p.

PettoreUi N., Hilborn A., Broekhuis F., Durant S.M. 2009. Exploring habitat use by cheetahs using ecological niche factor analysis // Journal of Zoology. Vol. 277(2). P. 141-148. DOI: 10.1111/j.1469-7998.2008.00522.x Phillips S.J., Anderson R.P., Schapire R.E. 2006. Maximum entropy modeling of species geographic distributions // Ecological Modelling. Vol. 190(3-4). P. 231-259. DOI: 10.1016/j.ecolmodel.2005.03.026 Posilico M., Meriggi A., Pagnin E., Lovari S., Russo L. 2004. A habitat model for brown bear conservation and land use planning in the central Apennines // Biological Conservation. Vol. 118(2). P. 141-150. DOI: 10.1016/j. biocon.2003.07.017 Préau C., Trochet A., Bertrand R., Isselin-Nondedeu F. 2018. Modeling Potential distributions of three european amphibian species comparing ENFA and MaxEnt // Herpetological Conservation and Biology. Vol. 13(1). P. 91-104.

Qi D., Hu Y., Gu X., Li M., Wei F. 2009. Ecological niche modeling of the sympatric giant and red pandas on a mountain-range scale // Biodiversity and Conservation. Vol. 18(8). P. 2127-2141. DOI: 10.1007/s10531-009-9577-7 R Development Core Team. 2020. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Available from http://www.R-project.org Reutter B.A., Helfer V, Hirzel A.H., Vogel P. 2003. Modelling habitat-suitability using museum collections: an example with three sympatric Apodemus species from the Alps // Journal of Biogeography. Vol. 30(4). P. 581590. DOI: 10.1046/j.1365-2699.2003.00855.x Rinnan D.S., Lawler J. 2019. Climate-niche factor analysis: a spatial approach to quantifying species vulnerability to climate change // Ecography. Vol. 42(9). P. 1494-1503. DOI: 10.1111/ecog.03937 Rinnan D.S. 2020. Package «CENFA», version 1.1.0. Climate

and Ecological Niche Factor Analysis. 36 p. Robin X., Turck N., Hainard A., Tiberti N., Lisacek F., Sanchez J.-C., Müller M., Siegert S., Doering M. 2020. Package «pROC», version 1.16.2. Display and Analyze ROC Curves. 95 p. Rotenberry J.T., Preston K.L., Knick S.T. 2006. GIS-based niche modeling for mapping species' habitat // Ecology. Vol. 87(6). P. 1458-1464. DOI: 10.1890/0012-9658(2006)87[1458:GNMFMS]2.0.œ;2 Santos X., Brito J.C., Sillero N., Pleguezuelos J.M., Llorente G., Fahd S., Parellada X. 2006. Inferring habitat-suitability areas with ecological modelling techniques and GIS: A contribution to assess the conservation status of Vipera latastei // Biological Conservation. Vol. 130(3). P. 416-425. DOI: 10.1016/j.biocon.2006.01.003 Sattler T., Bontadina F., Hirzel A.H., Arlettaz R. 2007. Ecological niche modelling of two cryptic bat species calls for a reassessment of their conservation status // Journal of Applied Ecology. Vol. 44(6). P. 1188-1199. DOI: 10.1111/j.1365-2664.2007.01328.x Silva L.D., Costa H., Azevedo E., Medeiros V., Alves M., Elias R.B., Silva L. 2017. Modelling native and invasive

woody species: a comparison of ENFA and MaxEnt applied to the Azorean forest // Modeling, Dynamics, Optimization and Bioeconomics II. DGS 2014 / A. Pinto, D. Zilberman (Eds.). Vol. 195. P. 415-444. DOI: 10.1007/978-3-319-55236-1_20 Strubbe D., Matthysen E. 2008. Predicting the potential distribution of invasive ring-necked parakeets Psittacula krameri in northern Belgium using an ecological niche modelling approach // Biological Invasions. Vol. 11(3). P. 497-513. DOI: 10.1007/s 10530-008-9266-6 Thiebot J.-B., Lescroel A., Pinaud D., Trathan P.N., Bost C.-A. 2011. Larger foraging range but similar habitat selection in non-breeding versus breeding sub-Antarctic penguins // Antarctic Science. Vol. 23(2). P. 117-126. DOI: 10.1017/S0954102010000957 Thioulouse J., Chessel D. 1992. A method for reciprocal scaling of species tolerance and sample diversity // Ecology. Vol. 73(2). P. 670-680. DOI: 10.2307/1940773 Traill L.W., Bigalke R.C. 2007. A presence-only habitat suitability model for large grazing African ungulates and its utility for wildlife management // African Journal of Ecology. Vol. 45(3). P. 347-354. DOI: 10.1111/j.1365-2028.2006.00717.x Tsoar A., Allouche O., Steinitz O., Rotem D., Kadmon R. 2007. A comparative evaluation of presence-only methods for modelling species distribution // Diversity and Distribution. Vol. 13(4). P. 397-405. DOI: 10.1111/j.1472-4642.2007.00346.x Unger D.E., Fei S., Maehr D.S. 2008. Ecological niche factor analysis to determine habitat suitability of a recolonizing carnivore // Proceedings of the 6th Southern Forestry and Natural Resources GIS Conference. Athens, GA: University of Georgia. P. 237-250. Valle M., Boija A., Chust G., Galparsoro I., Garmendia J.M. 2011. Modelling suitable estuarine habitats for Zostera noltii, using Ecological Niche Factor Analysis and Bathymetric LiDAR // Estuarine Coastal and Shelf Science. Vol. 94(2). P. 144-154. DOI: 10.1016/j. ecss.2011.05.031 Xuezhi W., Weihua X., Zhiyun O., Jianguo L., Yi X., Youping C., Lianjun Z., Junzhong H. 2008. Application of ecological-niche factor analysis in habitat assessment of giant pandas // Acta Ecologica Sinica. Vol. 28(2). P. 821-828. DOI: 10.1016/S1872-2032(08)60030-X Zajec P., Zimmermann F., Roth H.U., Breitenmoser U. 2005. The return of the Brown bear to Switzerland - Suitable habitat distribution, corridors and potential conflicts // KORA Bericht. Vol. 28e. P. 1-31. Zaniewski A.E., Lehmann A., Overton J.McC. 2002. Predicting species spatial distributions using presence-only data: a case study of native New Zealand ferns // Ecological Modelling. Vol. 157(2-3). P. 261-280. DOI: 10.1016/S0304-3800(02)00199-0 Zimaroeva A.A., Zhukov A.V., Ponomarenko A.L., Matsyura A. 2015. Ecological niche modelling of Fringilla coelebs Linnaeus, 1758 (common chaffinch) using GIS tools // Romanian Journal of Biology - Zoology. Vol. 60(2). P. 135-146.

Zimmermann F. 2004. Conservation of the Eurasian lynx (Lynx lynx) in a Fragmented Landscape - Habitat Models, Dispersal and Potential Distribution. PhD Thesis. Lausanne: University of Lausanne. 193 p.

References

Andryuschenko A.Yu., Zhukov A.V. 2016. Scale-dependent effects in structure of the wintering ecological niche of the Mute Swan during wintering in the Gulf of Sivash. Biological Bulletin of Bogdan Khmelnitskiy Melitopol State Pedagogical University 6(3): 234-247. [In Russian] Austin M. 2007. Species distribution models and ecological theory: a critical assessment and some possible new approaches. Ecological Modelling 200(1-2): 1-19. DOI: 10.1016/j.ecolmodel.2006.07.005 Basille M., Calenge C., Marboutin É., Andersen R., Gaillard J.-M. 2008. Assessing habitat selection using multivariate statistics: some refinements of the ecological-niche factor analysis. Ecological Modelling 211(1-2): 233240. DOI: 10.1016/j.ecolmodel.2007.09.006 Boyce M.S., Vernier P.R., Nielsen S.E., Schmiegelow F.K.A. 2002. Evaluating resource selection functions. Ecological Modelling 157(2-3): 281-300. DOI: 10.1016/S0304-3800(02)00200-4 Braunisch V, Bollmann K., Graf R.F., Hirzel A.H. 2008. Living on the edge - Modelling habitat suitability for species at the edge of their fundamental niche. Ecological Modelling 214(2-4): 153-167. DOI: 10.1016/j.ecolmodel.2008.02.001 Broennimann O. 2018. Package «ecospat», version 3.0.

Spatial Ecology Miscellaneous Methods. 107 p. Brotons L., Thuiller W., Araujo M.B., Hirzel A.H. 2004. Presence-absence versus presence-only modelling methods for predicting bird habitat suitability. Ecography 27(4): 437-448. DOI: 10.1111/j.0906-7590.2004.03764.x Browning D., Beaupré S., Duncan L. 2005. Using partitioned Mahalanobis D2(K) to formulate a GIS-based model of timber rattlesnake hibernacula. Journal of Wildlife Management 69(1): 33-44. DOI: 10.2193/0022-541X( 2005)069<0033:UPMDTF>2.0.œ;2 Calenge C. 2011. Exploratory analysis of the habitat selection by the wildlife in R: the «adehabitatHS» package. 60 p. Calenge C. 2020. Package «adehabitatHS», version 0.3.15.

Analysis of Habitat Selection by Animals. 67 p. Calenge C., Basille M. 2008. A general framework for the statistical exploration of the ecological niche. Journal of Theoretical Biology 252(4): 674-685. DOI: 10.1016/j. jtbi.2008.02.036 Calenge C., Darmon G., Basille M., Loison A., Jullien J.-M. 2008. The factorial decomposition of the Mahalanobis distances in habitat selection studies. Ecology 89(2): 555-566. DOI: 10.1890/06-1750.1 Carretero M.A., Sillero N. 2016. Evaluating how species niche modelling is affected by partial distributions with an empirical case. Acta Oecologica 77: 207-216. DOI: 10.1016/j.actao.2016.08.014

Caruso N., Guerisoli M., Luengos Vidal E.M., Castillo D., Casanave E.B., Lucherini M. 2015. Modelling the ecological niche of an endangered population of Puma concolor: First application of the GNESFA method to an elusive carnivore. Ecological Modelling 297: 11-19. DOI: 10.1016/j. ecolmodel.2014.11.004 Cayuela L. 2004. Habitat evaluation for the Iberian wolf Canis lupus in Picos de Europa National Park, Spain. Applied Geography 24(3): 199-215. DOI: 10.1016/j. apgeog.2004.04.003 Chefaoui R.M., Hortal J., Lobo J.M. 2005. Potential distribution modelling, niche characterization and conservation status assessment using GIS tools: a case study of Iberian Copris species. Biological Conservation 122(2): 327-338. DOI: 10.1016/j.biocon.2004.08.005 Cherednichenko O., Gorik V., Borodulina V. 2016. Herb vegetation diversity in the north of Central Forest Reserve (Tver' province, Russia). In: 25th Meeting of European Vegetation Survey. Book of Abstracts. Italy, Roma: Sapienza University. P. 26. Clark J., Dunn J., Smith K. 1993. A multivariate model of female black bear habitat use for a geographic information system. Journal of Wildlife Management 57(3): 519-526. DOI: 10.2307/3809276 Corsi F., Dupré E., Boitani L. 1999. A large scale model of wolf distribution in Italy for conservation planning. Conservation Biology 13(1): 150-159. DOI: 10.1046/j.1523-1739.1999.97269.x De Angelo C., Paviolo A., Di Bitetti M. 2011. Differential impact of landscape transformation on pumas (Puma concolor) and jaguars (Panthera onca) in the Upper Paraná Atlantic Forest. Diversity and Distributions 17(3): 422-436. DOI: 10.1111/j.1472-4642.2011.00746.x Demidov A.A., Kobets A.S., Gritsan Yu.I., Zhukov A.V. 2013. Spatial agricultural ecology and soil recultivation. Dnepropetrovsk: A.L. Svidler Press. 560 p. [In Russian]

Dettki H., Löfstrand R., Edenius L. 2003. Modeling habitat suitability for moose in coastal northern Sweden: empirical vs. process-oriented approaches. Ambio 32(8): 549-556.DOI: 10.1639/0044-7447(2003)032[0549:MH SFMI]2.0.CO;2 Dormann C.F., Elith J., Bacher S., Buchmann C., Carl G., Carré G., García Marquéz J.R., Gruber B., Lafourcade B., Leitâo P.J., Münkemüller T., McClean C., Osborne P.E., Reineking B., Schröder B., Skidmore A.K., Zurell D., Lautenbach S. 2013. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography 36(1): 27-46. DOI: 10.1111/ j.1600- 0587.2012.07348.x Dray S., Chessel D., Thioulouse J. 2003. Co-inertia analysis and the linking of ecological data tables. Ecology 84(11): 3078-3089. DOI: 10.1890/03-0178 Elith J., Graham C.H., Anderson P.R., Dudík M., Ferrier S., Guisan A., Hijmans R.J., Huettmann F., Leathwick J.R., Lehmann A., Li J., Lohmann L.G., Loiselle B.A., Manion G., Moritz C., Nakamura M., Nakazawa Y., Overton J.McC.M., Townsend Peterson A., Phillips S.J.,

Richardson K., Scachetti-Pereira R.K., Schapire R.E., Soberón J., Williams S., Wisz M.S., Zimmermann N.E. 2006. Novel methods improve prediction of species' distributions from occurrence data. Ecography 29(2): 129-151. DOI: 10.1111/j.2006.0906-7590.04596.x Falcucci A., Ciucci P., Maiorano L., Gentile L., Boitani L. 2009. Assessing habitat quality for conservation using an integrated occurrence-mortality model. Journal of Applied Ecology 46(3): 600-609. DOI: 10.1111/j.1365-2664.2009.01634.x Farber O., Kadmon R. 2003. Assessment of alternative approaches for bioclimatic modeling with special emphasis on the Mahalanobis distance. Ecological Modelling 160(1-2): 115-130. DOI: 10.1016/S0304-3800(02)00327-7 Fonderflick J., Azam C., Brochier C., Cosson E., Quékenborn D. 2015. Testing the relevance of using spatial modeling to predict foraging habitat suitability around bat maternity: a case study in Mediterranean landscape. Biological Conservation 192: 120-129. DOI: 10.1016/j. biocon.2015.09.012 Gallego D., Cánovas F., Esteve M.A., Galían J. 2004. Descriptive biogeography of Tomicus (Coleoptera: Scolytidae) species in Spain. Journal of Biogeography 31(12): 2011-2024. DOI: 10.1111/j.1365-2699.2004.01131.x Galparsoro I., Borja Á., Bald J., Liria P., Chust G. 2009. Predicting suitable habitat for the European lobster (Homarus gammarus), on the Basque continental shelf (Bay of Biscay), using Ecological-Niche Factor Analysis. Ecological Modelling 220(4): 556-567. DOI: 10.1016/j. ecolmodel.2008.11.003 Guisan A., Thuiller W., Zimmermann N.E. 2017. Habitat suitability and distribution models. Cambridge: Cambridge University Press. 462 p. Guisan A., Zimmermann N.E. 2000. Predictive habitat distribution models in ecology. Ecological Modelling 135(2-3): 147-186. DOI: 10.1016/S0304-3800(00)00354-9 Halstead B.J., Wylie G.D., Casazza M.L. 2010. Habitat Suitability and Conservation of the Giant Gartersnake (Thamnophis gigas) in the Sacramento Valley of California. Copeia 4: 591-599. DOI: 10.1643/CE-09-199 Hanley J.A., McNeil B.J. 1982. The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology 143(1): 29-36. DOI: 10.1148/ radiology. 143.1.7063747 Hemery L., Galton-Fenzi B., Améziane N., Riddle M., Rintoul S., Beaman R., Post A., Eléaume M. 2011. Predicting habitat preferences for Anthometrina adriani (Echinodermata) on the East Antarctic continental shelf. Marine Ecology Progress Series 441: 105-116. DOI: 10.3354/meps09330 HirzelA.H., Hausser J., Chessel D., Perrin N. 2002. Ecological-niche factor analysis: how to compute habitat-suitability maps without absence data? Ecology 83(7): 2027-2036. DOI: 10.1890/0012-9658(2002)083 [2027:ENFAHT]2. 0.CO;2

Hirzel A.H., Posse B., Oggier P.A., Crettenand Y., Glenz C., Arlettaz R. 2004a. Ecological requirements of reintroduced species and the implications for release policy: the case of the bearded vulture. Journal of Applied Ecology 41(6): 1103-1116. DOI: 10.im/j.0021-890L2004.00980.x Hirzel A.H., Hausser J., Perrin N. 2004b. Biomapper 3.1. Division of Conservation Biology. Bern: University of Bern. Available from http://www.unil.ch/biomapper Hirzel A.H., Le Lay G., Helfer V, Randin C., Guisan A. 2006. Evaluating the ability of habitat suitability models to predict species presences. Ecological Modelling 199(2): 142-152. DOI: 10.1016/j.ecolmodel.2006.05.017 Hirzel A.H., Le Lay G. 2008. Habitat suitability modelling and niche theory. Journal of Applied Ecology 45(5): 1372-1381. DOI: 10.1111/j.1365-2664.2008.01524.x Hutchinson G.E. 1957. Concluding remarks. Cold Spring Harbour Symposium on Quantitative Biology 22: 415427. DOI: 10.1101/SQB.1957.022.01.039 Isachenko T.I. 1980. South taiga forests. In: Vegetation of the European part of the USSR. Leningrad: Nauka. P. 93-96. [In Russian] Jackson D.A. 1993. Stopping rules in principal components analysis: a comparison of heuristical and statistical approaches. Ecology 74(8): 2204-2214. DOI: 10.2307/1939574 Kassambara A., Mundt F. 2020. Package «factoextra», version 1.0.7. Extract and Visualize the Results of Multivariate Data Analyses. 84 p. Keitt T. 2015. Package «colorRamps», version 2.3. Builds

color tables. 9 p. Knick S., Dyer D. 1997. Distribution of black-tailed jackrabbit habitat determined by GIS in southwestern Idaho. Journal of Wildlife Management 61(1): 75-85. DOI: 10.2307/3802416 Lobo J.M, Jiménez-Valverde A., Hortal J. 2010. The uncertain nature of absences and their importance in species distribution modelling. Ecography 33(1): 103114. DOI: 10.1111/j.1600-0587.2009.06039.x Manly B.F., McDonald L.L., Thomas D.L., MacDonald T.L., Erickson W.P. 2002. Resource selection by animals. Statistical design and analysis for field studies. London: Kluwer Academic Publisher. 222 p. DOI: 10.1007/0-306-48151-0 Martin J., Revilla E., Quenette P.-Y., Naves J., Allainé D., Swenson J.E. 2012. Brown bear habitat suitability in the Pyrenees: transferability across sites and linking scales to make the most of scarce data. Journal of Applied Ecology 49(3): 621-631. DOI: 10.1111/j.1365-2664.2012.02139.x Mertzanis G., Kallimanis A.S., Kanellopoulos N., Sgardelis S.P., Tragos A., Aravidis I. 2008. Brown bear (Ursus arctos L.) habitat use patterns in two regions of northern Pindos, Greece - management implications. Journal of Natural History 42(5-8): 301-315. DOI: 10.1080/00222930701835175 Mestre F.M., Ferreira J.P., Mira A. 2007. Modelling the distribution of the european polecat Mustela putorius in a Mediterranean agricultural landscape. Revue d Écologie (Terre Vie) 62(1): 35-47.

Milanesi P. 2014. Landscape ecology and genetics of the wolf in Italy. PhD Thesis. Bologna: Università di Bologna. 126 p.

Milanesi P., Caniglia R., Fabbri E., Galaverni M., Meriggi A., Randi E. 2015. Non-invasive genetic sampling to predict wolf distribution and habitat suitability in the Northern Italian Apennines: implications for livestock depredation risk. European Journal of Wildlife Research 61(5): 681-689. DOI: 10.1007/s10344-015-0942-4 Nawaz M.A., Martin J., Swenson J.E. 2014. Identifying key habitats to conserve the threatened brown bear in the Himalaya. Biological Conservation 170: 198-206. DOI: 10.1016/j.biocon.2013.12.031 Nazeri M., Kumar L., Jusoff K., Bahaman A.R. 2014. Modeling the potential distribution of sun bear in Krau wildlife reserve, Malaysia. Ecological Informatics 20: 27-32. DOI: 10.1016/j.ecoinf.2014.01.006 Neupane D., Kwon Y., Risch T.S., Williams A.C., Johnson R.L. 2019. Habitat use by Asian elephants: Context matters. Global Ecology and Conservation 17: e00570. DOI: 10.1016/j.gecco.2019.e00570 Nielsen S.E., Boyce M.S., Stenhouse G.B. 2004. Grizzly bears and forestry: I. Selection of clearcuts by grizzly bears in west-central Alberta, Canada. Forest Ecology and Management 199(1): 51-65. DOI: 10.1016/j. foreco.2004.04.015 Ogurtsov S.S. 2016. Brown bear habitat suitability modelling based on remote sensing data. In: Theriofauna of Russia and adjacent areas. Moscow: KMK Scientific Press Ltd. P. 294. [In Russian] Ogurtsov S.S. 2018. The Diet of the Brown Bear (Ursus arctos) in the Central Forest Nature Reserve (West-European Russia), Based on Scat Analysis Data. Biology Bulletin 45(9): 1039-1054. DOI: 10.1134/ S1062359018090145 Ogurtsov S.S. 2019. Brown bear (Ursus arctos) habitat suitability and distribution modelling in the southern taiga subzone using the method of maximum entropy. Nature Conservation Research 4(4): 34-64. DOI: 10.24189/ncr.2019.061 [In Russian] Peterson A.T., Soberon J., Pearson R.G., Anderson R.P., Martinez-Meyer E., Nakamura M., Araùjo M.B. 2011. Ecological Niches and Geographic Distributions. Princeton and Oxford: Princeton University Press. 329 p. Pettorelli N., Hilborn A., Broekhuis F., Durant S.M. 2009. Exploring habitat use by cheetahs using ecological niche factor analysis. Journal of Zoology 277(2): 141148. DOI: 10.1111/j.1469-7998.2008.00522.x Phillips S.J., Anderson R.P., Schapire R.E. 2006. Maximum entropy modeling of species geographic distributions. Ecological Modelling 190(3-4): 231-259. DOI: 10.1016/j.ecolmodel.2005.03.026 Posilico M., Meriggi A., Pagnin E., Lovari S., Russo L. 2004. A habitat model for brown bear conservation and land use planning in the central Apennines. Biological Conservation 118(2): 141-150. DOI: 10.1016/j. biocon.2003.07.017 Préau C., Trochet A., Bertrand R., Isselin-Nondedeu F. 2018. Modeling Potential distributions of three european

amphibian species comparing ENFA and MaxEnt. Herpetological Conservation and Biology 13(1): 91-104. Puzachenko Yu.G., Kozlov D.N. 2007. Geomorphological history of the development of the territory of the Central Forest State Nature Reserve. Proceedings of the Central Forest State Nature Reserve 4: 125-159. [In Russian]

Qi D., Hu Y., Gu X., Li M., Wei F. 2009. Ecological niche modeling of the sympatric giant and red pandas on a mountain-range scale. Biodiversity and Conservation 18(8): 2127-2141. DOI: 10.1007/s10531-009-9577-7 R Development Core Team. 2020. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Available from http://www.R-project.org Reutter B.A., Helfer V, Hirzel A.H., Vogel P. 2003. Modelling habitat-suitability using museum collections: an example with three sympatric Apodemus species from the Alps. Journal of Biogeography 30(4): 581-590. DOI: 10.1046/j.1365-2699.2003.00855.x Rinnan D.S., Lawler J. 2019. Climate-niche factor analysis: a spatial approach to quantifying species vulnerability to climate change. Ecography 42(9): 1494-1503. DOI: 10.1111/ecog.03937 Rinnan D.S. 2020. Package «CENFA», version 1.1.0. Climate

and Ecological Niche Factor Analysis. 36 p. Robin X., Turck N., Hainard A., Tiberti N., Lisacek F., Sanchez J.-C., Müller M., Siegert S., Doering M. 2020. Package «pROC», version 1.16.2. Display and Analyze ROC Curves. 95 p. Rotenbeny J.T., Preston K.L., Knick S.T. 2006. GIS-based niche modeling for mapping species' habitat. Ecology 87(6): 1458-1464. DOI: 10.1890/0012-9658(2006)87[1458:GN MFMS]2.0.C0;2 Santos X., Brito J.C., Sillero N., Pleguezuelos J.M., Llorente G., Fahd S., Parellada X. 2006. Inferring habitat-suitability areas with ecological modelling techniques and GIS: A contribution to assess the conservation status of Vipera latastei. Biological Conservation 130(3): 416-425. DOI: 10.1016/j. biocon.2006.01.003 Sattler T., Bontadina F., Hirzel A.H., Arlettaz R. 2007. Ecological niche modelling of two cryptic bat species calls for a reassessment of their conservation status. Journal of Applied Ecology 44(6): 1188-1199. DOI: 10.1111/j.1365-2664.2007.01328.x Silva L.D., Costa H., Azevedo E., Medeiros V., Alves M., Elias R.B., Silva L. 2017. Modelling native and invasive woody species: a comparison of ENFA and MaxEnt applied to the Azorean forest. In: A. Pinto, D. Zilberman (Eds.): Modeling, Dynamics, Optimization and Bioeconomics II. DGS 2014. Vol. 195. P. 415-444. DOI: 10.1007/978-3-319-55236-1_20 Strubbe D., Matthysen E. 2008. Predicting the potential distribution of invasive ring-necked parakeets Psittacula krameri in northern Belgium using an ecological niche modelling approach. Biological Invasions 11(3): 497513. DOI: 10.1007/s10530-008-9266-6

Thiebot J.-B., Lescroel A., Pinaud D., Trathan P.N., Bost C.-A. 2011. Larger foraging range but similar habitat selection in non-breeding versus breeding sub-Antarctic penguins. Antarctic Science 23(2): 117-126. DOI: 10.1017/S0954102010000957 Thioulouse J., Chessel D. 1992. A method for reciprocal scaling of species tolerance and sample diversity. Ecology 73(2): 670-680. DOI: 10.2307/1940773 Traill L.W., Bigalke R.C. 2007. A presence-only habitat suitability model for large grazing African ungulates and its utility for wildlife management. African Journal of Ecology 45(3): 347-354. DOI: 10.1111/j.1365-2028.2006.00717.x Tsoar A., Allouche O., Steinitz O., Rotem D., Kadmon R. 2007. A comparative evaluation of presence-only methods for modelling species distribution. Diversity and Distribution 13(4): 397-405. DOI: 10.1111/j.1472-4642.2007.00346.x Unger D.E., Fei S., Maehr D.S. 2008. Ecological niche factor analysis to determine habitat suitability of a recolonizing carnivore. In: Proceedings of the 6th Southern Forestry and Natural Resources GIS Conference. Athens, GA: University of Georgia. P. 237-250. Valle M., Boija A., Chust G., Galparsoro I., Garmendia J.M. 2011. Modelling suitable estuarine habitats for Zostera noltii, using Ecological Niche Factor Analysis and Bathymetric LiDAR. Estuarine Coastal and Shelf Science 94(2): 144-154. DOI: 10.1016/j. ecss.2011.05.031 Xuezhi W., Weihua X., Zhiyun O., Jianguo L., Yi X., Youping C., Lianjun Z., Junzhong H. 2008. Application of ecological-niche factor analysis in habitat assessment of giant pandas. Acta Ecologica Sinica 28(2): 821-828. DOI: 10.1016/S1872-2032(08)60030-X Zajec P., Zimmermann F., Roth H.U., Breitenmoser U. 2005. The return of the Brown bear to Switzerland - Suitable habitat distribution, corridors and potential conflicts. KORA Bericht 28: 1-31. Zaniewski A.E., Lehmann A., Overton J.McC. 2002. Predicting species spatial distributions using presence-only data: a case study of native New Zealand ferns. Ecological Modelling 157(2-3): 261-280. DOI: 10.1016/S0304-3800(02)00199-0 Zheltukhin A.S., Puzachenko Y.G., Kotlov I.P., Shironiya I.I., Sandlerskiy R.B. 2016. Spatial-temporal dynamics of the pine marten (Martes martes L.), mountain hare (Lepus timidus L.) and red squirrel (Sciurus vulgaris L.) trail activity in the European southern taiga. Zhurnal Obshchei Biologii 77(4): 262-283. [In Russian] Zimaroeva A.A., Zhukov A.V., Ponomarenko A.L., Matsyura A. 2015. Ecological niche modelling of Fringilla coelebs Linnaeus, 1758 (common chaffinch) using GIS tools. Romanian Journal of Biology -Zoology 60(2): 135-146. Zimmermann F. 2004. Conservation of the Eurasian lynx (Lynx lynx) in a Fragmented Landscape - Habitat Models, Dispersal and Potential Distribution. PhD Thesis. Lausanne: University of Lausanne. 193 p.

BROWN BEAR (URSUS ARCTOS) ECOLOGICAL NICHE AND HABITAT SUITABILITY MODELING IN THE SOUTHERN TAIGA SUBZONE

USING THE METHOD OF GNESFA

Sergey S. Ogurtsov

Central Forest State Nature Biosphere Reserve, Russia e-mail: [email protected]

The article presents the results of ecological niche and habitat suitability modelling of Ursus arctos, performed in the Central Forest State Nature Reserve and its buffer zone (Tver region, Russia), using the GNES-FA factor analysis method. Briefly discussed are the basic rules for performing such a study, approaches and modelling techniques. The ecological niche was modelled using FANTER, ENFA and MADIFA. The habitat suitability was modelled using Mahalanobis distances analysis. As ecogeographical variables, vegetation indices, morphometric characteristics of the relief, proximity rasters and types of landscape cover were used. The presence points of the species were recorded on permanent routes using a GPS navigator during 2008-2018. As a final result, six models with a different combination of input data (presence points and environmental variables) were selected. As a set of independent test data, we used U. arctos' presence points recorded by the forest rangers and research staff. The scenarios of the influence of utilisation weights on presence points and multicollinearity on variables are separately modelled. The global marginality of the best model was 1.17, while the tolerance was 0.69. The distance to anthropogenic food sources, the NDVI and the projective cover of grasslands made the highest contribution along the marginality axis. In the first three axes of specialisation, coniferous forests, Sphagnum pine (Pinus sylvestris) forests, and the distance to the rivers made the largest contribution. AUC values ranged from 0.66 to 0.73. Values of the continuous Boyce index were from 0.85 to 0.93. The area of the partially violated buffer zone turned out to be more suitable for U. arctos than the inner area of the Central Forest State Nature Reserve. Utilisation weights and the number of ecogeographical variables had a significant impact on the results of all analyses.

Key words: biotope, Brown bear, Carnivora, ENFA, GIS, HSM, MADIFA, Mahalanobis distances, spatial modelling, Ursidae