穿山甲是哺乳动物中一个较小的类群, 隶属于哺乳纲(Mammalia)、鳞甲目(Pholidota)、鲮鲤科(Manidae)、穿山甲属(Manis)^[1]。目前, 全世界现有8种穿山甲：中华穿山甲(Manis pentadactyla)、印度穿山甲(M. crassicaudata)、马来穿山甲(M. javanica)、大穿山甲(M. gigantea)、树穿山甲(M. tricuspis)、南非穿山甲(M. temminckii)、长尾穿山甲(M. teyradactyla)和巴拉望穿山甲(M. culionensis)。我国分布的穿山甲主要是中华穿山甲^[2-4]。该物种2014年被IUCN红色名录(Red List)列为“极危”(CR)^[5]；2017年被列入CITES附录Ⅰ(www.cites.org)；2020年6月5日, 我国将其从国家二级保护动物提升为国家一级保护动物。在过去几十年间, 滥捕滥猎和非法贸易使中华穿山甲种群数量急速下降；据估计, 其国内种群数量已减少90%以上^[6-7]。尽管中华穿山甲的濒危程度不断提高, 但近年来我国一直没有开展过全国性的野生中华穿山甲专项调查, 也未见相关报道, 其种群现状尚不明显, 分布位置尚不清晰, 在很大程度上限制了对该物种的有针对性的保护管理。因此, 基于已知中华穿山甲的分布位点预测其潜在适宜区对于其种群保护和栖息地管理具有十分重要的意义。

掌握濒危物种的种群状况和生态需求是开展物种保护的基础^[8]。掌握目标物种的空间分布和生境需求是濒危物种保护的一项重要工作, 但是有限的野外数据有时会限制对濒危物种分布区的掌握。生态位模型通过目标物种的地理分布数据和物种所在地的环境数据, 在物种已知的地理分布区域内进行测量或估算, 进而预测该区域内哪些区域满足该物种生态位的要求, 从而成为该物种潜在分布区的一部分^[9]。主流的生态位模型有GARP、ENFA、DOMAIN、BIOCLIM和MaxEnt。其中, 最大熵模型(the maximum entropy model software, MaxEnt)预测的结果精确度更高^[10]；即使在目标物种分布点比较少的情况下, MaxEnt模型也可以获得较准确的预测结果^[11-14]。自Phillips提出该模型以来, MaxEnt在外来物种入侵^[15-17]、动植物保护^[18-20]和全球气候变化对物种分布的影响等方面得到广泛应用^[21-23]。

在此背景下, 本研究对近年来我国各省市有关中华穿山甲的救护记录、救助新闻进行了整理, 收集了67个中华穿山甲分布数据, 选取地形地貌、地表类型、土地覆被、气候变量、生物学变量和人类影响变量等25个全球尺度的环境变量, 结合地理信息系统(GIS)、统计分析软件SPSS和最大熵模型等软件预测了中华穿山甲在我国的潜在适宜分布区, 以期能够为全国中华穿山甲种群专项调查和种群保护管理提供科技支撑。

1 材料与方法 1.1 数据收集与处理

利用MaxEnt模型预测物种适宜分布区需要物种的分布数据和环境数据。中华穿山甲的分布数据主要来源于互联网信息与各地的救助记录。我们通过百度搜索引擎查找关键词“中华穿山甲+省份”搜寻近十年来我国关于中华穿山甲的救护新闻, 筛选出有关中华穿山甲的救护报道, 提取位置信息, 借助百度地图输入分布点地理位置信息, 确定其十进制格式的经纬度。将上述途径所得到的分布点进行仔细筛选、整理, 获取分布点共67个(表 1)。

使用ArcGIS软件邻域分析中的缓冲区模块对分布数据点进行剔除, 保证每个栅格仅有一个分布点, 避免人为因素导致分布点过多而出现过度拟合现象。本研究所用环境数据空间分辨率为2.5 arc-minutes(约4.5 km), 设置缓冲直径为3 km, 当两个分布点之间的距离小于3 km时, 缓冲区面积会重叠, 此时只保留其中一个分布点^[24], 最后保留有效分布点65个(图 1)。分布数据按species、longitude和latitude格式依次输入Excel单元格, 并最终保存为MaxEnt模型支持的CSV格式文件用于后续分析。

本研究选取的环境数据包括1.地形地貌变量：海拔、坡度和坡向；2.地表类型变量：土地覆被类型；3.气候变量：19个气温和降水因子, 如年平均气温、最冷月最低温度和最冷季度降水量等；4.生物学变量：归一化植被指数(NDVI)；5.人类影响变量：人口密度。这些数据来源及生物学意义详见表 2。使用ArcGIS软件将上述25个全球尺度的环境变量掩膜提取出中国区域的环境数据, 并重采样统一精度2.5 arc-minutes, 通过栅格转换工具将栅格数据转换为*.ASCII格式数据以进行MaxEnt模型运算。

19个气候变量是以温度和降水量为基础数据计算得出, 因此它们之间不可避免地存在自相关和多重共线性等问题, 在模型运行时带入冗余信息, 导致输出的AUC值偏高^[24]。为此, 我们对所选的25个环境变量进行了主成分分析, 以剔除贡献率低的环境变量, 仅保留贡献率大于1.0%的环境变量^[25]。同时, 使用SPSS软件对所有环境变量进行Spearman相关性分析(附录), 如果两个环境变量之间的相关系数绝对值大于0.8, 说明两个变量间存在显著相关性, 剔除两者中对物种影响较小的环境变量^[26-28]。最终, 我们选择了6个环境变量对中华穿山甲进行适宜分布区预测, 包括最冷季度降水量(Bio_19)、人口密度(PD)、年降水量(Bio_12)、坡度(Slope)、坡面(Aspect)、海拔(Alt)。

1.2 模型运算与评价

将处理好的分布数据与环境数据导入MaxEnt模型软件, 创建响应曲线并设置75%的分布点作为训练集(training data), 25%为随机测试集(Random test percentage)；勾选“Random seed”选项, 重复次数(Replicates)设为10次^[29-30]；高级设置中勾选“Write plot data”选项, 方便后期制图与数据统计, 其他选项保持默认, 进行模型运算。根据最大训练敏感度和特异度法^[31-33] (Maximum Training Sensitivity Plus Specificity, MaxSS)确定阈值为0.166, 然后通过重分类, 将物种适宜区划分为4个等级：适生值＜0.166为非适宜区、0.166—0.332为低适宜区、0.332—0.498为中适宜区、适生值＞0.498为高适宜区, 接着计算4个适宜区在中国的分布面积。

采用刀切法(Do jackknife)来评价各变量对模型的重要性, 使用受试者工作特征曲线(ROC, Receiver Operating Characteristic Curve)与横坐标围成的面积即(Area Under Curve, AUC)来评价模型预测结果的准确性。AUC值在0—1之间：0.5—0.6表现为失败, 0.6—0.7表现为较差, 0.7—0.8表现为一般, 0.8—0.9表现为好, 0.9—1表现为极好；值越接近1表示模型预测准确度越高。

2 结果与分析 2.1 主导环境因子贡献率

通过MaxEnt模型分析了6种环境变量对中华穿山甲分布的贡献率, 结果表明气候变量中的降水因子是限制中华穿山甲分布的主导因素, 累计贡献率达85.4%；其中, 最冷季度降水量单因子贡献率高达68.9%。人类影响因子和海拔因子的累计贡献率达14.2%, 坡面和坡度的贡献率相对较小, 两者累计仅为0.4%(图 2)。

2.2 刀切法分析单一环境变量的重要程度

我们用刀切法进一步分析了单一环境变量对中华穿山甲分布结果的影响率(图 3)。蓝色表示当仅用此环境变量用于MaxEnt模型建模时, 该变量对模型的增益效果；绿色表示除此变量, 其他变量对模型的总体增益效果；红色表示所有环境变量应用模型所产生的增益效果。结果表明最冷季度降水量和年降水量对中华穿山甲的分布起主导作用, 训练增益分别为1.85和1.83；人口密度和海拔是影响中华穿山甲分布的次要因素, 训练增益分别为0.86和0.71；影响最小的是坡面和坡度, 训练增益不超过0.10。这与相对贡献率分析的结果是一致的, 表明所选的环境变量能够较好地反映中华穿山甲真实的生活环境^[25]。

2.3 主要环境变量的响应分析

中华穿山甲的分布与主要的环境变量之间具有较强的响应关系(图 4), 当存在概率大于0.498时, 可以认为此环境变量值最适合中华穿山甲生存。其中, 最冷季度降水量、年降水量、人口密度、海拔四个主要环境变量对中华穿山甲分布的最适宜范围分别为：最冷季度降水量141.22—439.46 mm、年降水量1471.67—2386.56 mm、人口密度390—28894人/km²、海拔＜316.98 m。

最冷季度降水量小于29.80 mm时, 中华穿山甲的存在概率几乎为零；小于93.85 mm时, 中华穿山甲的存在概率小于0.166, 为非适宜区；93.85—120.54 mm时, 存在概率大于0.166且小于0.332, 为低适宜区；120.54—141.22 mm或大于439.46 mm时, 存在概率大于0.332且小于0.498, 为中适宜区；141.22—439.46 mm时, 存在概率大于0.498, 为高适宜区。

年降水量小于1218.53 mm时, 为非适宜区；1218.53—1365.09 mm或大于3234.83 mm时, 为低适宜区；1365.09—1471.67 mm或2386.56—3234.83 mm时, 为中适宜区；1471.67—2386.56 mm时, 为高适宜区, 且在1560.50—1675.97 mm之间时, 中华穿山甲的存在概率达到峰值。

人口密度在30—6795人/km²时, 存在概率随人口密度的增加而增加, 在6795—36742人/km²时, 存在概率随人口密度增加而减少。低适宜区范围：0—30人/km²；中适宜区范围：30—390人/km²或28894—32457人/km²；高适宜区范围：390—28894人/km²。

海拔范围大于3351.13 m时, 中华穿山甲存在概率为零。随着海拔的降低, 中华穿山甲的分布概率增加。海拔范围在602.92—1166.86 m时, 为低适宜区；在324.92—563.20 m, 为中适宜区；海拔＜316.98 m, 为高适宜区。

2.4 中华穿山甲在中国的适宜区预测

MaxEnt模型运行十次后, 筛选AUC值最高的模型作为最优模型。按照重分类标准将中华穿山甲适宜性指数重分类成四类(图 5)。结果表明, 中华穿山甲适宜分布区总面积约为74.27×10⁴ km², 占中国国土面积的7.73%, 其中低、中、高适宜区面积分别为32.53×10⁴ km²、26.86×10⁴ km²和14.88×10⁴ km², 占总适宜区面积的43.80%、36.16%和20.04%。中华穿山甲在江西省和湖南省的高适宜区面积最大, 分别为6.87×10⁴ km²和3.13×10⁴ km², 占该省面积的44.96%和16.14%；其次是广东省和广西省, 面积分别为2.96×10⁴ km²和2.07×10⁴ km², 占该省面积的19.15%和9.91%, 具体各适宜区在各省(市、自治区)所占的面积及比例见表 3。

2.5 ROC曲线和AUC值对模型预测的准确性检验

基于所有环境变量与主导环境变量构建MaxEnt模型获得ROC曲线和AUC值(图 6), 获得的训练集AUC值为别为0.982和0.961(误差范围±0.014), 模型预测结果达到极好标准, 表明该模型对中华穿山甲在中国的潜在适宜区的预测结果非常好, 可信度比较高。

3 讨论

为了更有效地保护濒危野生动物, 首先要明确的是该物种的分布情况, 还要确定该物种的生境需求与所处环境之间的关系, 从而为该物种提供最合理、最有效、最具体的保护方案^[34]。本研究基于近年来我国中华穿山甲救护记录, 利用MaxEnt模型对其潜在适宜分布区进行了预测。研究结果揭示了中华穿山甲在中国的适宜分布区主要在我国长江以南地区, 主要集中在江西、广东、湖南和浙江省, 预测值与观测值高度吻合, 预测结果与中国兽类野外手册^[35]比较吻合, 预测结果比较准确, 预测结果可为我国中华穿山甲的种群保护和栖息地管理提供有力的科技支撑。

分布数据和环境数据的过度拟合通常会导致模型输出的AUC值偏高(图 6)。为了最大限度地避免数据点的偏差和模型的过度拟合, 我们对已获取的分布数据点进行了筛选和剔除, 保证每个栅格里只保留一个数据分布点。同时, 对所有的环境变量进行了相关性分析来确保选取的环境数据与中华穿山甲的分布密切相关, 以此减少冗余信息对模型结果的影响, 使模型输出后的AUC值更加真实。MaxEnt模型输出的AUC平均值0.961, 表明该模型的预测结果具有较高的可信度。

我们也分析了影响中华穿山甲空间分布的关键环境因子。结果表明, 气候变量中的降水因子是限制中华穿山甲分布的主要因素, 原因可能是降水因子与中华穿山甲的食物——白蚁的地理分布密切相关。1只中华穿山甲每年能吃700万只蚂蚁和白蚁, 这些蚁类在我国的活动分布主要在淮河以南的广大地区(热带、亚热带等温暖多湿地区), 向北逐渐稀少, 向南逐渐递增。在这些地区, 白蚁的分布与雨量关系极为密切, 尤其是江西及江西周围省份温度适宜、雨量充沛且白蚁种类繁多, 这为中华穿山甲提供了很好的食物来源。同时, 白蚁的分布区也与本研究预测的中华穿山甲潜在适宜区高度吻合^[36]。

物种与环境的关系是研究物种生态学需求和空间分布的重要方面。本研究认为, 虽然中华穿山甲的栖息地和食物来源会受到气候变化的影响, 但物种变化的过程是一个极其缓慢的过程^[37-38], 物种的生态位压缩主要还是由人类影响所驱动。栾晓峰等研究发现短时间内的气候变化对于物种分布呈现中性的效果, 或在某种情况下甚至可能是积极的, 人类活动干扰是影响中华穿山甲分布变化的最大因素^[39]。国家林业局1998年前后开展的中国重点陆生野生动物资源调查结果表明, 我国野生穿山甲的数量约为64000只^[3]。到2008年, 张立等人认为中国野生穿山甲的数量在25100—49450只之间, 即10年内其种群就减少了约1/2的数量^[3]。虽然中华穿山甲的濒危等级不断提高, 但近十年来我国一直没有开展全国性的野生穿山甲专项调查, 也未见相关报道, 在很大程度上限制了对该物种有效保护。穿山甲是地球上最古老的物种之一, 在地球上至少已延续了5000年。因其进化程度低、对新的环境适应能力差、繁殖力低下、人工难于驯养等原因, 一旦遭到大量捕杀导致种群数量下降后就难以恢复^[6]。如若不及时开展对中华穿山甲的生物学研究, 潜在的灭绝危险将进一步增加。

MaxEnt模型在物种分布数据很少的情况下预测结果也比较准确, 但是不充分的分布信息也会导致模型预测产生偏差^[40]。在避免过拟合的情况下, 物种的分布数据越多越好。因此, 基于更多野外活动点的中华穿山甲分布区研究仍然有待深入。本研究所获得的部分物种分布数据是通过百度地图提取其救助地点所获取的, 少数新闻救助点的具体位置只能精确到镇或乡。因此, 该救助地点的环境状况无法完全反应物种实际分布点的环境状况, 导致结果会有一定偏差。譬如，救护地点大多接近人类聚集区，周边人口密度较高，可能会导致模型分析过高估计穿山甲适宜生境中的人口密度；本研究中MaxEnt模型预测穿山甲适宜生境的人口密度为390—28894人/km²，虽然在一定程度上反映出穿山甲适宜生境距离人类活动区较近，但可能过高地估计了穿山甲适宜生境中的人口密度。另外, 虽然本研究选取的土地覆被类型和植被归一化指数有助于精确模拟物种适宜生长的微环境, 但是两者对于模型的贡献度过小, 使得我们无法从区域或微生境尺度分析环境因子对其分布的作用关系, 只能从较大的空间尺度预测其潜在分布区。因此, 多尺度和高精度上的中华穿山甲生境利用和空间分布适宜性预测有待进一步开展。