潘少安,李旭华,冯秋红,刘兴良,孙建新
1 北京林业大学生态与自然保护学院,北京 100083 2 四川省林业科学研究院,森林和湿地生态恢复与保育四川省重点实验室,成都 610081 3 四川卧龙森林生态系统国家定位观测研究站,成都 610081
多年前,科学家们就已证实,全球气候变化已经无可争辩[1]。当前的气候变化趋势和气候预测均表现出温度上升随气候变化的一致趋势[2]。陆地生态系统对温度变化高度敏感,且温室气体的增加使陆地生态系统面临发生重大转变的危险[3]。近几十年来,气候变暖导致亚高山针叶林分布向更高海拔转变[4],而水分的减少限制了亚高山树木的分布范围,同时还可能导致林木数量的减少[5]。有研究表明峨眉冷杉在海拔3000 m处分布面积逐渐减少最终消失并被其它物种取代,而在3600 m处以分散分布出现,最终成为优势物种[6]。亚高山森林分布的动态变化对气候变化具有十分重要的指示作用[7],尤其是冷杉属物种,对气候变化十分敏感,对温湿度稳定性有着严格要求[8],因此,研究岷江冷杉对气候的响应对于理解全球变化具有极其重要的意义。
岷江冷杉(Abiesfaxoniana)为松科冷杉属常绿乔木,是我国特有树种,也是川西亚高山地区森林的优势树种之一,通常在海拔2800 m以上的区域形成纯林,或与红杉、方枝柏等形成混交林[9],是岷江上游地区阴坡林线的主要建群树种[10]。岷江冷杉林是四川省面积最大的原始森林类型和大熊猫的重要栖息场所,在涵养水源、保持水土方面具有十分重要的功能,亦维持着大熊猫栖息地生态系统的稳定与安全[11—12]。目前关于岷江冷杉的研究较为丰富,涉及岷江冷杉林的群落动态变化[13]、凋落物分解[14]、年轮-气候响应研究[15—17]、林地水源涵养[18]以及天然更新[19—20]等多个方面。在全球气候变化背景下,岷江冷杉种群的潜在分布格局及其对未来气候变化的响应方面的研究却鲜有报道。
近年来,研究表明森林群落随着气候变暖而向上迁移,但很少提及这些变化的具体方式和具体方法,而植物在应对气候变化的过程中,迁移或扩散能力是植物生存的先决条件[21]。我国学者在对目标物种的潜在分布区预测,以及物种分布对气候变化的响应等方面已开展了诸多研究[22],而物种分布模型即为最常用的研究方法,该方法的基本原理是利用观测到的物种分布数据及其对应的环境变量信息,关联性构建二者的相关模型,进而确定目标物种的实际分布并对物种的潜在分布进行预测[23]。目前常用的物种分布模型有MaxEnt模型、Garp模型、Bioclim模型、Dumain模型、广义线性模型、广义可加模型等[24],但MaxEnt模型以其运行样本量要求低、易操作、运行时间短、数据处理能力强、预算精度高等优点[25—26]而得到广泛应用。因此,本研究采用物种分布模型—最大熵模型(MaxEnt)和ArcGIS软件,对岷江冷杉在当代气候条件和未来气候变化情景下的潜在分布区和生境适宜性进行评估和分析,并确定影响岷江冷杉空间分布格局的主要影响因子和分布阈值,分析岷江冷杉对未来气候变化的响应策略,为川西山地岷江冷杉天然更新、大熊猫栖息地恢复、退化森林质量精准提升等提供理论依据。
岷江冷杉的地理分布数据主要从3个方面获取:包括岷江冷杉样地调查数据、公开发表文献资料、权威的标本信息网络数据库,如全球生物多样性信息网GBIF(https://www.gbif.org),中国数字植物标本馆等(http://www.cvh.org.cn)。对搜集到的物种分布数据进行整理,剔除重复的和研究区域范围外的坐标点,共得到93个参与模型运算的物种分布点,并按要求将其另存为csv格式。
本文所用的环境数据包括气候数据、土壤数据和地形数据,其中气候数据来源于世界气候数据库WorldClim(https://www.worldclim.org/),分别选取当代(1970—2000年)气候数据和未来两个时期(2050s和2070s)3种气候变化情景RCP2.6,RCP4.5和RCP8.5的生物气候变量,3种气候情景分别表示未来时期CO2排放浓度低、中、高的3种可能变化路径[27],每个时期的气候数据变量均为19个。土壤数据从寒区旱区科学数据中心提供的世界土壤数据库(Harmonized World Soil Database,HWSD)中提取出研究区域与岷江冷杉生长相关的土壤因子,如表层土壤容重、pH值、沙粒含量、粉粒含量、粘土含量、碎石百分比、有机碳含量、基本饱和度、碳酸盐含量等。地形数据包括海拔、坡度和坡向,其中海拔数据(空间分辨率1 km)从全球海拔数据库(https://globalmaps.github.io/)中下载并提取,坡度和坡向则分别根据海拔数据运用ArcGIS中的空间分析工具计算得到。对上述所有环境因子进行研究区掩膜提取和空间处理,使其具有统一的空间坐标系(GCS-WGS84)和相同的空间分辨率(30 arc sec×30 arc sec,约1 km)。
研究指出环境因子间的多重共线性会导致物种分布模型的过度拟合[28—29],因此在运行模型前,对环境因子分别进行Pearson相关性分析,进而筛选出相关性较低且对物种分布有重要影响的环境因子(相关系数|r|<0.8)来参与模型运算。最终共选取了16个环境因子参与岷江冷杉分布区预测与生境评价,包括6个气候因子,7个土壤因子和3个地形因子。具体指标情况见表1。
表1 用于MaxEnt模型的环境变量
将符合MaxEnt模型格式要求的岷江冷杉物种分布数据,分别与当前和未来两个时期不同气候情景下的环境变量数据导入模型中,随机选取25%的分布点作为测试集,75%的物种分布点作为训练集,模拟共得到7个岷江冷杉适生分布区的预测结果。模型预测结果的精确度则通过受试者工作特征曲线(Receiver Operating Characteristic curve,ROC曲线)下面积(Area Under Curve,AUC)进行评价[30],该值大小与模型的预测精度呈正相关[31],当AUC值小于0.6时模型预测失败,当AUC值介于0.8—0.9表明模型预测精度较好,达到0.9—1.0时则模型预测精度很高[32]。同时,模型通过贡献率、置换重要值和刀切法检验(Jackknife test)3个方面来综合评估个各环境变量在影响岷江冷杉地理分布中的重要性。
模型的预测结果是通过物种在待预测地区的存在概率P(取值0—1)来表示物种的分布适宜性[33],将模型运行得到的分布区预测结果进行适宜性等级划分和可视化表达,结合岷江冷杉的生态学特点,采用重分类中的人工(Manual)分级方法,将岷江冷杉的潜在分布区划分为4个等级,分别是:适宜性指数在0—0.1为非适生区,0.1—0.3为低适生区,0.3—0.5为中适生区,0.5—1.0为高适生区。为了进一步分析岷江冷杉在未来气候变化情景下适宜分布区的空间变化格局,对比当代适生区分布结果,我们分别提取了新增适生区、丧失适生区、保留适生区和非适生区4种岷江冷杉适生区的变化类型。
图1 岷江冷杉分布区预测的受试者工作曲线 Fig.1 Receiver operating characteristic (ROC) of Abies faxoniana AUC: 曲线下面积 Area under curve
模型运行结果显示,训练集和测试集的AUC值分别为0.955和0.853(图1),达到了模型预测精度要求,由此表明该模型对岷江冷杉潜在分布区预测结果的可信度较高。
对预测的岷江冷杉潜在适生区等级划分的结果如图2所示,经空间分析模块统计各适生区的分布面积,结果表明岷江冷杉中适生分布区和高适生分布区的面积约为33176.86 km2和31416.84 km2,分别占全省总面积的5.91%和5.60%,低适生区占全省总面积的12.30%,而非适生区的占比最高为76.19%。岷江冷杉的中、高适生区集中分布于阿坝藏族羌族自治州汶川县、理县、黑水县、松潘县等地;雅安市的天全县、宝兴县、绵阳市的平武县、北川县等地,除此之外在成都市大邑县、甘孜藏族自治州的丹巴县、康定县等地的零星的面积也是岷江冷杉的潜在适生分布区。
图2 当代气候条件下岷江冷杉适生等级地理分布格局Fig.2 Geographic distribution pattern of suitable grade for Abies faxoniana in contemporary China
在影响岷江冷杉地理分布的环境因子中,气候因子的累计贡献率为74.16%,地形因子和土壤因子的累计贡献率分别为13.33%和12.51%,而这三类因子对应的置换重要性则分别为84.79%,12.25%和2.96%(表2),由此表明岷江冷杉地理分布主要受气候因子的影响,而地形对岷江冷杉地理分布的影响次之。
具体来讲,对岷江冷杉地理分布贡献率较高的前4个环境变量分别是降水季节性变异系数、气温年变化幅度、昼夜温差月均值和海拔,其累计贡献率为80.77%;置换重要值较高的环境因子分别是气温年变化幅度、年降水量、海拔和降水季节性变异系数,其累计值达到了84.99%(表2)。此外,根据刀切法检验结果,在仅使用单变量时的正规化训练增益最高的3个变量分别是气温年变化幅度、年降水量、海拔,而在使用除单个环境变量外的其他变量时,气温年变化幅度和降水季节性变异系数的增益值降低幅度最大(表2),由此可见,这二者具有其他环境变量中没有的有效信息。综合上述3类评估结果,影响岷江冷杉地理分布的气候因子主要是降水季节性变异系数、气温年变化幅度和年降水量,地形因子则主要是海拔。
根据环境因子响应曲线,计算影响岷江冷杉地理分布的各主导因子的阈值(即物种存在概率>0.5的范围)。结果显示(图3),岷江冷杉适宜生长的海拔范围是2310—3757 m;气温年变化幅度的适宜范围为29.3—32.5 ℃;降水季节性变异系数的适宜范围为76.5—85.4,最适宜为81.7左右;而适宜岷江冷杉生长的年降水量的范围是694.1—791.7 mm。
表2 影响岷江冷杉分布的环境因子贡献率
图3 主要环境因子的响应曲线Fig.3 Response curves of four major environmental factors
在未来两个时期(2050s和2070s)的不同气候变化情景下,岷江冷杉潜在中适生区和高适生区的面积较当代气候条件下适生区面积均有所增加,且潜在高适生区面积的增加幅度总体上高于中适生区面积的增加幅度。此外,未来两个时期的中适生区的面积占比均随着排放浓度的增加而逐渐减少,而高适生区的面积占比则表现出先增加后降低的变化(表3),即在中等排放浓度时的潜在高适生区面积最高。
表3 不同气候模式下岷江冷杉潜在分布区的面积占比/%
岷江冷杉适生分布区的空间变化格局显示,在未来不同气候变化情景下,岷江冷杉的潜在适生分布区总体向南移,新增适生区多分布于盆周山地,而丧失适生区集中分布于原适生区北缘和西北缘(图4)。其中在2050s三种排放情景下,岷江冷杉潜在适生区的新增率在RCP8.5排放情景下最高,为34.27%,其次是RCP2.6排放情景下的34.20%(表4);而在2070s,岷江冷杉新增适生区面积在RCP8.5排放情景下最少,新增率为31.62%,在RCP4.5排放情景下的新增适生面积最大,新增率为39.79%。岷江冷杉潜在适生区在未来两个时期不同气候变化情景下,丧失面积最大的是2070s的RCP8.5排放情景,丧失面积最小的是在2070s的RCP4.5排放情景。
表4 未来气候变化情景下岷江冷杉适生区的变化率/%
图4 不同气候变化情景下岷江冷杉适生区空间变化格局Fig.4 Spatial variation patterns of suitable areas of Abies faxoniana under different climate change scenariosRCP: 典型浓度路径 representative concentration pathway
当前,利用物种分布模型来预测物种的地理分布已开展了很多研究,在众多物种分布模型中,最大熵模型(MaxEnt)是使用最为广泛且具有较好预测能力的模型[22]。本研究基于MaxEnt模型和ArcGIS空间分析,定量展示了当代气候变化条件下岷江冷杉在四川省境内的适生分布区,预测结果显示岷江冷杉潜在分布区要比实际已知的生长区域范围大,主要分布于阿坝藏族羌族自治州、雅安市中部和北部以及甘孜藏族自治州丹巴县、康定县等地的零星地区,这些地区所在的森林分区为川西亚高山云冷杉林区,位于青藏高原东南缘,也是金沙江、大渡河、岷江、涪江等大江大河主要的水源涵养区或源头。适生分布区的气候受高原地形的决定性影响,形成冬寒夏凉的高山气候特点,这也符合岷江冷杉耐阴性强,喜冷湿气候的生物学特性。研究结果表明,当代气候条件下,岷江冷杉潜在适生区(中适生区和高适生区)的面积占研究区总面积的11.51%,说明岷江冷杉的生态位较为狭窄。此外,模型模拟结果经ROC曲线精度检验,AUC值高于0.8,且预测的潜在分布区与实际调查分布基本一致,由此表明了MaxEnt模型对岷江冷杉的预测精度较好,较为真实地反映了岷江冷杉潜在适生区的分布情况,但模拟的精度尚未达到很高,这可能与物种分布数据量不足有关,模型的表现与物种分布数据的数量正相关[34]。
在不同的环境尺度下, 影响物种分布的各类因素的作用程度是不同的。Soberon和Peterson[35]将影响物种分布的因素总结为4类:(1)非生物因素,包括气候、土壤条件、地形、水文等;(2)生物因素,包括种间竞争、取食、互利共生等;(3)该地物种自身迁移能力和地理区域的特性;(4)物种或种群对新环境的适应能力。通常,非生物因素主要在大尺度空间影响物种的分布,这些因素很大程度上决定了物种的分布范围和格局,包括生理制约、物种对气候和生境梯度的响应和选择等。刘晓彤等[22]文中指出,统计自2000年以来有关物种分布模拟的研究中,有50%以上的研究只考虑了气候因素,涉及地形和土壤及其他因子的研究较少,本研究综合考虑了气候、地形及土壤因子对岷江冷杉潜在分布区的影响。结果表明,气候因子对岷江冷杉地理分布预测的贡献最大,土壤因子和地形对分布预测的贡献率较小。后二者可能在小尺度范围内影响植物的分布,但在区域尺度上的影响较小,这可能是由于地形和土壤因子对植物的影响是间接的[36]。影响岷江冷杉分布的降水因子主要是降水季节性变异系数,这一气候因子反映了岷江冷杉在不同生长期对降水的需求不同,研究表明降水量的季节性变异系数在76.5—85.4的区域为岷江冷杉的适宜生长范围。影响岷江冷杉分布的温度因子是年均气温变化幅度,其适宜范围是29.3—32.5 ℃,即在年均气温变幅较高的地区适合岷江冷杉分布。上述两个气候因子的主导性反映了气候条件的变异性决定了岷江冷杉的分布格局。在地形因子中岷江冷杉受海拔的影响最大,在海拔2310—3757 m区间是岷江冷杉的适宜生长范围,而《四川松杉植物植物地理》[37]指出岷江冷杉垂直分布于海拔2400—4000 m,坡度>25°的阴坡、半阴坡的高山地带,本文结果符合岷江冷杉的实际生长环境;虽然海拔适宜范围的最低阈值低于2400 m,说明在海拔2300—2400 m的亚高山地带也是岷江冷杉的潜在适生分布区,在未来气候变化或人工种植的情况下岷江冷杉是可以存活的。本研究中参与模型构建的环境数据只涉及了非生物因素,未考虑岷江冷杉自身适应能力、种间相互作用及人类干扰数据,因而预测结果与实际情况也会存在一定的偏差。
在全球气候变化背景下,预测物种在未来气候变化下的分布变化趋势,对物种多样性保护具有重要意义[27]。岷江冷杉对气候变化的响应表现为:在未来两个时期的不同气候变化情景下,其潜在适生区面积均较当代气候条件下的适生区面积有所增加,但中等适生区的面积占比均随着排放浓度(低、中、高)的增加而逐渐减少,而潜在高适生区面积在中等浓度时最高,表明不同生境中的岷江冷杉对不同CO2浓度以及温度的变化产生不同的响应,在同一个时期内随着CO2浓度的升高,一定温度范围内,温度升高很大程度上将抑制低海拔岷江冷杉的生长,从而导致岷江冷杉部分生境丧失;但在高浓度排放情景下,高海拔地区的温度升高明显,使得岷江冷杉有向高海拔低温区扩张的可能[38],高海拔地区升高的温度有利于岷江冷杉形成层更早的活动,表明高海拔地区将由不适宜生境逐渐向适宜生境转变;在未来更长的时间内,中等CO2浓度排放情景更有利于岷江冷杉林生境的扩张,低CO2浓度排放情景和高CO2浓度排放情景下,升温幅度低于2 ℃和升温幅度太高都对岷江冷杉的生长不利。
未来气候变化情景下岷江冷杉的潜在适生分布区总体向南迁移,这与Liao[8]等人的研究结果一致,该研究表明在未来气候变化情景下,岷江冷杉的潜在适生区会沿着岷江流域向南扩张至云南省的东北部。综合各个排放情景,生境丧失区域主要分布在阿坝藏族羌族自治州的若尔盖县北部、阿坝县,以及雅安市、乐山市、眉山市以及凉山彝族自治州交界的马尔康县、天全县、荥经县、泸定县、石棉县等。根据生境适宜分布区预测结果(图2),这些分布区主要位于非适生区和低适生区的交界处,处于岷江冷杉分布的林线下缘。全球变暖背景下,中低海拔岷江冷杉存在季节性水分胁迫,抑制岷江冷杉生长,尤其是对大龄树木的影响更加强烈[15]。此外,物种生长过程中的局部环境的生存竞争压力也发挥着不可忽视的作用。高海拔地区气温逐渐升高,原本不利生境条件得到改善,林线下缘岷江冷杉面临更大的资源竞争压力[39]。较大的生存压力对岷江冷杉种群幼苗的存活极为不利,以至于岷江冷杉种群幼苗死亡率较高,不利于种群的更新与扩张[10]。四川省地形复杂多样,海拔垂直梯度较大,复杂的地形地貌条件导致局部气候变化多样,四川省生物气候类型未来转变趋势由温度低的类型像温度高的类型转变[40],低海拔一些原本不适宜岷江冷杉分布的区域,可能转化成潜在适宜分布区。
对于物种的新增适生区,表明未来气候变化使得这些区域的气候条件更加适宜该物种的生长,物种分布的可能性增加,但实际的分布格局还要取决于物种向新增分布区迁移的可能性[41]。未来气候变化条件下的新增适生区能够为岷江冷杉未来的更新与恢复提供指引,可以通过人为的方式在这些区域引种或培育岷江冷杉,增加其迁入的可能性。而对于丧失适生区,可以通过采取迁地保护、提高岷江冷杉幼苗繁育等措施来提高其应对气候变化的能力。未来气候变化条件下,保留适生区是物种应对气候变化的安全地和避难所[42],因而更应该注重和加强该区域岷江冷杉的保护。可以看到,未来不同气候变化情景下岷江冷杉适生区的保留率基本保持在80%左右,说明岷江冷杉应对气候变化的能力较强。
(1)最大熵模型MaxEnt能够较好地预测岷江冷杉的潜在适生分布(AUC值>0.8)。在当代气候条件下,岷江冷杉的中、高适生区集中分布于阿坝藏族羌族自治州、雅安市中部和北部、绵阳市的平武县、北川县等地,分别占全省总面积的5.91%和5.60%。
(2)影响岷江冷杉地理分布的气候因子主要是降水季节性变异系数、气温年变化幅度、年降水量、海拔和土壤碎石体积百分比。其中气候因子占主导作用,地形因子次之。适宜岷江冷杉分布的海拔范围2310 —3757 m,气温年变化幅度的范围为29.3—32.5 ℃;降水季节性变异系数的适宜范围为76.5—85.4,以及年降水量的适宜范围为694.1—791.7 mm。
(3)在未来气候变化情景下,岷江冷杉的适生分布区南移,新增适生区多分布于盆周山地,而丧失适生区则主要分布于原适生区北缘和西北缘。其中在岷江冷杉潜在适生区在未来两个时期不同气候变化情景下,丧失面积最大的是2070s的RCP8.5排放情景,丧失面积最小的是在2070s的RCP4.5排放情景。