基于SSPs气候场景的濒危植物银杉潜在分布区预测

罗楚滢,佘济云*,唐子朝

(1.中南林业科技大学林学院,湖南 长沙 410004;2.湖南省农林工业勘察设计研究总院,湖南 长沙 410007)

气候变暖影响着全球水热分布与土壤微生物活动,而植物分布格局主要受光照、水分、土壤、地形等因子制约[1-2]。因此,在全球气候变化的条件下,对植物潜在分布的预测已成为林学、生态学等相关学科的研究热点[3]。世界气候研究计划(WCRP)耦合模拟工作组(WGCM)组织了第6次国际耦合模式比较计划(CMIP6)。相较于CMIP5中的气候模式,引入了共享社会经济路径(shared socio-economic pathways, SSPs)的CMIP6,能够实现对各种减排措施成本效益的评估,同时提供了比RCPs更多元化的空气污染物排放情景,这些改进弥补了CMIP5中RCPs情景对社会经济活动考虑的不足[4-5]。SSPs包含了5种经济路径,可用于预测未来不同气候政策下温室气体的排放情景,评估政策和发展路径对植物生长与繁殖所带来的影响,模拟气候变化对气象灾害、农业生产、生态系统等领域所带来的影响[6-7]。

银杉(Cathayaargyrophylla)属松科银杉属,起源古老,是第三纪孑遗的珍稀树种,为我国一级濒危保护植物,具有生长慢、成年树高大且树形优美的特征[8]。目前,银杉仅在我国亚热带山地零星分布。已有相关学者对银杉的苗木生长特性、濒危因素、潜在分布区预测等进行了一系列的研究[9-10],运用SSPs气候情景预测银杉未来的分布变化趋势鲜有报道。

目前,预测与分析物种分布范围、变化情况的模型运用较多的有:生物气候分析和预测系统模型(BIOCLIM)、生态位比较模型(CLIMEX)、最大熵模型(MaxEnt)[11-13]等,其中MaxEnt模型已经被广泛运用于物种空间分布的预测、全球气候变化对物种分布的影响评估、珍稀濒危物种保护[14-17]等方面。如韩淑敏等[18]基于MaxEnt模型预测水土保持、防风固沙树种白榆(Cerasusconradinae)的潜在分布格局及未来变化,为白榆的科学保护和合理开发提供理论依据;周炳江等[19]基于MaxEnt模型对国家二级保护树种云南榧树(Torreyayunnanensis)进行潜在生境分析,为云南榧树的保护区建立、引种栽培等提供理论依据。本研究基于MaxEnt模型对银杉的潜在分布格局以及未来变化进行预测分析,以期为银杉的科学保护和合理开发提供理论依据。

1 材料与方法

1.1 银杉分布数据的收集及处理

银杉分布数据来自中国数字植物标本馆(CVH,www.cvh.ac)、中国植物图像库(PPBC,https://ppbc.iplant.cn/)、全球生物多样性信息网络(GBIF,www.gbif.org)等信息平台及文献检索。通过软件图新地球4(LocaSpace Viewer, LSV)获取相应的经纬度坐标,并通过坐标反查法验证其准确性。将分布数据加载进ArcGIS 10.8,为排除空间自相关性较大造成过拟合,根据环境因子数据的空间分辨率(1 km × 1 km),设置缓冲区半径为1 km,每个缓冲区内仅保留1个样本点[20]。去除重复、不精确及人工栽培样本,最后获得参与建模样本共41个。

1.2 环境因子的获取及筛选

根据相关学者的研究结果,筛选对银杉分布影响较为显著的环境因子,选用了4种环境变量,包括气候、土壤、地形和辐射[21-22],共计31个环境因子(表1)。其中气候数据(bio1—bio19)来源于世界气候数据库(Worldclim,www.worldclim.org),气候数据基于1970—2000年世界各地气象站所收集到气候因子,采用插值法生成的全球气候栅格数据,其空间分辨率为30′(1 km × 1 km)[23]。未来3个时期(2050s、2070s、2090s)气候数据选用共享社会经济路径(shared socio-economic pathways, SSPs)情景下的3种不同发展路径,分别为可持续发展路径(SSP1-2.6)、中间发展路径(SSP2-4.5)、以化石燃料为主发展路径(SSP5-8.5)[4]。地形因子包括海拔、坡度和坡向,其中高程数据来自地理空间数据云(DEM,www.gscloud.cn)。辐射数据来自第5代ECMWF全球气候和天气再分析资料ERA5(https://cds.climate.copernicus.eu),包含了逐小时大气、海浪和路面数据。选用的8种土壤因子(土壤参考深度、顶层土壤质地、淤泥含量、沙含量、黏土含量、酸碱度、土壤有效含水量、有机碳含量)数据均来源于世界土壤数据库(HWSD,www.fao.org)。

表1 环境变量及基于最大熵模型的各环境变量贡献率

在构建初始MaxEnt模型时选取上述31个环境因子进行建模,找到对预测结果贡献率较小或没有贡献的环境因子,同时运用SPSS软件中的Pearson分析,找到自相关性过高的环境因子[24],将二者实验结果相结合,剔除自相关性过高且贡献率较低的环境因子,最终得到17个环境因子参与建模(表1)。

1.3 银杉适宜分布区预测模型构建与评价

当前气候情景下,银杉潜在分布区预测方法如下:将筛选后的分布数据与环境因子导入MaxEnt模型,建模时采用MaxEnt 3.4.4版本;设置随机选取25%样本点作为测试数据集,75%样本点作为训练数据集[25];通过刀切法来检验环境因子的权重,最大迭代次数设为500;选取AUC值来评估模型精度[AUC(area under the ROC curve)的取值范围在[0,1],越接近1就说明预测结果越好,其模型预测结果就越准确]。具体评估标准为:AUC值为(0.5,0.6]时预测结果无效,(0.6～0.7]时预测结果较差;(0.7～0.8]时预测结果一般,(0.8～0.9]时预测结果良好,(0.9～1.0]时预测结果极好[26]。AUC值最大的一组作为最终预测结果[27]。同样,采用上述方法预测银杉在2050s (2041—2060年)、2070s (2061—2080年)、2090s (2081—2100年)3个不同时期的潜在地理分布。

1.4 银杉适生生境等级划分与地理空间变化预测

将MaxEnt输出的预测结果导入ArcGIS,转换为栅格格式后进行重分类,ArcGIS提供了多种重分类方法,参考相关学者等级划分的方法[28],经尝试发现手动分类(manual)、分位数分类(quantile)、几何间隔分类(geometrical interval)、自然断点分类(natural breaks)、标准差(standard deviation)等方法均会导致生态适宜区过大或过小,这与银杉在我国的分布现状严重不符。最终采用最为合适的平均间距法(equal interval),将其按适宜性指数划分:[0.00,0.190 593]时为非适生区;(0.190 593,0.381 186]、(0.381 186,0.571 779]时为一般适生区;(0.571 779,0.762 372]、(0.762 372,0.952 965]时为核心适生区。得到结果后计算各适宜生境等级的面积。未来银杉潜在分布区适宜性指数的划分也采用相同的等级,再通过3D Analyst工具对未来3个不同时期的重分类栅格数据进行加减法计算,得到银杉未来不同时期分布的缩减区、稳定区和扩张区,以探讨其分布的总体变化趋势。

2 结果与分析

2.1 当前环境条件下银杉在中国的潜在分布

银杉潜在分布预测模型训练集数据和测试集数据的AUC值分别为0.971和0.996,表明模型预测结果可靠,可信度较高。由分析结果可知,银杉潜在分布范围为92°57′～121°54′E,21°34′～32°51′N,主要为我国的南方地区(图1)。其中核心适生区面积约为1.325×105km2,约占我国国土面积的1.38%,主要分布于我国湖南省中部娄底市、南部永州市、郴州市、西南部邵阳市,湖北省西南部恩施土家族苗族自治州,广西壮族自治区北部桂林市、柳州市,台湾东部宜兰县、花莲县,重庆市南部和东南部,贵州省北部的遵义市、铜仁市、东南部的黔东南苗族侗族自治州。此外,在四川省东南部泸州市,江西省西部九江市、吉安市、东部上饶市、抚州市,安徽省南部黄山市,江苏省南部无锡市,浙江省西南部丽水市,福建省西北部南平市亦有少量核心适生区。一般适生区面积约7.443×105km2,约占研究区面积的7.73%,主要分布在我国云南省中部、西部、南部、东部,四川省东南部,贵州省全域,重庆市全域,湖南省全域,湖北省西南部,广西壮族自治区北部,广东省北部,江西省部分地区,福建省北部、西部,浙江省西南部、东北部,江苏省南部,安徽省南部,台湾东部。非适生区面积约为8.757 2×106km2,约占研究区面积的90.89%。银杉在研究区内的适宜分布面积较小,尤其是最适合银杉生长的核心适生区,这与目前银杉在我国的实际分布相符。

底图审图号为GS(2020)4630号。下同。The same below.图1 当前银杉在我国的潜在分布区Fig. 1 Potential distribution area of Cathaya argyrophylla in China under current climatic condition

2.2 银杉适生区环境因子分析

由最大熵模型中各变量的贡献率可知,前4名的因子累计贡献率为91.9%,其包括:最干月份降水量(66.1%)、下行紫外线辐射(14.6%)、海拔(9.1%)、最冷月最低温(2.1%)。由此可见,上述4个环境因子对银杉的潜在分布影响较大;而最干月份降水量、下行紫外线辐射可以被认定是影响银杉潜在分布的主要因子。

环境因子响应曲线能够反映环境变量与物种存在概率之间的相关性。若以生存概率大于0.5作为适应范围,则适宜银杉生存的主要环境条件为:最干月份降水量18.03～215.63 mm,下行紫外线辐射1 070 728～1 437 806 W/m2,海拔493.68～1 731.10 m,最冷月最低温-1.01～4.05 ℃(图2)。

图2 主要环境因子响应曲线Fig. 2 Response curves of major environmental factors

2.3 未来气候变化下银杉的核心适生区预测

由模型预测结果可知,银杉在不同年份和发展路径下的核心适生区均有所增加,银杉核心适生区面积增加量与发展路径有关。虽然银杉核心适生区总体呈增加趋势(表2),但核心适生区面积占比仍然很小,分布区主要以黔、湘、桂、渝为中心(图3)。银杉未来变化趋势:扩张区主要分布于云南、浙江、福建、贵州四省;保留区主要分布于湖南省与广西壮族自治区交界一带、贵州省与重庆市交界一带、湖北省西南部、贵州省中部地区;缩减区主要分布于湖南省中部、贵州省中部地区。

SSP1-2.6、SSP2-4.5、SSP5-8.5分别指共享社会经济路径下可持续发展路径、中间发展路径和以化石燃料为主的发展路径。SSP1-2.6,SSP2-4.5,SSP5-8.5 respectively refered to SSPs scenario, namely the sustainable development path, the inter-mediate development path, and the fossil flue based development path.图3 不同气候情景下银杉核心适生区地理分布变化Fig. 3 Changes of geographic distribution of core suitable areas of C. argyrophylla under different climate scenarios

在SSP1-2.6情景下,银杉核心适生区面积表现出低扩张率、低缩减率,总体面积较为稳定,在2090年大量增加。相比当前核心适生区,2050s整体变化量为+1.20×103km2(表2);2070s整体变化量为+1.93×104km2;2090s整体变化量为+7.49×104km2。在SSP2-4.5情景下,银杉表现出高扩张率、低缩减率,但扩张面积逐渐减少。2050s整体变化量为+15.95×104km2,2070s整体变化量为+12.83×104km2,2090s整体变化量为+9.72×104km2。在SSP5-8.5情景下,银杉表现出较高扩张率,低缩减率。2050s整体变化量为+5.27×104km2,2070s整体变化量为+5.03×104km2,在2090s整体变化量为+6.50×104km2。

3 讨 论

3.1 预测结果合理性与差异性分析

本研究预测银杉核心适生区主要分布在重庆东南境内的金佛山及周边地区、湘桂二省交界处的越城岭腹地、湖南省境内雪峰山一带、贵州北部大娄山和东南部苗岭地区、台湾地区内雪山山脉与中央山脉北部,银杉潜在分布区总面积较小,且核心分布区较为狭窄,这与以往的研究结果较为符合[22],说明本次预测结果具有一定的可靠性。同时,在环境因子选取方面,冉巧等[22]以19个生物气候因子作为影响因子进行研究,发现最干月份降水量是影响银杉潜在分布的主要影响因子,贡献率达到58.6%。本研究中,最干月份降水量贡献率居首位,达到66.1%。

此外,本研究与前人的研究结果也存在一定的差异。刘倩等[21]选用3种辐射因子、4种气温降水因子、2种生物气候因子、5种土壤因子作为影响因子进行研究,结果指出,太阳辐射是影响其潜在分布最主要的环境因子,贡献率达到67.8%。而本研究在19种生物气候因子的基础上,结合紫外线辐射因子(1种)、土壤因子(8种)、高程因子(3种),共31种环境因子,更加全面系统地参与模型构建与结果分析。为了避免高度自相关的变量在模型预测过程中干扰预测结果,使用SPSS软件进行Pearson分析并与第1次MaxEnt模型预测结果相结合,剔除自相关性过高的环境因子,减少了冗余信息对模拟结果的影响,提高了预测结果的准确性。

在银杉未来潜在分布区预测方面,相关学者研究采用IPCCAR5定义的4种不同浓度路径(RCPs)[22],而本研究采用政府间气候变化专门委员会于2010年颁布的共享社会经济路径(SSPs),该情景是在典型浓度路径(RCPs)情景基础上发展而来,弥补了RCPs情景中的考虑不足,预测结果更为合理[6]。

3.2 关键环境因子对银杉分布的影响

银杉集中分布于我国南方山地地区,该地区受亚热带季风气候影响,夏季高温多雨,冬季温暖无严寒,光照充足[29]。这与本研究对银杉适生区环境因子分析结论相符。

研究结果显示,排名前4名的因子分别是最干月份降水量、下行紫外线辐射、最冷月最低温、海拔。最干月份降水量与太阳辐射因子对银杉分布的影响更大。当最干月份降水量达到196.313 6 mm以上时,银杉达到最大分布概率。充足的水分对银杉种群发生至关重要,是影响银杉分布的最关键环境因子;当最冷月最低温在1.46 ℃左右时,银杉分布概率也较高。研究结果符合前人报道的银杉适宜生长在冬季温和、降水量充沛的潮湿区域结论[22]。同时,下行紫外线辐射也是限制银杉地理分布的主要因素之一。紫外线辐射对植物生长发育具有重要影响,其双重性在分子、细胞、器官、植物个体甚至生态系统水平上影响着植物的生长发育,干预植物生理生化的代谢过程及相关基因的表达与调控[30],紫外线辐射强度可能通过影响银杉繁殖与生长来限制其分布范围,张旺锋等[31]的研究显示,适度遮阴有利于银杉抵御冬季的光抑制。本研究预测结果表明海拔对银杉也存在一定的影响,根据前人采集的标本以及实地调查发现,天然银杉种群大多分布于海拔800 m以上的山脉地区,例如沙角洞自然保护区银杉分别位于海拔1 050 m与950 m处[32],高海拔地区温度较低,降水丰沛且人为干扰较少,更适宜银杉喜雾、喜凉、怕高温[33]的生长特性。另外,海拔较高的区域有利于保护银杉野生种质资源,这也与核心适生区分布大多在海拔较高区域的结果相一致。笔者认为,除了气候变化的主要影响,低海拔也不利于银杉形成竞争优势,而较高海拔的核心适生区与较低海拔的其他适生区间断出现,也可能是当前银杉分布区破碎的主要原因之一。

3.3 不同气候情景对银杉未来分布的影响

由本研究结果可知,未来不同气候情景下,银杉的核心适生区面积均呈增加趋势。根据相关学者的研究结果,到21世纪末,全球气候将持续变暖[34],长期的气象观测数据表明未来气候变化下中国区域平均年降水量将有所增加(0～20%)[35]。而最干月份降水量、最冷月最低温都是影响银杉分布的主要因子,也是促进银杉未来核心适生区面积扩张的原因。同时,未来不同气候情景下,银杉的核心适生区面积扩张程度不同。在SSP2-4.5情景下,核心适生区扩张面积较多,在SSP1-2.6、SSP5-8.5情景下,核心适生区扩张面积较少。这说明不同气候情景模式下,银杉动态变化有一定的差异性。根据气候变化影响评估情景工作组的最新成果可知,到21世纪末,在可持续发展的情景(SSP1-2.6)下,全球平均气温将上升2.0 ℃;在中等发展情景(SSP2-4.5)下,全球平均气温将上升3.0 ℃;在常规发展情景(SSP5-8.5)下,全球平均升温4.6 ℃[36]。由本研究预测结果可知,SSP2-4.5情景下银杉核心适生区扩张面积最多,说明未来全球平均气温上升3.0 ℃将更有利于银杉的生长与繁殖。

3.4 银杉种质资源保护与发展建议

我国在应对气候变化方面采取了多项措施,针对林业发展方面,国家林业和草原局制定了林业发展和林业应对气候变化等相关纲要[37]。根据预测结果,银杉在未来具有很大的发展潜力,建议相关部门加强对银杉种质资源的保护与培育,以应对未来气候变化。

银杉在我国黔、湘、桂、渝4个省份具有高度的适生性,首先建议在银杉高度集中的地方实行就地保护,建立并完善银杉自然保护区网络体系,目前已有的保护区包括花坪银杉自然生态保护区、沙角洞银杉自然保护区等。同时通过银杉种子的收集与选育、对银杉苗木实行多途径培育,并在银杉现有分布区附近选择适宜生境进行迁地保护,以缓解银杉濒危现状。目前,云南、浙江、福建等省鲜有银杉分布的相关报道,但其部分区域属于银杉的核心适生区,由此可推断,这些区域在过去曾有过银杉的出现,后因第四纪冰川而消失。建议在这些区域对银杉进行引种栽培,以扩大银杉在全国的种植区域。其次,在后期人工培育银杉的过程中,尽可能在前期对银杉幼苗的生长采取一定的人为干扰,为银杉的培育创造条件。此外,银杉的生长可能还与人类活动因子、地理阻隔、扩散能力以及遗传适应环境变化的速度[1]等因素有关,今后的研究中可从不同角度进行探索,以期完善对银杉这一濒危植物的系统认识和有效保护。