气候变化情景下基于最大熵模型的青海云杉潜在分布格局模拟

曹雪萍,王婧如,鲁松松,张晓玮,*

1 甘肃农业大学林学院, 兰州 730070 2 兰州大学生命科学学院, 兰州 730000

青海云杉(Piceacrassifolia)为松科云杉属多年生常绿针叶乔木,是我国青藏高原东北缘特有树种,主要分布于中国甘肃、宁夏、内蒙古、青海等海拔1600—3800米地带[1]。青海云杉因其树形高大通直,材质轻软、纹理直、有弹性,是重要的工业用材,同时其抗旱性较强,生长迅速,是青海东部、甘肃北部山区和祁连山区优良的造林树种[2-3]。贺兰山区是青海云杉的遗传多样性中心,而其分布中心主要位于祁连山北坡[4]。在祁连山北坡,青海云杉多呈带状、块状与草原镶嵌形成森林-草原复合生态系统。上述山区是我国西北地区重要的生态安全屏障,作为其森林的优势树种,青海云杉在维系我国西北地区生态平衡、水土保持、水源涵养和生物多样性等方面发挥着重要作用[5- 8]。目前的研究主要集中在青海云杉生理学、生态学、育苗造林、病虫害防治,以及青海云杉林的林学特征、水源涵养和生态效益等领域[9- 11],但对该物种地理分布格局的研究主要集中于祁连山地区[8,12-13],而未见有关该物种整体分布区的相关模拟研究。且有关研究表明,祁连山中部正处于相对干旱和相对高温时期,并呈现出向暖干化发展的趋势[14],而降水减少和气候变暖均会严重影响青海云杉的生长[15-17]。

植物地理分布是由气候、水文、土壤、人类活动等多种因素长期共同作用的结果,其中,气候是决定植物在区域尺度上地理分布的最主要因素[18-20]。通过研究气候与物种分布的关系,可找出影响物种分布的主导因子,可探索物种分布形成的原因、确定物种的潜在分布区以及分析未来气候变化下物种分布区的变化情况,因此,物种的空间分布模拟一直是生物地理学的重要研究内容之一[12,21]。目前,因相关研究的侧重点不同,已建立了多种物种分布模型[21-24],其中以最大熵模型(Maxent)的预测能力较高,被广泛应用于濒危物种、入侵物种和重要经济物种的适生区预测和物种保护方面[25-26]。例如,独叶草(Kingdoniaunioflora)、毛红椿(Toonaciliatevar.pubescens)、红豆杉(Taxuschinensis)和水葫芦(Eichhorniacrassipes)等植物的研究[27- 30]。

本研究以青海云杉为研究对象,基于气候相似性原理,利用最大熵模型筛选影响青海云杉地理分布的主导因子,给出其潜在的地理分布范围,并预测其在未来气候变化背景下潜在分布区的变化,从而深化未来气候变暖对其潜在分布区变化影响的认识,以期为青海云杉的经营管理及应对气候变化提供依据。

1 材料方法

1.1 收集物种现有分布数据

青海云杉的现有分布样点通过在线查阅中国数字植物标本馆(http://www.cvh.ac.cn/)中已有记载的标本样点和收集实验室已有野外采集样点而得,并参考《中国植物志》中所述该物种的分布范围,对分布点进一步筛选,去除不确定和重复样点,最终确定69个样点,用于后续分析(见表1)。

表1 青海云杉样点地理信息

1.2 当前及未来潜在分布区预测

1.2.1获取气候数据

当前气候图层下载自WorldClim 网站(http://www.worldclim.org/),该网站提供了1970—2000年30间与降雨和气温有关的19个气候变量：年均温(bio1),月均温度变幅(bio2),等温性(bio3),温度季节性变化(bio4),最暖月最高温(bio5),最冷月最低温(bio6),温度年变幅(bio7),最湿季均温(bio8),最干季均温(bio9),最暖季均温(bio10),最冷季均温(bio11),年降雨(bio12),最湿月降雨(bio13),最干月降雨(bio14),降雨季节性变化(bio15),最湿季降雨(bio16),最干季降雨(bio17),最暖季降雨(bio18),最冷季降雨(bio19)。

未来20世纪50年代和80年代(2050s和2080s)的气候图层下载自CCCFS网站(http://www.ccafs-climate.org),气候增温情景则参照了Ren等人的方法[31],从CMIP5计划中选取了3种大气环流模型(General circulation models, GCMs: BCC-CSM1- 1、 CCCma_CanESM2和CSIRO-Mk3.6.0)和3种代表性浓度路径情景(Representative concentration pathways, RCPs: RCP2.6、RCP4.5和RCP 8.5))。这3种情景代表了不同的温室气体排放策略：RCP2.6下通过提高能源使用率和使用可再生新能源替代化石燃料等方法降低温室气体排放量,到2100年控制中国平均气温增幅低于2℃(缓和情景)；RCP4.5和RCP8.5则分别为中度稳定和高度温室气体排放情景,两种情景下,到2100年中国平均气温分别增加2.88℃和5.51℃[32-33]。以上数据图层均采用2.5′的空间分辨率。

1.2.2潜在分布区预测

图1 基于Maxent 模型预测青海云杉当前潜在分布的受试者工作特征曲线(ROC)曲线Fig.1 The receiver operating characteristic curve (ROC) curve predicting the current potential distribution of Picea crassifolia based on Maxent图中AUC表示受试者工作特征曲线下面积Area under the receiver operating characteristic curve

物种潜在分布区的预测基于最大熵模型(Maxent)来完成。Maxent通过物种现有分布数据可构建物种当前潜在分布模型,并能预测物种在气候变化下未来的潜在分布,是目前应用最广泛的物种分布模型[34-35]。将以上收集的青海云杉现有分布数据和当前及未来2050s和2080s气候图层导入软件Maxent,使用软件默认参数,计算得到青海云杉当前和未来潜在分布模型。在此过程中,从数据基中抽取75%作为训练数据,剩余25%作为交互数据以评估模型的有效性。模拟结果的精确度通过受试者工作特征曲线(Receiver operating characteristic curve,ROC)进行评估,用曲线下面积(Area under the curve, AUC)表示,AUC值介于0—1.0,其值为0.5—0.7时模型预测精度较差,0.7—0.9时为一般,值高于0.9时预测精度较好。各气候因子对青海云杉物种潜在分布的贡献率通过Jackknife模块进行评估。

此外,结合植物志描述的实际分布区和野外调查结果,本研究将青海云杉当前和未来潜在分布区域划为4个等级：<0.1为非适生区；0.1—0.3为低度适生区；0.3—0.5为中度适生区；0.5—0.8为高度适生区。通过分别计算各等级对应分布区面积,以明确气候变化对青海云杉实际分布的影响。另外为对比青海云杉未来潜在分布区总面积与当前的差异,本研究在进行模型运算时采用了第10分位训练存在阈值(10thpercentile training presence threshold),以获取物种分布的二元分布图像(仅指示适生与非适生两种状态)。以上分布区面积的计算和图形绘制皆在ArcGIS 10.3中完成。

2 结果

2.1 影响我国青海云杉地理分布的主要气候因子

基于Maxent 模型预测青海云杉当前潜在分布的ROC曲线显示,其训练集和测试集的曲线下面积(AUC)均高于0.99,表明该模型对青海云杉潜在分布预测的准确度很高(图1)。通过Jackknife模块对各气候因子的贡献率进行评估可知,等温性(bio3)和最冷季降雨(bio19)的训练得分值为最高,其次为最干季降雨(bio17)、最干月降雨(bio14)、降雨季节性变化(bio15)、最冷月最低温(bio6)、最干季均温(bio9)和最冷季均温(bio11)(图2)。绘制以上8个气候因子对青海云杉适生度的响应曲线可以发现,这些曲线均存在显著峰值,且最干月、最干季和最冷季降雨的适生范围显著低于其他因子,表明一年中的最低降雨量是影响青海云杉分布的主导因子(图3)。

图2 各气候因子对青海云杉当前潜在分布模型的Jackknife检验得分Fig.2 Jackknife test for climate factor contributions in the current potential distribution of Picea crassifolia图中All variables为全部环境变量; bio1: 年均温Annual mean temperature; bio2: 月均温度变幅Mean diurnal range; bio3: 等温性Isothermality; bio4: 温度季节性变化Temperature seasonality; bio5: 最暖月最高温Max temperature of warmest month; bio6: 最冷月最低温Min temperature of coldest month; bio7: 温度年变幅Temperature annual range; bio8: 最湿季均温Mean temperature of wettest quarter; bio9: 最干季均温Mean temperature of driest quarter; bio10: 最暖季均温Mean temperature of warmest quarter; bio11: 最冷季均温Mean temperature of coldest quarter; bio12: 年降雨Annual precipitation; bio13: 最湿月降雨Precipitation of wettest month; bio14: 最干月降雨Precipitation of driest month; bio15: 降雨季节性变化Precipitation seasonality; bio16: 最湿季降雨Precipitation of wettest quarter; bio17: 最干季降雨Precipitation of driest quarter; bio18: 最暖季降雨Precipitation of warmest quarter; bio19: 最冷季降雨Precipitation of coldest quarter

图3 主要气候因子的响应曲线Fig.3 Response curves of main climate factors

2.2 青海云杉未来潜在分布区变化

在3种情景、3种大气环流模型(GCMs)和两个时间段(2050s和2080s)下,通过模拟得到18个青海云杉的未来潜在分布模型,且所有模型的平均AUC测试值均高于0.99,表明模型的预测精度极高。对比潜在分布区面积变化发现(表2),尽管在3种增温情景下,青海云杉在2050s和2080s的总分布区面积与当前相比变化并不大(约1%),但不同适生等级的潜在分布面积变化显著。其中,0.1—0.3和0.3—0.5(低度和中度适生区)两个等级的分布受气候增温影响显著,在低度和高度增温情景(RCP2.6和RCP8.5)下这两个适生等级在2080s的分布区面积大幅下降(>30%),而中度增温情景(RCP4.5)下,低度和中度适生区面积则显著增加(>10%)。与此相对,青海云杉高度适生区的面积在所有情景下,与当前相比均呈减少趋势(表2)。从分布图上看,3种增温情景下,与当前相比,在2050s青海云杉在青藏高原中部及南部的潜在分布区(低度适生区)大幅减少,而增加了其在中国东部的潜在分布(低度和中度适生区),呈现出潜在分布区向东迁移的趋势；在2080s,青海云杉在青藏高原的潜在分布(低度适生区)进一步减少,而东部边缘地区(低度和中度适生区)依然呈扩张趋势(图4)。尽管如此,青海云杉的核心分布区域(高度适生区)仍然以青海东部、甘肃北部为主,无明显变迁趋势(图4)。

3 讨论

植物在区域尺度上的地理分布主要受气候因子的制约,其中水热条件起主导作用。本研究通过Maxent模型的Jackknife模块对影响青海云杉分布的气候因子进行评估后,证明了年最低降雨量是影响青海云杉分布的主要气候因子。与此不同,Xu等人利用最大熵模型模拟青海云杉在祁连山的分布后得出影响其分布的主要气候因子是最暖月最高温和最湿季均温[36]。这种结果的出现可能与预测分布尺度的差异有关：在本研究中,对青海云杉潜在分布的预测是基于全球尺度,远大于祁连山尺度。前人研究表明,青海云杉对气候因子的响应受海拔因素的影响,在森林上限(海拔3300 m)温度是影响青海云杉径向生长的主要因素,而在森林下线(海拔2770 m)其径向生长主要受降雨量的影响[37]。因此,我们推测在区域水平上,青海云杉的分布受到温度的影响,但其总体分布仍然受限于最低年降雨量。

Maxent模型预测结果显示,当前气候条件下,青海云杉的潜在分布区主要集中于：青海东部、甘肃东南部、宁夏大部分地区、西藏东部、四川西部山区以及陕西、新疆和内蒙古部分地区。此预测分布范围大于青海云杉的实际分布范围。原因可能来自于多方面：1)随着人类活动的加剧导致以青海云杉为优势树种的森林生态系统遭受了严重的破坏[38],仅1989年到1998年期间祁连山保护区的森林植被就减少了33.5%[39]。别强等[8]通过分布1960年至2000年期间祁连山保护区青海云杉林分布的动态变化发现,自20世纪60年代至本世纪初青海云杉面积显著减少,减少率最大时达到11.7%；2)青海云杉作为青藏高原东北边缘特有树种,主要集中分布于海拔2500—3200 m的范围[38],而本研究重点探讨气候因子与物种分布的关系,没有考虑分布海拔、坡向等地理空间因素,可能会造成预测的潜在分布区域大于实际已知生长区域的结果；3)Maxent模型基于气候相似性进行模拟,但并未考虑到物种现实分布时受到的其他因素,如物种的生长特性、扩散和迁移能力、自然更新能力、生物间的相互作用等[40-41]。尽管如此,利用Maxent模型得出的青海云杉潜在分布区域均具有相似的气候特点,说明青海云杉林具有很高的恢复潜力[42],可为上述地区进行森林生态系统修复与恢复提供一定的树种参考价值。

气候变化背景下,尽管青海云杉的总潜在分布面积波动较小(约1%),但其在不同适生等级的潜在分布面积及分布范围均发生了改变。这种差异的产生可能与本研究中对潜在分布的等级划分有关,总适生面积是通过Maxent第10分位训练存在阈值生成的二元图计算而来,而适生等级则是人为划分,与前者的适生范围存在一定差异。在不同适生等级,虽然中度及以下的适生区在未来3种增温情景下两个时间段(2050s和2080s)的面积变化不一致,但高度适生区的面积均呈下降趋势,这表明未来气候增温仍然影响了青海云杉的分布。前人对祁连山青海云杉过去近30年对气候响应的调查分析发现,青海云杉在全球变暖背景下增强了对温度的敏感性,这种敏感性与区域内水分含量密切相关,干旱会增加青海云杉与温度负相关的敏感性[43]。且有研究表明青海云杉的分布地未来气候将向温暖干燥转化[14],这也间接支持了与本研究的结果,即青海云杉的核心分布范围呈下降趋势。梁和延对我国21世纪末期气候变化预测发现,在未来气候显著变暖情景(RCP8.5)下,中西部地区的夏季降雨量明显下降[33]。这也解释了本研究中青海云杉在该情境下本世纪末期(2080s)在中国中西部的分布范围显著低于其他情景及时间段的现象(图4)。此外,在不同增温情境下,不同适生等级的青海云杉适生区在2050s和2080s的潜在分布区面积变化趋势并不一致,这可能与大气环流模型存在地区气候变幅差异有关[44-45]。

表2 青海云杉不同适生等级未来潜在分布与当前相比面积变化百分比/%

Table 2 The percentage change in potential environmental area forPiceacrassifoliain the future when compared with the current at different level of fitness

适生等级Level of fitnessRCP 2.6RCP 4.5RCP 8.52050s2080s2050s2080s2050s2080s<0.1-8.10+26.62-19.23-9.46-9.24+43.560.1—0.3+5.83-13.81+2.20+10.19-1.83-22.500.3—0.5+20.54-32.14+13.22+20.764.74-37.890.5—0.8-9.48-16.98-21.73-3.72-12.80-6.25总面积变化 Percentage change in total area-0.87-0.71-0.53-0.77-1.13-0.34

“+”和“-”分别指示“增加”和“减少”

图4 在3种增温情景下,青海云杉当前及未来不同适生等级的潜在分布区变化Fig.4 The change of potential distribution for Picea crassifolia in the future when compared with the current at different fitness under three greenhouse gas emission scenarios

4 结论

本研究基于最大熵模型证明了年最低降雨量是影响青海云杉分布的主要气候因子；在未来气候变化情景下,青海云杉分布范围虽然有减小趋势,但其高度适生区仍以祁连山、贺兰山区为主。因此在未来西部高山地区,特别是作为我国重要生态屏障的祁连山、贺兰山区,青海云杉依然具有重要的经济价值并将持续其生态服务功能。