气候变化对生物系统多样性的影响是近年来研究的热点领域^[1]。联合国政府气候变化专门委员会(IPCC)发布的第五次评估报告(AR5)及《IPCC全球升温1.5℃特别报告》指出, 未来全球气温仍持续升高, 对发展中国家可能有更为严重的影响^[2-3]。气候变化直接或间接的影响着生物多样性^[4], 使物种和植被的水平和垂直分布发生了巨大变化。因此研究物种对气候变化的响应以及预测气候变化情景下物种的潜在地理分布和迁移, 对该物种的保护和利用具有非常重要的作用^[5-6]。

生态位模型(ecological niche models, ENMs)是以生态位理论为研究基础, 利用模型分析研究对象的已知分布点及其相关的环境数据, 从而来预测物种在不同因子下的实际分布和潜在分布^[7]。目前常用的生态位模型有：生物气候分析系统(bioclimate analysis and prediction system, BIOCLIM)^[8]、生态位因子分析模型(ecological niche factor analysis, ENFA)^[9]、基于规则集的遗传算法模型(genetic algorithm for rule-set production, GARP)^[10]和最大熵模型(maximum entropy modeling, MaxEnt)^[11-12]等。其中MaxEnt模型与其他模型相比, 运算准确性高, 运算时间短, 较少的有效分布点也可建模, 所以被广泛用于动物栖息地预测^[13]、病虫害防治^[14]、濒危物种预测与保护^[15]等领域的研究。该模型多次运用于预测干旱区植物的适宜分布区及变迁^[16-18], 说明该模型适合对干旱区物种分布进行研究。

绵刺(Potaninia mongolica Maxim.)为绵刺属(Potaninia)蔷薇科(Rosaceae)的强旱生小灌木, 分布于亚洲中部西鄂尔多斯荒漠地带狭小区域, 是第三纪古地中海孑遗植物, 国家首批二级濒危保护植物^[19-21]。绵刺为荒漠建群种植物, 耐旱能力极强, 能通过休眠的方式躲避极端干旱气候^[22]。繁殖方式分为有性和无性两种, 野生绵刺大多以劈裂生殖为主^[22-24]。在荒漠地区, 绵刺不仅发挥着重要的生态作用, 而且是一种天然饲料, 青鲜时牲畜喜食, 因此具有重要的科研价值和应用价值^[25]。目前关于绵刺的研究主要集中在生理生态学特性^[26]、胚胎发育^[27]、克隆生长^[28]、劈裂生长^[23]、繁殖对策^[29]、遗传多样性^[30]、绵刺根系形态对土壤水分响应的研究^[31]、光合特性的研究^[32]以及群落结构和物种多样性方面^[33]。近年来, 由于过度放牧及人为干扰, 绵刺生长缓慢, 繁殖力差, 分布范围日趋缩小, 数量锐减, 处于濒临灭绝的危险^[33-34]。目前国内外对于绵刺基于生态位模型在过去、当代和未来气候变化下潜在适生区预测的研究较少。因此, 本研究基于19个环境气候因子, 利用MaxEnt模型, 结合ArcGIS、SPSS等软件预测绵刺在末次间冰期(Last Interglacial)、末次盛冰期(Last Glacial Maximum)、当代(Current)和未来(2050年、2070年)五个不同气候情境下的潜在适生区及迁移路线, 分析影响我国绵刺分布的主导环境因子, 为我国荒漠地带绵刺保护、利用和维持生态环境提供一定的科学依据。

1 数据来源与研究方法 1.1 物种分布数据来源

绵刺当前在中国的地理分布数据来自于公开发表的相关文献和中国植物数字标本馆(CVH, http://www.cvh.org.cn), 利用奥维互动地图对收集的数据点进行筛选, 删除重复的或者经纬度不精确的点, 得到85条分布准确数据点。将85条分布数据点导入ENMTools软件, 并结合2.5弧分精度(5 km×5 km)的环境数据筛选分布数据, 每个栅格中只保留一个分布点, 最终获得73个有效分布点数据(图 1), 在excel中保存为.csv格式备用。

1.2 地理和环境数据

地理数据在国家基础地理信息系统(http://nfgis.nsdi.gov.cn/)中下载1 ∶ 400万的中国行政区划矢量图作为分析的底图。从世界气候数据库(http://www.worldclim.org/)下载末次间冰期(Last Interglacial)、末次盛冰期(Last Glacial Maximum)、当代、2050年、2070年五个时期包括19个有关气候和降水的环境气候数据(表 1)。空间分辨率为2.5arc-minutes, 采用的是WorldClim数据库中美国大气研究中心(NCAR)开发的全球气候模式(CCSM4)^[35]。未来2050年、2070年气候变化情景包括典型浓度路径4.5(RCP4.5)和典型浓度路径8.5(RCP8.5)两种温室气体排放情景的气候数据。RCP4.5和RCP8.5分别为未来温室气体排放中等等级和最高等级情景^[36], 可以对未来气候的变化进行较为准确的预测。

1.3 分布模型建立与研究方法

运用MaxEnt模型模拟绵刺在末次间冰期(Last Interglacial)、末次盛冰期(Last Glacial Maximum)、当代、2050年、2070年五个不同时期在中国的潜在适生区。首先, 利用ArcGIS 10.4软件将绵刺分布点数据和19个环境因子转化为ASC Ⅱ格式后导入MaxEnt 3.3.3k软件中, 随机选取25%的分布点作为测试数据(test data), 剩余75%的分布点作为训练数据集(training data), 迭代运算500次, 重复运行10次, 其他参数设置为默认, 删除贡献率为0的6个环境因子；其次, 在ArcGIS 10.4软件中利用多值提取模块, 将剩余的13个环境因子提取至样点；最后, 利用SPSS 23.0进行Pearson相关系数检验, 对于相关系数大于0.80的两个环境因子^[37], 保留贡献率较高的一个, 最终选取了8个环境因子参与建模(表 1)。

用ArcGIS 10.4软件进行可视化处理, 再利用自然断点法进行重分类划分, 将绵刺的分布区划分为四个等级：高适生区(0.4-1)、中适生区(0.2-0.4)、低适生区(0.05-0.2)、不适生区(0-0.05)四个等级, 并绘制出不同时期绵刺在中国的潜在适生区分布预测图^[38]。

用ArcGIS 10.4中的SDM工具箱分析绵刺不同时期在中国的潜在地理分布变化及适宜生境中心的变化。得出不同时期绵刺的稳定、扩张、收缩区域面积, 利用向量计算不同时期绵刺的适宜生境变化的位置及方向, 绘制出绵刺在不同气候情景下的适宜区的迁移路线图^[35]。

1.4 MaxEnt模型精度检验

MaxEnt模型采用受试者工作特征曲线(receive operating characteristic curve, ROC)分析法进行模型精度检验。ROC曲线与横坐标围成的面积为AUC值, AUC值用来判断预测结果的准确性。AUC值数值范围是0-1, 数值越大, 说明预测结果的准确性越高。通常, AUC值≤0.6, 认为预测结果失败；0.6 < AUC值≤0.7, 认为预测结果较差；0.7 < AUC值≤0.8, 认为预测结果一般；0.8 < AUC值≤0.9, 认为预测结果精准；0.9 < AUC值≤1.00认为预测结果极精准^[39]。

2 结果与分析 2.1 MaxEnt模型精度检验

利用MaxEnt 3.3.3k软件对绵刺当前73个分布点和8个环境变量在中国的潜在适生区进行模拟, 10次重复, 平均AUC值为0.988(图 2), 其他时期的AUC值也均在0.9以上, 表明MaxEnt模型预测结果极精准。

2.2 影响绵刺地理分布的主要环境因子

利用MaxEnt模型计算环境变量对物种分布的相对贡献率及SPSS 23.0进行Pearson相关检验。结果表明影响绵刺地理分布的主要环境因子有：最干月降雨量(bio14)、年均降雨量(bio12)、最热月的最高温(bio5)、等温性(bio3)、最冷季度平均降雨量(bio19)、最湿月降雨量(bio13)、最暖季度平均温(bio10)、最冷季度平均温(bio11), 累计贡献率达97.4%。其中贡献率最大的是最干月降雨量(bio14), 贡献率为47.9%；最热月的最高温(bio5)次之, 贡献率为24.8%；第三位的是年均降雨量(bio12), 贡献率为21.6%(图 3)。

Maxent模型刀切法被用于判断主导物种分布的主要环境变量因子^[40], 刀切法检验结果显示(图 4), 正规化训练增益最高依次是年均降雨量(bio12)、最湿月降雨量(bio13)、最冷季度平均降雨量(bio19), 单变量等温性(bio3)正规化训练增益最小, 表明该环境因子对绵刺的地理分布预测的重要性较小。综上分析表明, 影响绵刺地理分布的重要环境因子是降水(最干月降雨量、年均降雨量、最湿月降雨量、最冷季度平均降雨量), 温度(最热月的最高温)。

根据绵刺的适生等级划分, 当绵刺为高适生区时：最干月降雨量(bio14)最适范围是0.1-0.5mm, 最热月的最高温(bio5)最适范围是28-31℃, 年均降雨量(bio12)最适范围是80-183mm, 最湿月降雨量(bio13)最适范围是20-60mm, 最冷季度平均降雨量(bio19)为0.5-3.5mm。

2.3 绵刺潜在地理分布区及变迁 2.3.1 当前绵刺在中国的潜在地理分布区

当前绵刺在中国的适生分布范围大约在北纬32.97°-42.12°, 东经71.78°-123.92°, 主要分布在我国内蒙古自治区, 甘肃省、宁夏回族自治区、新疆维吾尔自治区有少量分布, 西藏自治区、青海省、陕西省、山西省也有零星分布, 包括了所有的绵刺有效分布点, 表明在当代气候情景下模拟得到的绵刺潜在地理分布区与其实际分布范围非常吻合。当前绵刺在中国的潜在适生区面积约是51.94×10⁴ km², 约占国土总面积的5.41%, 分布中心在内蒙古阿拉善盟。其中高适生区面积约为11.58×10⁴ km², 约占总适生区面积的22.29%；中适生区面积约为14.52×10⁴ km², 约占总适生区面积的27.96%；低适生区面积约为25.84×10⁴ km², 约占总适生区面积的49.75%。其中高适生区主要分布于内蒙古阿拉善盟中部偏南、鄂尔多斯市东北部与巴彦淖尔市接壤处、宁夏西北部, 新疆中部、西部有少量带状分布, 甘肃省张掖市、白银市也有少量破碎状分布(图 5, 表 2)。

2.3.2 不同气候情景下绵刺在中国的潜在生境变化

在过去及未来各种气候情景下, 绵刺的潜在适生区面积呈现先降低后增加, 再到缓慢变化的趋势。末次间冰期潜在适生区面积约是203.71×10⁴ km², 其中高、中、低适生区面积分别约为54.38×10⁴ km²、55.57×10⁴ km²、93.75×10⁴ km²。到末次盛冰期, 潜在适生区面积大幅度减少, 高、中适生区面积几乎可以忽略不计, 下降幅度达97.69%。从末次盛冰期到当代, 潜在适生区面积恢复到51.94×10⁴ km²。当代到未来气候情景下, 潜在适生区面积缓慢变化, 但在2070年RCP8.5情景下高适生区面积下降到3.94×10⁴ km², 与当代相比, 下降幅度为65.98%(图 5, 表 2)。

从末次间冰期到末次盛冰期, 绵刺潜在适生区面积大幅度的收缩, 具体表现为高、中适生区几乎丧失, 低适生区也从最初的93.75×10⁴ km²减少到4.71×10⁴ km², 新疆、甘肃、青海、宁夏、四川、云南东北部、贵州西北部、陕西、山西、河北南部、河南西北部、山东等地都有明显的收缩。末次盛冰期到当代, 绵刺高、中、低潜在适生区有所恢复, 潜在适生区面积增加了47.23×10⁴ km², 新增区域主要集中在内蒙古, 宁夏、甘肃、新疆也有分布(图 6, 表 3)。

在未来气候情景下, 绵刺的潜在适生区主要集中在内蒙古中西部, 宁夏北部、甘肃中东部也有部分分布, 与当代相比, 绵刺的总适生区面积相对稳定, 适生区边缘地带有微小变化, 总体表现为适生区向西北地区扩张, 东北、东南方向收缩。在RCP4.5情景下, 2050年、2070年总适生区面积增加, 但高适生区面积减少, 2050年收缩地区主要在内蒙古赤峰及鄂尔多斯, 宁夏吴忠及新疆部分地区, 2070年与2050年相比, 部分边缘地带有微小的扩张和收缩。在RCP8.5情景下, 2050年总适生区面积增加, 但高适生区面积减少, 2050年收缩面积主要在内蒙古赤峰及鄂尔多斯, 宁夏吴忠及新疆部分地区, 2070年与2050年相比, 适生区向西南方向扩张, 高适生区面积收缩, 主要是内蒙、新疆、甘肃的边缘地带(图 6, 表 3)。

2.3.3 不同气候情景下绵刺在中国的分布中心变化

MaxEnt模型模拟表明, 绵刺从过去、当代到未来气候情景下, 分布中心整体变化是先向东北方向迁移, 再向西北方向迁移, 最后向西南方向迁移。绵刺在末次间冰期的分布中心位于甘肃省武威市天祝县(37.19°N, 102.76°E)(图 7)。从末次间冰期到末次盛冰期, 向东北方向迁移到宁夏石嘴山市平罗县(38.88°N, 106.14°E)。当前气候情景下, 分布中心位于内蒙古阿拉善右旗(40.19°N, 104.26°E), 与末次盛冰期相比, 分布中心向西北方向迁移。未来气候情景下, 绵刺分布中心都位于内蒙古阿拉善右旗, 总体向西南方向偏移, RCP 4.5情景下, 2050年分布中心位于北纬40.08°, 东经103.28°, 2070年分布中心位于北纬40.08°, 东经103.11°；RCP 8.5情景下, 2050年分布中心位于北纬39.99°, 东经102.88°, 2070年分布中心位于北纬39.16°N, 东经101.24°(图 7)。与当前相比, 未来气候情景下绵刺分布中心向西南、低纬度地区迁移, 有返回冰期的趋势。

3 讨论

绵刺是单种属植物, 系古地中海孑遗植物, 是西鄂尔多斯等干旱半干旱地区重要的荒漠植物, 具有一定的科学研究价值, 同时它又是一种天然饲料, 青鲜时牲畜喜食, 绵刺荒漠为中国的东阿拉善至西鄂尔多斯地区的主要放牧场地^[41]。当前绵刺分布区狭小, 由于过度放牧和矿区开采, 致使其种群数量日益减少。植被的不断破坏, 环境条件日趋恶劣, 绵刺种群的正常生长发育受到影响, 处于日益衰退的状态。当前及未来全球增温的情景下绵刺分布区的研究, 对我国绵刺物种保护和荒漠化防治有极其重要的现实意义。本研究利用MaxEnt模型和ArcGIS软件, 通过绵刺当前有效分布点和最终筛选出来的8个环境因子, 预测绵刺在我国末次间冰期、末次盛冰期、当前及未来气候情景下潜在分布区和分布区迁移方向。各气候情景下AUC值均在0.9以上, 表明MaxEnt模型预测结果较为精准。当前预测结果与绵刺在我国的实际分布相符^{[30, 42]}。

西部荒漠地区由于其独特的气候特点、地形地貌及古地理环境, 是许多珍稀植物, 特别是古老孑遗植物的“避难所”。很多学者已对该地区的孑遗濒危植物进行了未来气候变化下潜在地理分布研究, 例如段义忠等^[43]基于MaxEnt模型预测不同气候情景下沙冬青属植物在我国的潜在地理适生区, 结果表明未来2050年(RCP4.5、RCP8.5)沙冬青植物的潜在分布总面积有所增加, 但在RCP8.5情景下沙冬青属植物的最佳及高适生区面积将减少。马松梅等^[44]利用MaxEnt模型分析过去及未来蒙古扁桃的分布范围及空间格局的变化, 结果表明, 蒙古扁桃未来将会面对潜在分布范围扩展, 但最适生的分布区不断破碎、收缩的挑战。气候变化使得物种的地理分布及适生区范围发生不同程度的改变。鲁客等^[45]模拟了黑沙蒿未来的潜在分布区变化, 研究表明黑沙蒿在未来不同情景适生区面积均会缩减。本研究结果表明, 未来气候情景下, 绵刺的潜在适生区面积除2070年RCP8.5情景下少量减少外, 2050年(RCP4.5、RCP8.5)、2070年RCP4.5情景下均呈现增加的趋势, 但是潜在高适生区均收缩, 一部分高适生区转变为低适生区, 尤其是在2070年RCP8.5情景下, 高适生区由当代11.58×10⁴ km²减少到3.94×10⁴ km², 这说明大气中适量的温室气体浓度增加对绵刺的生存威胁较小, 然而当温室气体浓度过高时, 各环境变量受到影响发生剧烈变化, 长期高浓度的温室气体引起的气候变化将造成其适生区的减少和丧失, 最终可能影响绵刺的生存^[46]。由于未来气候不断变化, 绵刺的高适生区将减小, 生境片段化现象日趋严重, 该物种适生的存活生境面临严峻考验, 该物种灭绝的风险将增大。对绵刺野生群体而言, 未来气候变化下模拟的高适生区对于绵刺保护尤为重要。应对绵刺的核心分布区进行重点保护, 对未来将会消失的中低度适生区应划分保护区, 避免或者减少开发。

目前, 大多数学者认为, 在未来气候变化情景下, 随着全球温度的升高, 物种将呈现向高纬度和高海拔地区迁移的趋势^[47]。李艳红等^[15]对全球变化背景下南方红豆杉地域分布变化进行预测发现未来南方红豆杉将出现海拔向上迁移的趋势。张华等^[48]利用MaxEnt模型对未来气候变化背景下孑遗植物桫椤在中国的潜在地理分布分析表明, 桫椤的高适宜生境中心向西北地区和高纬度地区转移。本研究中, 绵刺的分布中心较当代整体向西南方向高海拔地区迁移, 甘肃、西藏、青海、新疆将逐渐出现新的适生区。说明不同物种对环境的适应能力不同, 因此在未来气候变化条件下不同物种的变化趋势不同。利用MaxEnt模型计算环境变量对物种分布的相对贡献率得出, 限制绵刺地理分布的重要环境因子有降水(最干月降雨量、年均降雨量、最冷季度平均降雨量、最湿月降雨量)和温度(最热月的最高温), 其中最干月降雨量范围仅为0.1-0.5mm, 可能是绵刺植物在最干月的生长状况对水分极具敏感性和依赖性^[26]。高温、干旱的夏季, 绵刺以“假死”的休眠状态存在；降雨时, 绵刺会打破休眠, 正常生长、开花、结果^[26], 因此绵刺对最干月降水和最热月的高温极度敏感。相比同区域生长的另一种植物四合木, 研究表明影响四合木分布的主要环境因子是最冷季度的平均降水量和年温的变化范围, 其次是降水量变异系数和温度季节性变化的标准差^[49]。

4 结论

本研究模型预测结果准确性较高, AUC值均在0.9以上。当前绵刺的潜在适生区主要位于我国内蒙古中西部地区, 宁夏东北部和西北部、甘肃中东部, 新疆中部和西部少量狭长分布。未来气候情景下, 绵刺的分布中心向西南方向小幅迁移, 高适生区面积减少, 甘肃、西藏、青海、新疆将逐渐出现新的适生区。降水(最干月降雨量、年均降雨量、最湿月降雨量、最冷季度平均降雨量)和温度(最热月的最高温)是影响绵刺地理分布的重要环境因子, 其中降水的作用更明显。内蒙古自治区鄂尔多斯市, 乌海市已经建立了西鄂尔多斯国家级自然保护区, 主要保护对象为四合木、半日花等古老残遗濒危植物和荒漠生态系统。建议相关部门设立绵刺保护点, 进行绵刺的繁殖试验和引种栽培。本研究的结论将为绵刺的保护利用提供科学依据, 并为合理建立绵刺种群提供参考。