受气候变化及人类活动的影响, 生物多样性的丧失已成为全球面临的最严重的问题之一, 尤其是对扩散能力及对气候变化的适应均有限的植物类群^[1], 未来气候变化将进一步加剧这些类群的灭绝风险^[2]。目前气候的持续变暖, 导致环境因子发生不同程度的改变, 引起植物种群作出不同变化或响应^[3-4]。相比低海拔地区, 分布于高海拔区域的植物响应更为强烈, 如珙桐、附子、油松等, 由于温度升高、融雪提前、导致这些植物(珙桐、附子、油松)生长、授粉和繁殖的时间发生变化, 生态位发生显著变化, 自然分布面积急速缩减, 并且表现出向高纬度区域迁移的趋势^[5-7]。高山植物类群对气候变化的不同响应模式, 将引起高山生态系统的物种组成和群落结构发生改变, 导致区域高山本土植物或特有植物种群消失或消亡。因此, 本研究探究西北地区高山植物的生态适宜性及对未来气候变化情景的可能响应, 将指导区域高山植物多样性的科学管理与保护。

近些年, 随着气候变化情景的深入研究以及地理信息系统(GIS)技术的飞速发展, 利用物种分布模型预测不同生物类群在气候变化情景下的适宜生境变化成为解决景观生态学和保护生态学相关热点问题的主要途径之一。物种分布模型(SDMs), 基于物种已知分布点和环境因子数据模拟和预测物种在地理空间的分布状况^[8-9], 已经成功地应用于生态学、生物地理学、生物多样性保护等领域^[10-11]。目前, 被广泛应用的物种分布模型主要有：最大熵模型(MAXENT)、广义线性模型(GLM)、广义加性模型(GAM)、多元自适应回归样条(MARS)、随机森林(RF)等。其中, 最大熵采用机器学习算法, 是目前最流行的建模方法之一, 相比其它物种分布模型, 该方法已经被许多研究证实具有较优秀的预测能力和精度^[12-13]。近些年, MAXENT模型与GIS技术结合对很多高山植物开展了未来气候情景下适宜分布区变化的预测, 表现出较强的预测准确性与空间表现力, 预测结果对高山地区生物多样性的保护具有重要的指导意义^[14-15]。

天山花楸(Sorbus tianschanica), 蔷薇科(Rosacea), 花楸属(Sorbus), 落叶灌木或小乔木, 常生于海拔2000-3200米的高山溪谷或云杉林边缘^[16]。在我国, 天山花楸的自然分布主要集中在新疆阿尔泰山、天山以及青海与甘肃高海拔区域的针阔叶林下、林缘^[17]。天山花楸是极具发展前景的高山植物, 具有良好的观赏、药用、经济、生态及科研价值。但是, 近些年, 由于过量采伐以及种子萌发受限等原因, 该植物在森林的比重中大幅降低^[18], 原生地生境破碎化严重, 资源量显著下降, 亟待保护。但是, 目前对天山花楸已开展的研究主要体现在化学成分提取、种子育苗栽培等方面^[19-20]。该植物在我国西北地区的的生态适宜性、适宜分布特征及其环境驱动机制都并不清楚, 成为限制天山花楸资源利用和科学保护的瓶颈。

因此, 本研究利用天山花楸的地理分布信息和相关环境数据, 基于MAXENT模型、GIS工具、R程序包, 利用统计学和空间统计学方法分析不同气候情境下天山花楸在我国西北地区的适宜分布。主要解决以下科学问题：(1)不同气候情景下, 天山花楸在我国西北地区的适宜分布范围与空间分布特征；(2)限制天山花楸分布的关键驱动因子；(3)天山花楸的生态适宜性及其对未来气候变化情景的可能响应。

1 数据与方法 1.1 分布数据来源与处理

本研究通过检索中国数字植物标本馆、全球生物多样性信息网络结合已发表文献^[21-22], 共获得212个天山花楸的自然分布点。通过BIGMAP地图下载器和中国地名录确定详细地理坐标, 基于R中的raster和biomod2程序包进一步剔除重复和模糊的分布点, 矫正分布点偏差。最终, 获得51条分布记录用于模型的模拟, 覆盖了该植物在新疆、甘肃和青海的自然分布范围(图 1)。

1.2 环境数据来源与处理

本研究中, 19个气候因子来源于IPCC生物气候数据集(http://www.worldclim.org)。选取2050时段CCSM4模型的气候数据。该模式包括RCP2.6、RCP4.5、RCP6.0和RCP8.5四种排放情景, 本研究选择中等排放情景(RCP4.5)。坡度、坡向利用ArcGIS软件空间分析工具从中国科学院地理空间数据云平台提供的数字高程模型(DEM)中提取。归一化植被指数(NDVI)数据采用平均后的中国科学院地理空间数据云平台提供的所有中国地区500m分辨率月值。模型模拟使用的中国省级行政区划底图来源于国家基础地理信息系统网站。

考虑空间环境变量对物种的限制作用, 降低变量的空间纬度避免过拟合现象, 运用R中的Hmisc程序包对23个环境因子进行Pearson相关性分析^[23]。剔除相关性高(>80%)的变量, 最终得到10个环境因子, 包括7个气候因子：平均气温日较差(Bio2)、温度季节性变化标准差(Bio4)、最冷月最低气温(Bio6)、最干季平均气温(Bio9)、最湿月降水量(Bio13)、最干月降水量(Bio14)、降水量季节性变化(Bio15)；2个地势因子：坡度、坡向以及1个植被因子：NDVI。根据研究区域, 利用ArcGIS裁剪环境图层数据, 用于后续模型的模拟分析。

1.3 模型构建及环境因子分析

MAXENT模型构建过程中, 随机选取75%的分布点作为训练模型, 剩余25%的分布点用于验证模型, 进行10次重复^[24]。为了方便模型模拟结果的可视化表达, 参照已发表的文献^[25], 利用ArcGIS空间分析工具计算所有分布点在潜在分布图上的适宜性, 并计算这些值的标准差σ和平均值μ, 选择p=σ－μ作为阈值, 区分“适生区”和“不适生区”；选择p=σ+μ作为阈值, 区分“适生区”和“高度适生区”。进一步明确未来气候变化下天山花楸适宜分布区的迁移趋势, 利用ArcGIS计算不同气候下该植物适宜分布区的质心。

利用MAXENT模型生成的环境变量贡献率和置换重要值分析各环境因子对模型模拟的贡献以及模型对各环境因子的依赖；利用刀切法(Jackknife)分析各环境因子在建立模型时的重要性^[26]。利用ArcGIS提取裁剪后环境图层数据, 得到上述7个气候因子的生态位参数, 并通过R的ggplot2程序包绘制核密度图, 量化分析。

1.4 多元环境相似度面和最不相似变量分析

利用多元环境相似度面(MESS)和最不相似变量(MoD)分析未来时段相比基准时段研究区的气候变化及其主要变化因子。利用MAXENT.jar文件中的density.tools.Novel工具计算基准与2050时段研究区内7个气候因子的相似度(S), 得到多元环境相似度面(MESS)。当S值为正值时, S值越小表示该区域的气候因子的差异越大, S值为100.00时, 表示没有差异；当S值为负值时, 表示该区域至少有1个气候因子的不相似度超出了参考范围, 该区域的环境差异极大^[27]。同时, 提取相似度最低的气候因子, 得到最不相似变量(MoD), 用ArcGIS进行可视化处理。

2 研究结果与分析 2.1 天山花楸的适宜分布及其对未来气候变化的响应

基准及2050气候情境下, 天山花楸预测模型的AUC值分别为0.9336和0.938, 表明模型模拟效果较准确。基准气候下, 天山花楸的适宜分布较破碎, 占研究区总面积的13%, 适宜分布集中于高山及高海拔区域。基于本研究的阈值划分原则, 将研究区划分为不适生区(0-0.279)、适生区(0.279-0.734)和高度适生区(0.734-1)。具体看, 在新疆, 天山花楸的适生区主要分布于阿尔泰山中西段、准噶尔西部山地、天山中西段、博格达山以及帕米尔高原和喀喇昆仑山北段的片段化区域；在甘肃和青海, 该植物集中分布于河西走廊中东部、祁连山以及南山。其中, 高度适生区均为较小的破碎斑块, 占研究区总面积的2.39%, 分布于新疆阿尔泰山西段、准噶尔西部山地、天山西段及祁连山中段(图 1)。

2050时段, 天山花楸的适宜分布范围相比基准气候将略微增加(3.65%), 将在阿尔泰山中段、准噶尔西部山地、天山中西段及祁连山中段扩增, 且整体适宜分布呈向北迁移趋势。具体表现为：该植物在其适宜区北部的阿尔泰山及准噶尔西部山地的面积明显增加, 而在适宜区南部河西走廊东部的面积显著减少。此外, 天山花楸在新疆的适生区将向西北方向迁移；在甘肃和青海的适生区将向东北方向迁移。(图 1)。

2.2 研究区的气候波动及影响天山花楸的关键因子分析

多元环境相似度面显示, 2050时段比基准气候研究区绝大部分地区的气候波动较明显, 气候因子相似度低于20的区域占研究区总面积的72%。2050时段相比基准气候, 准噶尔西部山地、天山西段、塔里木盆地东南部均表现出强烈的气候波动, 气候因子相似度均低于0(图 2)。天山花楸的已知分布区域也表现出明显的气候波动, 气候因子的平均相似度为17。最不相似变量显示, 最冷月最低气温、温度季节性变化标准差和降水量季节性变化将发生显著变化, 是造成研究区内67%的区域发生气候波动的主要因子(图 2)。最干月降水量、最湿月降水量、平均气温日较差和温度季节性变化标准差是影响天山花楸82%的分布点气候波动的主要因子。

2.3 影响天山花楸分布的驱动因子及其变化范围

MAXENT模拟结果显示, 最干月降水量、最湿月降水量和平均气温日较差是影响天山花楸适宜分布的关键驱动因子。用于模型模拟的10个环境因子中, 归一化植被指数、最干月降水量和最湿月降水量对模型的贡献较大, 累积贡献率达71.7%。置换重要性值的结果表明, 模型对归一化植被指数、最干月降水量和坡向的依赖性更强(表 1)。刀切法检验结果也表明(图 3)：归一化植被指数对天山花楸模型模拟结果的影响较大；降水因子(bio13?#322;io14?#322;io15)在正则化训练增益、测试增益及受试者工作特征曲线下面积的贡献都高于气温因子(bio2?#322;io4?#322;io6?#322;io9)。其中, 降水因子的最湿月降水量贡献最高, 气温因子的平均气温日较差贡献最高；地势因子中, 坡度对模型的重要性较高, 表明这些因子对天山花楸适宜分布的模拟结果影响较大。

天山花楸适宜分布区7个关键气候因子的数值范围(图 4)分析结果表明：基准气候下, 天山花楸适宜分布区气候因子的数值范围为：最湿月降水量6-127mm；最干月降水量0-18mm；气温日较差8.2-16.3℃。2050时段, 相对基准气候, 天山花楸适宜分布区的降水量明显增加, 极端气温因子将略有升高, 最干季平均气温和最冷月最低气温将分别上升0.7℃和0.15℃。

3 讨论 3.1 天山花楸的适宜分布及其驱动因子

本研究表明, 基准气候下, 天山花楸的适宜分布面积较小且破碎, 主要集中在阿尔泰山西段、准噶尔西部山地、天山西段及祁连山中段的高海拔山地区域(图 1), 这与天山花楸主要分布于高山高海拔区域的记载相符^[28]。其中, 模型模拟的天山花楸在甘肃祁连山东中段区域具有高度适生生境, 文献记载也表明该植物在祁连山东中段具有大片的天然分布^[29]。相关研究表明天山花楸在天山一带未发现纯林, 均与天山野苹果、天山云杉等植物混生, 天山野苹果与天山云杉常分布于天山北麓的西坡、北坡、西北坡以及昆仑山北坡^[30], 而本研究的模型识别出该植物在上述区域均具有适宜和高度适宜分布, 说明本研究模型的预测具有较高的准确性。

本研究的环境因子分析均表明NDVI对天山花楸模型的模拟结果影响较大, 响应曲线表明植被指数与模拟的适生程度存在显著正相关关系(图 5)。本研究模型构建时加入NDVI, 显著提升了模型的模拟精度和空间表现, 模型识别的适宜分布区的关键气候变量的数值范围相比未考虑该因子的模拟结果更接近该植物实际的分布条件^[31]。相关研究表明植被遥感数据在区域尺度上能显示出较高的空间异质性, 和气候变量结合能大大提高生态位模型的应用效果^[32]。此外, 本研究模拟中加入坡度因子也进一步提高了天山花楸模型模拟的识别度(表 1), 这可能是因为地形起伏造成的小生境异质性进一步区分了研究区环境的异质性, 使获得的模型模拟结果更加精细、破碎。土壤因子、地形因子等非生物因子的加入也被表明能提高环境异质性, 从而提升模型对植物类群模拟的准确度^[33]。

本研究表明主要限制天山花楸适宜分布的气候因子为降水条件(表 1、图 3)。这与天山花楸自身的生物学特性有关, 天山花楸属浅根系树种, 根系多分布在0.3-0.6m的土层内, 不能深入土层汲取水分, 且叶片具皮孔, 水分蒸腾快, 不能很好地储存水分^[34-35]。相关文献记载, 天山花楸分布的生境多集中于永久性或季节性流水的小溪旁, 对含水量低的环境适应能力较差^[36]。此外, 天山花楸面积扩增区域的最不相似变量也为降水因子, 且核密度图显示最干月降水量和最湿月降水量明显增加(图 4), 表明降水量的增加利于天山花楸适宜分布区的局部扩展, 如准噶尔西部山地, 祁连山中段由于2050时段降水量的增加高度适生区显著扩张。这可能由于西北地区作为全球同纬度最干旱的地区之一^[37], 降水因子也是影响该地区植物生长的主要环境因素。此外, 温度因子对天山花楸的局部适宜生境也有一定的影响, 累积贡献率达11.4%, 如天山花楸在天山、河西走廊、帕米尔高原的片段化适宜分布区的温度2050时段将明显升高(图 2、图 4), 将导致该植物在天山北坡中段、河西走廊东部及帕米尔高原中段区域的适宜分布进一步破碎、减少。

3.2 天山花楸对未来气候变化的可能响应

2050时段, 天山花楸的适宜分布将向高纬度高海拔地区迁移, 适生范围仅在局部高海拔山地：阿尔泰山及准噶尔西部山地、天山及祁连山区略有扩展, 这可能与未来气候下降水将在中亚北部地区略有增加有关^[38]。在气候变化背景下, 北半球大多数高山植物的分布区将发生迁移, 如：对美国东北部134种乔木树种的适宜分布的研究结果表明, 受未来气候变化的影响, 超过3/4的物种的适生区将向北迁移^[39]；运用响应面曲线分析的未来气候变化情景下美国西北部乔木和灌木类群不仅向高纬度和高海拔地区扩散, 还将向多个方向迁移^[40]。

预测天山花楸对未来气候变化的可能响应, 能够有针对性的对其进行保护。本研究的结果表明未来气候情景下天山花楸适生分布区的面积有少量增加, 但物种的实际分布受到多种因子的制约, 实际增加面积未必能达到模型的预测值。此外, 将获得的该植物的适宜分布区与其适宜分布的海拔范围(2000-3200m)叠合、残余分析后的适生区更加破碎, 适宜面积和高度适宜面积仅占研究区总面积的4%和0.84%(图 5)。然而, 目前, 天山花楸植物资源十分稀少, 需要通过人工繁育扩大种植数量, 才能满足药用资源、城市绿化等的需要, 且天山花楸原生地无纯林, 多年生老树周围也极少发现幼树, 属于萌发困难、更新能力较差的树种^[41], 加之被过度采挖, 天山花楸自然资源量正迅速下降。尽管如此, 天山花楸目前尚未存在自然保护区, 仅在祁连山地区有研究者引种栽培, 但天山花楸的自然分布大多集中在新疆天山。因此, 为保护该植物资源, 本研究建议将模型识别的天山花楸在阿尔泰山西段、准噶尔西部山地、天山西段及祁连山中段的适宜分布区, 都应作为保护与管理的重点区域。并且, 该物种在河西走廊东部、青海南山的高度适宜区在2050时段气候情景下将显著减少, 这些敏感性区域更应作为对该植物进行保护的高风险区域, 很有必要开展原地和异地保护。因此, 本研究模拟预测的天山花楸的适宜分布区及其生态适宜性变化直接服务于该植物资源的科学管理和保护, 也利于高山区域的生态修复。