气候变化对中国梅花鹿潜在栖息地影响

刘 策 张 日 杜海荣 孙 悦 袁子奥 于晶晶 张明海

(东北林业大学野生动物与自然保护地学院，哈尔滨，150040)

由于全球的气候在不断发生变化，而这种变化又是影响全球生物的分布、生长、发育和繁殖的主要因素。并且已有研究表明，在未来几十年，气候变化还将成为损害生物多样性的主要原因。因此，在能够准确预测气候变化对未来生物发展影响的基础上从而提出有效的物种保护理论已经成为当前保护生物学研究的热点问题，也是保护生物学家面临的一个重大挑战[1-5]。

在评价气候变化对物种分布的影响时，物种具体分布信息是种群动态分析和种群未来发展趋势预测不可缺少的参数，能用于物种受胁程度评估和保护效绩评价[6]。而物种分布模型是一个很有效的工具[7]。物种分布模型的发展，始于BIOCLIM模型[8]的开发和应用，随后的20多年内，涌现了HABITAT、DOMAIN、生态位因子分析模型(ecological niche factor analysis，ENFA)、马氏距离(Mahalanobis distance，MD)、边界函数方法(border function，BF)、最大熵模型(maximum entropy，MaxEnt)、广义线性模型(generalized linear model，GLM)、广义加法模型(generalized additive model，GAM)、分类与回归树模型(classification and regression tree，CART)基于人工智能的模型[9]。

物种分布模型在预测种群分布动态及区域物种多样性变化等领域已得到广泛应用。如张薇等[10]应用MaxEnt模拟了驼鹿(Alcesalces)未来生境分布，并提出未来高排放情景下驼鹿潜在生境面积急剧下降，且潜在生境分布整体呈现向高海拔、高纬度迁移的趋势。罗翀等[11]应用ENFA和 MaxEnt预测林麝(Moschusberezovskii)在秦岭山系的生境分布。吴建国等[12]利用CART分析了在气候变化下大熊猫(Ailuropodamelanoleuca)的生境变化。

本研究选择全国地区作为研究区域，以中国梅花鹿(Cervusnippon)为研究对象。目前该物种被列为国家Ⅰ级重点保护动物，被《中国物种红色名录》列为极危物种(CR)。梅花鹿分布我国东北、四川、华南等地区，其栖息地变化对于梅花鹿的生物地理学、保护生物学研究极为重要。Aanes等[13]研究提出冬季降水量的增加使得鹿科(Cervidae)动物的种群数量减少。因此，本研究基于该物种当前分布现状，结合气候物种分布模型，预测分析了由气候改变引起的中国梅花鹿潜在生境分布的变化，从而能够更好地保护中国梅花鹿的生境。

1 研究区域及研究物种概况

研究区域包括我国的东北地区、华南地区和四川甘肃地区。其中东北地区包括黑龙江省、吉林省、辽宁省，地理位置为121°31′—134°44′E，38°43′—53°32′N，总面积约808 400 km2。东北地区的海拔变化幅度大，地形地貌复杂，包括大小兴安岭、松嫩平原、松辽平原、长白山和三江平原等自然区域。本区以温带季风气候和寒温带季风气候为主，水系密布。森林类型以针叶林居多，优势树种为落叶松(Larixgmelinii)和樟子松(Pinussylvestrisvar.mongolica)；其次多为针阔混交林，优势树种为红松(Pinuskoraiensis)、水曲柳(Fraxinusmandshurica)和蒙古栎(Quercusmongolica)[14]。本区为依靠森林生存的肉食及草食动物提供了广袤的栖息地。

华南地区包括江苏省、浙江省、江西省，地理位置为113°58′—122°56′E，24°29′—35°6′N，总面积约519 700 km2。梅花鹿主要在浙江清凉峰国家级自然保护区、江西桃红岭梅花鹿国家级自然保护区内，其分别具有明显的亚热带中山山地季风特征和温暖湿润的季风气候[15]。高海拔地区以针阔叶混交林为主，低海拔地区以常绿阔叶林和常绿落叶阔叶混交林为主。

四川甘肃地区，包括四川省以及甘肃省的陇南市、甘南藏族自治州、定西市、临夏回族自治州、兰州市、白银市、天水市，地理位置为97°24′—108°29′E，26°2′—37°35′N，总面积约519 700 km2。梅花鹿主要分布区域为四川省若尔盖县铁布梅花鹿自然保护区，该区域受西风环流及东南季风的影响，夏季温凉、冬季寒冷、干湿季明显，河岸及河漫滩地带以柳(Salixspp.)为主，其间分布有沙棘(Hippophaerhamnoides)、小檗(Berberisspp.)等的河滩灌丛；局部地带有油松(Pinustabuliformis)林、云杉(Piceaasperata)林、青杄(Piceawilsonii)林，该植被带是四川梅花鹿主要栖息生境[16]。

梅花鹿属于哺乳纲(Mammalia)，偶蹄目(Artiodactyla)，鹿科。梅花鹿是亚洲东部的特有种类，主要分布在俄罗斯东部、日本和中国。目前，梅花鹿的生存情况不容乐观，不仅在朝鲜和越南的种群面临着灭绝的风险，而且在韩国已经确定灭绝。

梅花鹿在历史上几乎遍布我国华东、华北、华中、华南、西南和青藏高原的东部[17]。由于气候环境的改变以及人为影响，至20世纪40年代，梅花鹿山西亚种(C.n.grassianus)、华北亚种(C.n.mandarinus)和台湾亚种(C.n.taiouanus)在野外灭绝。目前中国仅存3个亚种：东北亚种(C.n.hortulorum)、华南亚种(C.n.kopschi)和四川亚种(C.n.sichuanicus)，仅分布于黑龙江、吉林、江西、浙江、安徽、四川和甘肃等省狭窄的区域内[18]。

2 数据来源和研究方法

2.1 梅花鹿的空间分布数据

本研究搜集现有的文献及数据库数据，从中国知网、全球生物多样性信息网络(GBIF)数据库中收集梅花鹿物种分布数据进行筛选。本研究在开始前制定了规范的数据入选标准及处理方法：将在中国知网文献(学术期刊、硕博论文、会议、报纸)中有记载的我国野生梅花鹿的分布数据进行有效筛选并全部录入；将2000年后的全球生物多样性信息网络(GBIF)数据库中我国梅花鹿分布相关的数据全部录入；将梅花鹿分布数据中所记载的地点为城镇或林场、农场名称等面数据通过ArcGIS 10.3结合Google Earth提取这些面数据的中心经纬度坐标作为梅花鹿分布点，获得89个梅花鹿分布记录点。

最后将所获得的89个分布记录点与通过ArcGIS 10.3中根据梅花鹿的生境需求信息与1∶100万数字化植被图层进行叠加获得的梅花鹿适宜生境分布区图进行校对，去除不在梅花鹿适宜生境分布图中的分布记录点。

2.2 预测环境因子

本研究共选取与梅花鹿分布相关的23个预测环境因子(表1)。

生物气候变量：从世界气候数据库(WorldClim)版本2.1中下载(http：//www.worldclim.org)相关的生物气候变量19个，其包括降水、气温等相关气候因子，最热月/最冷月气温等(表1)。

空间分辨率为2.5 min(约等于5 km×5 km)[19-20]。选用的19个气候因子主要反映了温度与降水的特点及季节性变化特征，这些气候因子因具有较强的生物学意义已被广泛应用于物种栖息地的预测中。

当前气候条件下的生物气候变量为1970—2000年全球多个站点观测数据的平均值，未来生物气候变量来自北京气候中心参与CMIP6的中等分辨率的BCC-CSM2-MR模式输出数据，包括SSP2-4.5和SSP5-8.5两种气候情景下2021—2041年(2030s)、2041—2060年(2050s)、2061—2080年(2070s)3个时期的全球生物气候数据。

BCC-CSM2-MR采用IPCC第6次评估采纳的共享社会发展路径(SSPs)中的SSP1-2.6、SSP2-4.5、SSP3-7.0和SSP5-8.5情景，生物气候数据的空间分辨率为5 km×5 km。对BCC-CSM2-MR模式而言，其相对比CMIP5大气模式版本模式，虽然两者的大气物理过程方案相同，但BCC-CSM2-MR采用的大气分量模式对大气辐射、深对流过程和重力波方案进行了改进，同时考虑了气溶胶间接效应，使得模式能够合理地再现气候分布特征[21-22]。CMIP6中已开展的研究表明，BCC-CSM2-MR模式对全球气温和降水分布、大气辐射、海洋收支等方面的模拟相比CMIP5的早期版本有了显著改进[21]。

地形变量：DEM高程数据取自于中国地理空间数据云。根据DEM高程数据，使用ArcGIS 10.3以及R软件将DEM数据进行计算得坡度、坡向、海拔。

植被数据：植被覆盖指数(NDVI)，数据来源于中科院自然地理所提供的地表植被覆盖数据。

2.3 环境因子去相关性

利用ArcGIS 10.3软件将上述23个变量的栅格数据统一到相同坐标系、相同范围、数据精度统一重采样为5 km。各环境变量间存在相关性时，模型预测结果解释能力会相应较差。本研究对23个初始环境变量进行成对Pearson相关性检验。若相关系数|r|≥0.75，我们认为变量间存在相关性；再对其进行VIF相关性检验，若VIF相关系数大于4，则认为存在相关性，并保留VIF值较小的变量。最终，选择了8个环境变量用于梅花鹿潜在分布范围，分别为坡度、坡向、海拔、NDVI、等温性、最干季均温、最干月降水量、降水季节性。

2.4 环境预测模型

模型的模拟与评估选用基于最大熵理论的MaxEnt模型来模拟气候变化下梅花鹿潜在生境分布，该模型具有较高的预测精度与良好的稳定性，在物种分布模型中具有一定的竞争力。采用MaxEnt 3.4.0建模时，输入梅花鹿分布点数据和8个环境因子数据，设置参数，随机选择75%的物种分布点作为训练组，以获取建立最大熵模型的参数，25%作为测试组，以验证模型的准确性[23]。重复次数为10，最大迭代次数为500次，输出的数值为逻辑斯蒂值(logistic)，其他参数为软件默认参数。

2.5 模型评估

本研究选取受试者工作特征曲线(receiver operator characteristic curves，RO)下的面积(AUC)对模型精度进行评价。AUC值因其不受诊断阈值影响，且对物种发生率(prevalence)不敏感，目前是公认的最优的模型预测指标。AUC值评价模型的评估标准为：0.90—1.00，极好；0.80—0.90，好；0.70—0.80，一般；0.60—0.70，差；0.50—0.60，失败。模型模拟输出结果为物种存在概率栅格图，值在0—1内，值越接近1表示物种越可能存在。本研究采用灵敏度和特异度之和最大时对应的物种存在概率值P(P=0.32)作为阈值[24]。

2.6 模型优化

本研究采用空间筛除法和偏差栅格文件进行取样偏差的纠正。空间筛除法设定建模位点最小间距为10 km，即当2个分布点的距离小于10 km时，只取其中1点。通过上述方法对获得的梅花鹿分布记录点进行校对、筛选，最终确定89个梅花鹿分布点。偏差栅格文件采用所有89个位点绘制50 km缓冲区，在缓冲区内随机选取背景位点，将偏差栅格文件添加到MaxEnt软件高级设置中的bias file选项，软件根据取样点分布的栅格图层，生成与取样点布里格曼分歧(Bregman divergence)最小的背景位点进行建模[25]。

3 研究结果

3.1 梅花鹿当前潜在生境分布

输出结果中，东北亚种的潜在分布模拟结果的AUC值为0.992，华南亚种的潜在分布模拟结果的AUC值为0.938，四川亚种的潜在分布模拟结果的AUC值为0.933。其AUC值均大于0.9，模型模拟极准确，表明模型对梅花鹿潜在生境的预测效果较好。

当前气候条件下，为了充分研究潜在生境的适宜度，将梅花鹿适宜生境分为最适宜生境(P>0.5)和中适宜生境(0.32

梅花鹿华南亚种最适宜生境面积为18 208.41 km2，占研究总面积3.504%，中适宜生境的面积为60 286.3 km2，占研究总面积的和11.6%。其适宜生境为安徽南部、江西东北部、浙江西北部，这与现实已观测到的结果相符。其中最适宜生境集中分布在浙江省临安区、淳安县；安徽省歙县、绩溪县、黄山区、黟县、石台县、江西省彭泽县、浮梁县等。

梅花鹿四川亚种最适宜生境面积为48 610.50 km2，占研究总面积7.721%，中适宜生境的面积为108 850.80 km2，占研究总面积的和17.29%。分布于四川省若尔盖县的铁布、包座自然保护区以及与之相邻的甘肃省迭部县和四川省九寨沟县的部分区域，这与现实已观测到的结果相符。其中最适宜生境集中分布在四川九寨沟、铁布保护区、甘肃省迭部县、卓尼县等部分地区。

3.2 影响梅花鹿潜在生境分布的重要因子分析

对23个环境变量去相关性后得到的8个环境因子：坡度、坡向、海拔、NDVI、等温性、最干季均温、最干月降水量、降水季节性，对模型的贡献率是不相等的(表2)，东北亚种适宜生境预测中，对模型贡献率较大的环境变量为降水季节性(49.3%)，降水相关变量对模型的贡献率54.6%；温度相关变量对模型的累计贡献率为18%；其他变量：高程、坡度、坡向和NDVI的贡献率依次是19.5%、0.3%、4.6%和3.1%。华南亚种适宜生境预测中，对模型贡献率较大的环境变量为高程(37.9%)，温度相关变量对模型的累计贡献率为38.5%；降水相关变量对模型的贡献率14.1%；其他变量：坡度、坡向和NDVI的贡献率依次是7.6%、1.6%和0.3%。

表2 环境变量贡献率

四川亚种适宜生境预测中，对模型贡献率较大的环境变量为最干月降水量(61.1%)。降水相关变量对模型的贡献率79.8%；温度相关变量对模型的累计贡献率为8.5%；其他变量：高程、坡度、坡向和NDVI的贡献率依次为1.5%、2.6%、4.2%和3.4%。

梅花鹿东北亚种生活在高纬度地区，温度相关的环境因子对模型的贡献率远不如降水相关的环境因子，四川亚种也是降水相关的环境因子大于温度相关的环境因子对模型的贡献率，而华南亚种中，海拔因子对模型的贡献率最大，温度相关变量其次。由此可见，即使是同种不同亚种的物种，环境因子和气候因子对其适宜生境预测模型的贡献率也大不相同。

3.3 气候变化下梅花鹿分布区面积的变化

在气候变化影响下，当前梅花鹿适宜生境面积的减少和新梅花鹿适宜生境面积的增加引起了梅花鹿总适宜生境面积变化。利用ArcGIS 10.3对未来梅花鹿适宜生境(P>0.32)变化的面积和比例进行统计(表3)。表3显示，气候变化对梅花鹿的生境分布面积产生极大的影响。

表3 不同气候变化情景下梅花鹿不同亚种适宜生境面积变化

随着时间的推移，两种气候变化情景下研究区内梅花鹿当前潜在生境面积变化明显，梅花鹿东北亚种在ssp245情景下，随着时间推移，适宜生境面积在2030s、2050s和2070s阶段分别增加92.33%、138.15%和265.46%，适宜生境逐渐扩大，而在ssp585情景下，梅花鹿适宜生境面积在2030s、2050s阶段增加幅度将大于ssp245情景，其在相应阶段分别增加173.14%和524.79%，而在2070s阶段，面积增幅减小，增加了139.12%，减少了71.56%，分布区总体向北迁移。随着气候变化，未来梅花鹿东北亚种在我国适宜生境面积将不断增加，且在主要分布区如老爷岭、珲春等地分布有国家重点保护动物东北虎(Pantheratigrisaltaica)、东北豹(Pantherapardusorientalis)等大型食肉动物，对虎豹数量恢复极其重要。

梅花鹿华南亚种在ssp245情景下，随着时间推移，适宜生境面积在2030s、2050s和2070s阶段分别减少78.77%、86.89%和88.75%，而在相应阶段适宜生境增加面积最大仅为10%，生境破碎化严重，适宜生境仅为原有十分之一。而在ssp585情景下，梅花鹿适宜生境面积减小幅度将大于ssp245情景，其在相应阶段分别减少77.61%、96.29%和93.62%。适宜生境逐渐消失。至2070s阶段，仅在浙江清凉峰自然保护区千顷塘区域有高适宜生境。梅花鹿四川亚种鹿模型预测结果表明随着气候变化。

梅花鹿四川亚种在ssp245情景下，随着时间推移，适宜生境面积在2030s、2050s和2070s阶段分别减少20.03%、20.44%和15.49%，减少面积相对稳定，并且都有5%—10%的新适宜生境，而在ssp585情景下，梅花鹿适宜生境面积减小幅度将大于ssp245情景，在相应阶段分别减少29.58%、43.34%和40.23%。而新适宜生境的增加仅为5%—6%，面积在随时间变化而不断减少。未来梅花鹿四川亚种在我国适宜生境主要位于川西高原的北部边缘，光盖山，隶属于四川省阿坝藏族羌族自治州若尔盖县，涉及冻列、崇尔、热尔3乡。

在ArcGIS 10.3中，将未来梅花鹿潜在生境分布图处理，获得梅花鹿潜在生境分布区随时间变化的地理位置(图2，图3，图4)。由此可知，未来梅花鹿生境面积将急剧变化，东北亚种适宜生境面积在两种气候模式下，随时间变化不断增加，栖息地不断扩大，但在ssp585情景2070s阶段，适宜生境呈现破碎化，并向北迁移。华南亚种适宜生境将主要分布在清凉峰一带，并且适宜生境逐渐缩小、破碎，在ssp585情境下，在2050s阶段，梅花鹿的高适宜生境面积仅为当前高适宜生境的4.99%。四川亚种，随着时间的推移，适宜栖息环境变化不大，在ssp585情景下比ssp245情景面积减小明显。

4 讨论

4.1 梅花鹿分布与环境因子关系

MaxEnt模型预测分析结果表明：梅花鹿潜在适宜生境的重要环境因子为降水和海拔因子。影响东北亚种和四川亚种潜在适宜生境分布的主要是降水因子，贡献率最大的因子分别为降水季节性(bio15)和最干月降水量(bio14)，东北地区夏季降水百分率都在60%以上，超过全国平均水平，而冬季降水百分率却小于4%[26]。四川亚种分布于川西北高山高原高寒气候区。该区海拔高差大，气候立体变化明显，从河谷到山脊依次出现亚带、暖温带、中温带、寒温带、亚寒带、寒带和水冻带。总体上以寒温带气候为主，水热不足，干湿季明显，已有研究表明，鹿科动物更愿意生活在水资源丰富的区域，因此，降水变化对这两个亚种的适宜生境预测贡献最大。

对华南亚种潜在适宜生境分布影响最大的因子为海拔，温度和降水因子也起到50%以上的作用。但因其具有种群数量稀少、分布区域面积狭窄以及其栖息环境的复杂性和特殊性等特点，对环境因子的要求大于对气候因子的要求。

4.2 模拟结果准确性评价

目前，模型模拟预测法已成功应用到气候变化下动、植物的分布预测中[27]。AUC值对模型结果评价表明，MaxEnt模型预测的梅花鹿生境分布范围具有极高的准确性。然而梅花鹿的潜在生境分布由多种因素综合影响决定，除本研究考虑的气候和地形因素外，生物间相互作用(捕食、竞争等)、物种的生物学特点(物种生活史、迁移能力等)和人类活动对梅花鹿生境分布的影响也不能忽略[28-29]，若能充分考虑以上因素将其量化并加以计算，则模型预测结果将更加接近物种的现实分布。

4.3 保护建议

4.3.1 对保护区范围进行调整

我们建议主要在区内新建自然保护区，填补地区的保护空缺。同时可扩展原有国家级自然保护区的面积，如黄泥河、清凉峰和桃红岭自然保护区等。并通过提升省、市、县级自然保护区面积来使自然保护区面积满足梅花鹿的需求。

4.3.2 减少人类活动的影响

进一步加大保护梅花鹿宣传，丰富宣传手段，充分利用网络电视等公共信息获取平台。提高人们保护野生动物的意识，让社会共同关注，参与到保护野生动物事业中来，才能引起社会各阶层对梅花鹿的保护认识。

4.3.3 封山保护，扩大种群

随着保护区和生态公益林的建立，大规模地砍伐森林、拓垦土地、家畜放牧等行为得到了有效控制，生态公益林和保护区内植物长势良好，早期破坏过的植被也已恢复，但植物生长茂盛后形成的常绿落叶阔叶林并不适宜梅花鹿栖息。因此针对以上问题，应当适当在生境区内构建适宜梅花鹿生存的乔、灌、草相结合的三层复合林分，进一步改造适宜生境，扩大种群数量[30]。

4.3.4 适时投食

梅花鹿作为食草性动物，冬天大雪封山会导致梅花鹿因觅食等原因远离原先的活动区域而提高被猎杀的风险，因此，在冬季大雪封山之时，应该在梅花鹿活动区域及时投放供给一定的青草和干草，以便缓解梅花鹿食物的不足。

4.3.5 加大科研投入

为了有效保护梅花鹿这个物种，建立长期、有效的监测计划来对梅花鹿进行动态的、全方位的监测是必须的。政府部门应加大对梅花鹿研究项目的补助，鼓励科研机构开展梅花鹿的科学研究，推进梅花鹿的种群扩繁和种质资源保存。另外，确认优先保护区域有利于有效管理[31]。

4.4 国外保护策略

欧共体在1993年12月批准了生物多样性公约。欧盟生物多样性战略规划出了一个明确的行动框架，并规定了总目标和部门目标，列出了4项计划：保护和持续利用生物多样性，分享利用遗传资源所产生的利益，研究、鉴定、监测和信息交换，教育、培训和提高公众的保护意识。

澳大利亚在1999年通过了环境与生物多样性保护法(简称EPBC法)。根据该法，澳大利亚联邦政府(由环境与遗产部负责)、州及地方政府、非政府机构、慈善机构以及社区团体共同工作，对本国原产的动植物种进行保护。

在进行城市规划或重大基础设施建设前期，将是否会对范围内的野生动物造成影响纳入前期生态分析。例如意大利政府在规划建设跨海大桥时，委托野生动物保护研究所进行了“跨海大桥对鸟类迁徙的影响”的论证，从而判断该桥是否可以建设。

5 结论

本文基于最大熵模型和物种分布数据，研究了气候变化下我国梅花鹿潜在分布及其变化趋势，分析了保护建议预测。结果表明：当前气候条件下，预测的适宜生境都大于现阶段梅花鹿生境分布区域，在未来气候变化下，东北亚种高适宜生境分布范围将逐渐增加，四川亚种高适宜生境分布范围减少区域不大，主要分布区无变化，华南亚种高适宜生境分布范围将在ssp585情景2050s时趋向于消失，破碎化极其严重。当前及未来气候条件下，梅花鹿的优先保护区域分布变化较小，东北亚种主要分布于黑龙江东南部和吉林东部临近边境地区，华南亚种安徽南部、江西东北部、浙江西北部，四川亚种适宜生境集中分布在四川九寨沟、铁布岭保护区、甘肃省迭部县、卓尼县等部分地区。但保护空缺面积较大。因此，应加强梅花鹿适宜地区保护地体系建设，这将有益于当前及未来梅花鹿和本区域其他物种的长期保护。