基于组合模型预测蛇足石杉在中国的适生区

陈 勤，孔梦盈，孙奇芳，管毕财

(南昌大学 生命科学学院，南昌 330031)

气候是物种分布的决定因素之一[1]。气候变化会改变栖息环境，降低生物多样性，增加濒危植物灭绝的风险[2]。21世纪以来全球气候变化复杂多样，物种对气候的适宜性随之发生变化，从而影响物种分布格局，进而威胁物种的生存[3]，迫切需求将气候变化纳入物种保护的研究中[4]。

物种分布模型(species distribution models, SDMs)是基于生态位理论量化环境与物种之间关系的相关性模型[5]。它能够揭示物种对栖息地的适宜性与环境的关系，被广泛用于生态学研究中，特别是在预测物种的分布及其热点[6]、管理入侵物种[7]、保护濒危物种[8]及验证气候变化与物种分布的关系[9]等方面。用于构建物种分布模型的算法主要有：人工神经网络(artificial neural networks, ANN[10])、生态位因子分析(ecological niche factor analysis, ENFA[11])、广义线性模型(generalized linear models, GLM[12])及最大熵模型(maximum entropy, MaxEnt[13])等，但是基于单一算法构建的模型准确性有待提升[14]。基于R程序包Biomod2[15]构建组合模型，能够提升预测结果的可靠性，日益广泛应用于模型构建。Ren等[16]使用Biomod2调查了对中国云南亚热带常绿阔叶林的物种分布及生物多样性。Ray等[17]使用Biomod2预测了橡树林地的分布。Vacchiano等[18]使用Biomod2预测了意大利西北部入侵植物窄叶黄菀的入侵情况与潜在分布。Legal等[19]使用Biomod2研究了分布有限且不均匀的孑遗生物短角宝凤蝶的潜在分布及其保护，这些研究均得到了比较好的预测效果。

蛇足石杉[Huperziaserrata(Thunb. ex Murray) Trevis.]属于石松科(Lycopodiaceae)石杉属(Huperzia)多年生蕨类植物[20]。蛇足石杉喜阴湿[21]，主要生长在海拔350～1 700 m的林下、灌丛和路旁[22]，广泛分布于亚洲、大洋洲、美洲热带、亚热带及温带地区，在中国大部分地区都有分布但生物量普遍较少[23]。蛇足石杉具有重要的药用价值，从蛇足石杉中分离出的蛇足石杉碱及石杉碱甲等生物活性物质可用于治疗阿兹海默症、重症肌无力和记忆力下降等临床症状[24-26]。但目前蛇足石杉产业化繁育尚未取得成功，材料获取主要依赖野生资源，由于过度采集和生态环境恶化等因素影响，生物资源骤减[27]。因此，对蛇足石杉生物资源的保护刻不容缓，2021年9月，蛇足石杉被国家林草局正式列为国家二级重点保护野生植物[28]。迄今为止，对蛇足石杉的研究也主要集中在化学成分[29]、药理作用[30]、快速繁殖[31]及分子生物学[32]等方面，为了开展对该物种的保护工作，我们需要调查收集蛇足石杉的分布情况并进行适生区的预测。因此，本研究基于蛇足石杉在中国的分布数据和气候数据构建组合模型，预测蛇足石杉当前和未来气候情景下在中国的适生区并计算适生区面积及质心变化，为蛇足石杉及其生物资源的保护提供理论支持。

1 材料和方法

1.1 蛇足石杉分布数据

通过中国数字植物标本馆(http://www.cvh.ac.cn/)、中国国家标本资源平台(http://www.nsii.org.cn/)、全球生物多样性信息网络(https://www.gbif.org/)和文献查阅获取蛇足石杉分布数据，对有准确地址和经纬度的数据进行校对，对有准确地址无经纬度的数据进行经纬度定位，剔除重复以及过于密集的分布数据，最终总共得到186条蛇足石杉分布数据，从地理国情监测云平台(http://www.dsac.cn/)下载国家行政区矢量地图作为底图进行可视化表达(图1)。

图1 蛇足石杉在中国的分布

1.2 气候数据的获取与筛选

从世界气候数据库(WorldClim, https://www.worldclim.org/)下载20世纪中、后期(1960～1990)的气候数据作为当前气候数据和对中国气候模拟能力较好的CCSM4模式[33]下的21世纪下半叶(2070)的RCP2.6、RCP4.5、RCP6.0及RCP8.5这4种典型浓度路径(representative concentration pathways, RCPs)的气候数据作为未来气候数据，气候数据分辨率均为2.5 minutes。RCPs表示2100年相对于1750年的辐射强迫总和，向大气排放的温室气体和气溶胶会影响地球辐射平衡产生辐射强迫[34]。RCP2.6辐射强迫为2.6 W·m-2，排放处于极低水平，到2100年地球表面温度将上升0.3～1.7 ℃；RCP8.5为8.5 W·m-2，排放处于高水平，到2100年地球表面温度将上升2.6～4.8 ℃；RCP4.5和RCP6.0分别为4.5 W·m-2和6.0 W·m-2，排放水平介于二者之间，到2100年地球表面温度将分别上升1.1～2.6 ℃和1.4～3.1 ℃[35]。

气候数据包括19个生物气候变量(表1)，使用ArcGIS 10.6提取中国生物气候变量数据，再用SDM Toolbox[36]移除相关性大于0.7的生物气候变量，避免生物气候变量之间相关性过高影响适生区预测准确性，将等温性(Bio3)、温度季节性(Bio4)、最暖月最高温(Bio5)、最冷季均温(Bio11)、最湿月降水量(Bio13)和降水季节性(Bio15)6个变量用于适生区预测。

表1 19个生物气候变量

1.3 物种分布模型的构建与评估

基于Biomod2构建物种分布模型需要物种分布数据和不分布数据，不分布数据很难获得[11]，通常生成假缺失数据作为不分布数据[5]。本研究通过Biomod2随机生成1 000个假缺失数据作为不分布数据，使用Biomod2包含的10种算法(表2)分别构建物种分布模型。采用接收者操作特征曲线(receiver operating characteristic curve, ROC)和真实技巧统计法(true skill statistic, TSS)对模型进行评估。AUC(area under the curve)是ROC与横坐标围成的面积，值介于0.5～1之间[37]。TSS不依赖发生率，能很好评价预测结果[38]。使用TSS值大于0.8的算法用于构建物种分布组合模型，将AUC值作为组合模型的最终评价标准，其值越大说明预测效果越好。在模型构建过程中使用Biomod2内置评价标准对用于构建物种分布模型的6个生物气候变量进行重要性评估。将模型预测结果输入ERDAS IMAGINE 2015导出为TIFF格式，再使用ArcGIS 10.6对TIFF文件进行地理配准及可视化表达。

表2 10种用于模拟的算法

1.4 物种适生区及其变化预测

使用ArcGIS 10.6对组合模型预测结果进行归一化，将预测结果转化为0-1的形式，采用自然断点分级法(natural breaks)进行分类，将预测结果分为非适生区和低、中、高适生区，再对分类结果进行重分类，采用两种分类方法探究其适生区变化情况：1.将分类结果分为适生区和非适生区，适生区包括低、中、高适生区，即蛇足石杉的总适生区；2.将非适生区和低适生区合为一类，中适生区和高适生区合为一类。使用SDM Toolbox计算不同时期各适生区面积变化和总适生区质心变化情况，通过ArcGIS 10.6对结果进行可视化表达。

2 结果与分析

2.1 模型预测精度评估

基于AUC值和TSS值评估单一模型和组合模型的预测精度(表3)，除SRE算法构建的模型外，所有模型均AUC值大于0.8且TSS值大于0.7，说明单一算法构建的模型能较好预测蛇足石杉适生区，但组合模型AUC值和TSS值大于所有单一模型，预测精度更高。组合模型AUC值大于0.9，说明蛇足石杉适生区预测结果极好。

表3 组合模型及单一模型AUC值与TSS值

2.2 生物气候变量重要性评估

由生物气候变量重要性评估结果(图2)可知，最湿月降水量(Bio13)是对蛇足石杉适生区影响最大的生物气候变量，此外重要性由高到低依次是温度季节性(Bio4)、最冷季均温(Bio11)、最暖月最高温(Bio5)、降水季节性(Bio15)和等温性(Bio3)。总体上看，降水和温度都对蛇足石杉适生范围有重要影响。

Bio3.等温性;Bio4.湿度季节性;Bio5.最暖月最高温;Bio11.最冷季均温;Bio13.最温月降水量;Bio15.降水季节性

2.3 蛇足石杉适生区预测

蛇足石杉当前适生区主要位于中国南部(图3)。从适宜性来看，蛇足石杉在贵州、重庆、湖南、福建及浙江等地区适宜性较高，在新疆、青海、甘肃、黑龙江及宁夏等地区适宜性很低。从适生区来看(表4)，当前总适生区面积226.99×104km2，其中：低适生区占33.69%，主要位于四川东部、广西、湖南东部、广东东部及江西等地区；中适生区占32.97%，主要位于云南和四川东部，在其他省区少量散布；高适生区占33.34%，主要位于贵州、重庆、湖南西部、福建及浙江等地区。4种未来气候情景下适生区预测结果如图4，低适生区持续扩张，在RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下分别扩张31.47%、50.17%、55.01%和84.14%。中适生区和高适生区不断收缩，中适生区在RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下分别收缩6.03%、13.04%、15.55%和31.90%，高适生区在RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下分别收缩17.54%、28.65%、43.23%和58.69%。总适生区变化不大，在RCP2.6及RCP4.5情景下分别扩张2.77%和3.06%，在RCP6.0及RCP8.5情景下分别收缩1.00%和1.73%。

表4 当前及未来(2070)气候情景下的蛇足石杉适生区面积

图3 当前蛇足石杉适生区

图4 未来(2070)气候情景下的蛇足石杉适生区

2.4 蛇足石杉适生区及质心变化

蛇足石杉总适生区面积净变化不大，总变化与排放水平呈正比(图5，表5)。在RCP2.6及RCP4.5情景下，蛇足石杉适生区面积分别扩张5.80%和8.46%，收缩3.32%和5.90%，净扩张。在RCP6.0及RCP8.5情景下，蛇足石杉适生区面积分别扩张7.33%和10.69%，收缩8.82%和12.96%，净收缩。在适生区扩张方面，4种情景都向云南、四川、辽宁及吉林等地区扩张；在适生区收缩方面，4种情景下江苏和湖北等地区都存在明显收缩情况，此外，RCP6.0及RCP8.5情景下湖南及安徽等地区适生区也明显收缩。

图5 未来(2070)气候情景下蛇足石杉总适生区变化

表5 未来(2070)气候情景下蛇足石杉总适生区面积变化

蛇足石杉中、高适生区面积随排放水平由低到高持续收缩(图6，表6)。在RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下分别净收缩14.34%、23.31%、31.63%和47.08%，四川东部、广东、广西、湖南及福建等地区的中、高适生区向低适生区和非适生区转变，四川中部和云南北部等地区的少量低适生区和非适生区转变为中、高适生区。

表6 未来(2070)气候情景下蛇足石杉中、高适生区面积变化

图6 未来(2070)气候情景下蛇足石杉中、高适生区变化

蛇足石杉总适生区质心和中、高适生区质心均在湖南省西部，在RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下质心均向西北转移，与适生区变化情况一致(图7)。

A.蛇足石杉总适生区质心变化；B.蛇足石杉中、高适生区质心变化

3 讨 论

3.1 蛇足石杉的未来分布

蛇足石杉在中国分布广泛，但准确的分布记录极不均匀，在山东、江苏及湖北等地区的精确分布数据尤其少甚至没有，物种分布模型的准确性与物种分布数据的数量和精确程度密切相关[39]，本研究筛选精确到村级单位和有明确经纬度及地址的蛇足石杉分布数据，进行经纬度查询、校对和空间过滤后用于构建模型，尽可能确保准确性，组合模型AUC值均大于0.9且与实际分布相符也说明预测结果可信度高, 但是分布数据有待未来通过完善标本记录和准确目击记录进一步优化。此外，物种与环境之间的关系存在复杂性和不确定性，对气候变化的响应存在滞后性[40]，增加了物种分布模型预测的难度。因此，物种分布模型预测的是适合物种生长的区域，即物种的潜在分布区域。本研究主要基于对19个生物气候变量去相关性获得的6个变量构建物种分布模型，而在物种实际生存中，海拔、坡度及地质等环境变量都可能会影响分布[41]，这些变量在深入研究中需要纳入考量。

在二氧化碳与气溶胶排放水平逐渐升高的RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下，蛇足石杉总适生区面积变化不大，但蛇足石杉中、高适生区剧烈收缩，收缩面积与排放水平成正比，说明高排放对蛇足石杉中、高适生区有消极影响，迫切需要合理调控排放保护蛇足石杉适生区。在RCP2.6、RCP4.5、RCP6.0及RCP8.5情景下，蛇足石杉适生区质心均向高纬度转移，与当前大部分研究结果一致，随着温度升高物种适生区会向高纬度地区迁移[42-46]，与Shafer等[47]研究发现气候变化会使得物种向各个方向迁移相符。

3.2 蛇足石杉的保护

蛇足石杉因含有药用成分石杉碱甲，野生资源被大量采挖作为药用材料，现在中国将蛇足石杉列为国家二级重点保护野生植物，为蛇足石杉的保护提供了政策依据。结合研究结果来看，合理控制温室气体排放也有利于蛇足石杉的生存及其生物资源的保护。此外，蛇足石杉高适生区主要在贵州、重庆、湖南及浙江等地区，其高适生区质心在湖南靠近贵州方向，且存在向贵州方向转移的趋势，建议重点加强湖南及贵州的蛇足石杉野生资源的保护工作。

4 结 论

通过构建组合模型预测珍稀植物蛇足石杉当前和未来(2070)气候情景下在中国的适生区并计算适生区及其质心变化，结果显示，蛇足石杉当前总适生面积为226.99×104km2，中、高适生区面积154.51×104km2，质心均位于湖南省西部。蛇足石杉未来总适生区面积变化不大，中、高适生区剧烈收缩，收缩面积与温室气体和气溶胶排放水平成正比，预测未来蛇足石杉中、高适生区将向非适生区和低适生区转变。温度和降水都对蛇足石杉适生区有重要影响，最湿月降水量(Bio13)、温度季节性(Bio4)、最冷季均温(Bio11)和最暖月最高温(Bio5)是对其适生区影响最大的气候变量。