基于MaxEnt优化模型的闽楠潜在适宜分布预测

叶兴状,张明珠,赖文峰,杨淼淼,范辉华,张国防,陈世品,刘 宝,*

1 福建农林大学林学院,福州 350002

2 福建师范大学地理科学学院,福州 350007

3 福建省林业科学研究院特色资源研究所,福州 350012

气候变暖影响全球水热分布和土壤微气候,而植物分布格局主要受水份、光照、土壤、海拔、坡向和CO2浓度等因子的制约[1]。近年来,气候变暖对植物分布的影响是植物学、生态学及其交叉学科的研究热点之一。第四次工业革命以来,全球气候变暖趋势加剧,导致其中部分物种的适生区面积大大减小和生境的破碎化,致使许多物种濒危甚至灭绝,严重影响生态安全和生物多样性[2- 3]。据统计,目前世界上约六分之一的物种受不同程度的威胁,有研究表明,随着未来全球气温的升高,濒危物种的灭绝风险将大大增加[4]。因此,在全球气候变暖的背景下,探索濒危物种的分布格局如何受未来气候变化的影响,将有助于制定濒危物种的保护策略及种质资源的保存计划,以期规避灭绝风险,维持人类赖以生存的生态环境安全。

物种分布模型是研究物种适生区域受气候环境变化影响的重要工具[5]。目前,国内外学者采用不同的算法开发了如DOMAIN(领域模型)、GARP(遗传算法模型)和MaxEnt(最大熵值模型)等物种分布模型,其中,MaxEnt模型用于处理群体性温度因子及降水因子的失真度小、稳定性好,是一种理想的预测模型,也是目前应用最普遍的物种分布模型之一[6],已被广泛应用于濒危植物适生区变化预测,如红豆树[7](Ormosiahosiei)、半枫荷[8](Semiliquidambarcathayensis)和四合木[9](Tetraenamongolica)等。因此,采用MaxEnt模型研究气候变化对濒危植物适生区的影响,进而制定相应的保护策略具有重要参考价值。

闽楠(Phoebebournei)属樟科(Lauraceae)楠属(Phoebe)植物,是我国特有二级濒危珍稀树种,其树干高大通直,材质致密,是不可多得的珍贵用材树种,同时还是实现碳中和的理想树种,具有很高的经济和生态价值[10-11]。天然闽楠林因历年遭受严重的滥砍滥伐,导致其野外生境破碎化现象异常严重,现有天然林资源几近枯竭[10- 12]。与历史记载相比,其分布范围严重收缩,除福建沙县罗卜岩自然保护区和湖南平江等少数地点的天然林外,野外已难见集中成片的闽楠天然林[13- 14],因此,加大对天然闽楠资源的保护力度已迫在眉睫。此外,我国南方林区多为杉木或马尾松人工纯林,林分结构单一、病虫害日趋严重、森林生态系统功能稳定性差,我国战略储备林规划指出,应对战略储备树种采用培育复层异龄混交林等措施,以期提高我国森林资源质量,而闽楠是理想的阔叶造林树种之一[15]。然而,目前闽楠的研究主要集中在育苗造林、群落生态、良种选育和遗传多样性等方面[15- 19],目前还尚未有关于闽楠适生区预测的相关报道,而适生区模拟预测将有利于揭示闽楠未来分布变化规律及影响其分布的关键环境因子,对闽楠核心种质资源库完善、资源保育和引种栽培具有重要意义。

本研究在野外实地调查和前人对闽楠资源分布调查的基础上,应用基于Enmeval数据包辅助选择MaxEnt最优参数设置来模拟推测闽楠的地理分布,同时分析影响闽楠分布的关键环境变量,并进一步分析其适生区的环境特征,以期为闽楠的资源保护、天然种群恢复和引种栽培提供科学依据。

1 材料与方法

1.1 研究区域

闽楠是我国特有珍贵树种,其天然分布北界自湖北沿神农架、大别山直至江西北部山丘,南界为南岭山脉的广东、广西丘陵至花坪国家自然保护区,东界为浙江雁荡山、乌岩岭至福建东北部丘陵,西界自重庆西部山区、贵州中西部山区至广西西北部山区[17]。闽楠天然分布区的地理坐标范围为23°—30°N,105°—122°E[16],考虑到要为科学引种提供参考,因此以100°—130°E,10°—35°N为预测区域。

1.2 样本数据的搜集与筛选

课题组在1998—2019年陆续对福建、浙江、江西、湖南、广西和广东等6省部分闽楠天然种群进行了实地调查,搜集分布记录133条；查阅已出版的文献资料,搜集闽楠分布记录46条[11- 19]；检索全球物种多样性信息库(http://www.gbif.org)、国家标本平台(http://www.nsii.org.cn)、教学标本资源共享平台(http://mnh.scu.edu.cn/main.aspx)、中国数字植物标本馆(http://www.cvh.org.cn)分别搜集闽楠记录170条、635条、113条和653条,合计1750条。参照李垚[20]的方法对上述全部记录进行筛选,删掉非天然分布记录及重复点,并删除2000年之前的记录点,为降低群集效应带来的误差,每个网格(10 km×10 km)只保留一个分布点,最终得到有效样本186个。

1.3 环境变量筛选和数据处理

气候数据选择来自于WorldClim数据库(http://worldclim.org)的第六次国际耦合模式比较计划(CMIP6)模式的共享经济路径。为提高预测准确性,选取现代及未来两个时期(2050s和2090s)各3个气候情景(SSP126、SSP245和SSP585)进行等模拟,空间分辨率为2.5 arc-minutes。本研究气候数据采用第二代国家气候中心中等分辨率气候系统模式(BCC-CSM2-MR)生成,该模式对中国降水与气温的模拟较为可靠[21]。采用Spearman Correlation Analysis和VIF方差膨胀因子,选择相关系数小于0.8的因子。生态极限和植物的实际分布密切相关,野生闽楠多见于山地沟谷阔叶林中,海拔在40—1 220 m,且多分布于海拔600 m以下丘陵及低山,结合植物生存环境需求及其生理特性,在相关系数大于0.8的变量中保留最重要的生态因子[17]。最终,保留9个变量用于建模,包括年均温(bio1)、昼夜温差日均值(bio2)、等温性(bio3)、最暖月最高温(bio5)、最冷月最低温(bio6)、最湿季度平均温度(bio8)、最干月降水量(bio14)、降水量变异系数(bio15)和最暖季度降水量(bio18)。

1.4 模型建立、优化和评估

为使闽楠出现的概率接近正态分布,选择75%的数据用于模型训练,剩下的数据用于模型测试,其他值为默认值。使用Rv3.61中的Emneval数据包优化MaxEnt模型[22],将调控倍频设置为0.5—5,每次间隔0.5,一共10种调控倍频；MaxEnt模型提供了5种特征,分别是线性特征、二次型特征、片段化特征、 乘积型特征和阈值性特征。采用8个特征组合。Enmeval 数据包将上述 80种参数组合进行测试,最终采用最小信息准则AICc(The minimum information criterion AICc value, delta.AICc)模型,采用训练AUC和测试AUC值之差(Difference between the AUC values, avg.AUCDIFF)检验模型的拟合度与复杂度[23]。

优化完成后,采用优化参数进行模拟和预测闽楠不同时期适生区。利用受试者曲线下的面积评估MaxEnt预测的准确性。AUC值的取值范围为0—1,值越大意味着模拟可信度越高,AUC=0.5表示该模型判断完全失败,小于0.8表示可信度低,0.9—1.0表示极准确[7- 9]。采用平均逻辑值将适宜度分为(Ⅰ：0—20%；Ⅱ：20%—40%；Ⅲ：40%—60%；Ⅳ：60%—100%)4级,分别表示非适生区、一般适生区、中度适生区和高度适生区。

1.5 物种适生区空间格局变化

本研究将逻辑值Q≥40%作为潜在适生区,即适生区包括了中度适生区和高度适生区,用“1”表示；将逻辑值Q<40%作为非适生区,用“0”表示。本研究定义了4种物种适生区变化的类型：新增适生区、丧失适生区、保留适生区和非适生区。未来面积变化均基于当前闽楠适生区面积来比较计算。矩阵值0→1为新增适生区,1→0为丧失适生区,1→1为保留适生区,0→0为非适生区。具体方法参考张殷波[24]的方法,最后将矩阵变化值加载入Arc GIS10.4中,实现闽楠适生区空间格局变化的可视化表达。

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

参考李垚[19]的方法,裁剪闽楠现代潜在适生区的环境变量作为参考图层,采用多元环境相似度面和最不相似变量分析探讨未来6个温室气体排放情境下的气候异常区域及引起潜在适生区改变的关键气候因子。多元相似度能够反映特定时期某点的气候状况与参考图层气候状况的相似程度,负值表明该点至少有一个环境变量的数值超出了参考图层对应值的范围,称之为气候异常点,最大值100代表该点的气候完全正常。某点的最不相似变量指的是相似度最低(即异常程度最高)的变量,异常程度最高的环境变量很可能是引起地理分布变迁的关键因素[24]。

2 结果与分析

2.1 模型优化及准确性评价

本研究基于闽楠186个分布点和9个生态因子,运用Enmeval 数据包调用MaxEnt对闽楠现代及未来的潜在分布区进行预测。当FC(特征组合)=PT(乘积型特征),RM(调控倍频)=0.5时,delta.AICc=0,该参数下模为最优模型(表1)。优化后的最小信息准则AICc(the minimum information criterion AICc value, delta.AICc)和训练集AUC与测试集AUC差值(Difference between the AUC values, avg.diff.AUC)均低于默认设置,表明优化后的参数有利于减低模型的拟合度和复杂度,故本研究选取FC=PT,RM=0.5作为建模的参数设置。在该参数下,MaxEnt预测结果显示受试者工作曲线AUC值为0.9846±0.0037,表示预测结果极准确。

表1 MaxEnt模型Enmeval优化评价指标Table 1 Evaluation metrics of MaxEnt model generated by Enmeval

2.2 环境变量的重要性

从表2可知,贡献率最重要的6个环境因子累计占比达86%,依次分别为最干月降水量(bio14,37.4%),适宜范围为大于17.93 mm；降水量变异系数(bio15,15.8%),适宜范围为大于46.75—68.35；最冷月最低温(bio6,12%),适宜范围为-3.32—6.31℃；最暖季降水量(bio18,11.1%),适宜范围为大于72.32 mm；湿季度平均温度(bio8,10.2%),适宜范围为40.69—81.86；昼夜温差日均值(bio2,9.6%),适宜范围为小于8.89 ℃。置换重要值位居前6位的最冷月最低温(bio6)、最干月降水量(bio14)、最湿季度平均温度(bio8)、海拔(elev)、等温性(bio3)和昼夜温差日均值(bio2)。

表2 闽楠主要环境因子参数Table 2 Environmental variables and their main parameters of P. bournei

2.3 现代和未来潜在适生区

186条分布记录的平均逻辑值为0.49,最高值为0.76(湖南蓝山),最低值为0.11(广西田林),分别属于高度适宜区、中度适宜区、一般适宜区和不适宜区的分布点比例分别为21.76%、57.41%、16.66%和4.17%。从表3可知,当前潜在适生区的面积为54.32×104km2,其中高度适生区面积和中度适生区面积分别为3.37×104km2和21.60×104km2。如图1,当前气候条件下闽楠高度适生区主要位于广西东北部、湖南南部、广东北部和福建中北部局部。中度适生区主要分布于福建中东部、江西中南部、贵州东南部、湖南西部及南部、广东北部和广西东北部。

表3 不同时期闽楠适生区面积变化/(×104 km2)Table 3 Changes of suitable area of P. bournei in different periods

图1 闽楠现代适生区Fig.1 Modern suitable habitats area of P. bournei

2.4 未来潜在适生区空间变换格局

将6个未来不同气候情景下闽楠适生区的空间格局与现代闽楠适生区进行对比分析(表4、图2和图3),结果表明：新增适生区在2050s-SSP585排放情景下分布面积最大,新增率达到81.32%,新增面积主要集中在湖南中西部、贵州东部；中高度适生区约为42.17×104km2,为现代的1.69倍,保留面积21.78×104km2,保留率达87.21%；潜在适生区丧失面积最少,丧失3.01×104km2,丧失率仅为12.06%,主要丧失区域为广西北部、贵州中部和江西中部。

表4 不同气候情景下闽楠适生区空间变化Table 4 Changes in suitable distribution area of P. bournei different climate scenarios

图2 MaxEnt模型预测的不同气候情景下闽楠潜在适生区Fig.2 MaxEnt model predicted potential suitable growth areas of P. bournei under different climate scenarios

图3 不同气候情景下闽楠适生区空间变换格局 Fig.3 Suitable distribution of P. bournei under different climate change

2090s-SSP126排放情景下,新增面积最小,新增率为29.32%,新增面积主要集中于湖南南部、江西西部以及福建中西部,中高度适生区为22.42×104km2；该时期潜在丧失面积最大,达到10.07×104km2,丧失率为40.32%,主要丧失区域包括贵州中南部、广西东北部、湖南西部、江西南部和浙江南部。

2090s-SSP245和2090s-SSP585情景下,保留适生区面积分别为第二及第三大,保留面积分别为17.22、16.73×104km2,保留率分别为68.96%、66.99%。两个情景的保留区均主要位于闽东北、桂东北和湘南。两个情景的新增适生区均主要位于黔东北及湘西,丧失区均主要位于黔东南、桂东北及浙东南。前者适生区较现代增加3%,而后者较现代减少3.1%。

2050s-SSP126和2050s-SSP245情景下,保留适生区面积几乎相同,均较小,保留面积分别到16.45×104、16.57×104km2,保留率分别为65.89%、66.35%。这两个时期丧失区较大,前者丧失区主要位于闽浙交界、赣中南部及湘西,后者丧失区主要位于黔东、赣中南部及湘西。前者新增适生区主要位于闽西、湘西和黔东,而后者新增适生区为闽西和湘西,前者适生区比现代略有增加,后者适生区比现代增长6%。

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

如图4所示,在2050s-SSP126、2050s-SSP245、2050s-SSP585、2090s-SSP126、2090s-SSP245和2090s-SSP585等6个排放条件下,闽楠186个现代分布点的平均多元相似度分别为2.53、19.17、10.03、15.46、7.46和-2.06,多元相似度为负值的点比例分别为18.98%、5.56%、18.98%、6.48%、26.85%和61.11%,表明2090s-SSP585排放条件下,气候异常程度最高,2050s-SSP126情景次之,2050s-SSP245情景下气候异常程度最低。在2050s-SSP126情景下,现代适生区内的气候异常区域主要出现在中北部,最不相似变量为bio3和bio5。在2090s-SSP126情景下,现代适生区内的气候异常区主要出现在中北部,最不相似变量为bio5、bio13和bio15。在2050s-SSP245情景下,现代适生区内的气候异常区主要出现在中北部,最不相似变量为bio2、bio5和bio15。在2090s-SSP245情景下,现代适生区内的气候异常区主要出现在中北部和南部,最不相似变量为bio5和bio2。在2050s-SSP585情景下,现代适生区内的气候异常区主要出现在中部、南部和北部,最不相似变量为bio5和bio15。在2090s-SSP585情景下,现代适生区内除西部和东部少数区域外,其他均为气候异常区,最不相似变量为bio5(图5)。

图4 MaxEnt模型预测的不同排放情景下闽楠多元相似面(Mess)Fig.4 Multiple similarity surfaces (Mess) of P. bournei under different emission scenarios predicted by MaxEnt model

图5 MaxEnt模型预测的不同排放情景下闽楠最不相似变量(Mod)Fig.5 The most dissimilar variable (Mod) predicted by MaxEnt model in different emission scenarios of P. bournei

3 讨论

3.1 模拟结果的可靠性

生态位模型模拟结果的可靠性主要取决于三个方面：首先是模型的选择,不同模型预测精度差异较大,经验证,MaxEnt模型预测可靠性较其他模型更高[26- 28]；其次是样本采集范围及样本覆盖均匀度,采样点的选择越接近物种真实分布情况,预测精度越高；第三是环境因子的来源、种类及其精度,模拟时选取因子种类和数量越多,预测精度就越高。调用Enmeval数据包对MaxEnt进行优化可降低过拟合程度和抽样偏差、从而提高预测精度[29]。本研究经优化后,AUC =0.9846±0.0037,大于0.9,表明预测精度很高。此外,本研究预测的现代闽楠潜在适生区范围与葛永金得到的结论相近[23],因此模拟结果是可靠的。

3.2 制约闽楠分布的主要环境因子

本研究的8种主要评估途径中,降水因子重要性列入第一位6次、第二位3次,温度因子4次列入第一位、第二位6次,因此,降水和温度共同制约着闽楠现代分布格局。虽然单从重要值来看,分布格局主导因子为最干月降雨量(bio14),但闽楠适宜的温暖指数为141.7—197.1℃·月[30],且本研究表明温度因子的置换重要值(Permutation importance)累计占比达86.7%,说明温度也是制约闽楠分布的重要因子。本研究在进行潜在适生面积预测时,虽然未将城镇建设和工农业用地排除,但这并不影响本次模拟的意义,对发掘潜在闽楠天然分布区和指导闽楠人工林营建均具有一定的参考价值。相关研究证明,MaxEnt模拟预测存在潜在适生区偏大的现象[7- 8],导致这一现象是由于多种因素综合作用的结果,物种的空间分布不仅与降水和温度有关,还受紫外线、生物、人为活动和物种扩散等因素制约,未来的研究中,增加上述变量可能会提高预测的准确性。

贡献率、置换重要值、单独使用该因子的正则化训练增益和单独使用该因子的测试增益均表明,最干月降水量是影响闽楠潜在地理分布的首要因素,这与葛永金的研究结论相一致[31]。吴显坤等[32]研究发现影响浙江楠(Phoebechekiangensis)地理分布的首要因素是最干季降水量,而闽楠和浙江楠同为楠属濒危植物,亲缘关系较近,而且最干月降水量月与最干季降水量两因子相关性很高。因此本研究得出最干月降水量(bio14)是制约闽楠分布格局的关键因子的结论是可靠合理的。此外,本研究发现最冷月最低温、最冷季降水量、年平均气温、降水季节性和昼夜温差月均值也是制约闽楠分布格局的主要环境因子。

气候因子并非相互独立地影响闽楠的生长与分布。多元环境相似性面和最不相似变量分析表明,6个情景闽楠分布区北侧气候异常程度均比其他地区强烈,北侧最不相似变量主要是最暖月最高气温和最湿季平均气温,除2050s-SSP585情景外,与此对应的闽楠未来其他5个情景新增潜在适生区均与丧失区面积无显著性差异,温度因子(主要指最冷月最低温和最热月最高温)是制约闽楠适生区南扩的重要因子。然而,闽楠分布变迁是由多个因子共同作用的结果,这一现象尤其在2050s-SSP245和2090s-SSP126情景下表现尤为突出,温度和降水因子共同制约其地理分布的变迁。

3.3 闽楠潜在适生区的变化格局

在未来6种温室气体排放情景下,MaxEnt预测的闽楠潜在分布区都将北移, 其中在2050s-SSP585情景下,其潜在适生区面积扩大较明显,这可能是由于闽楠为喜暖湿树种,其现代适生区北部毗邻区气温和降水的增加,可能为其提供了更适宜的生长环境。相关研究也表明楠木[33](Phoebezhennan)在2050s的潜在适生区面积比现代更大；桫椤[34](Alsophilaspinulosa)潜在适生区在2090s将逐渐向西北迁移,且4个情景面积也比现代更大。然而,预测植物潜在分布还需结合土壤、地形、植被、人为干扰和种间竞争等因素,综合生物与非生物因子,才能更有效地揭示闽楠的未来潜在适生区。未来各情景下,闽楠生境破碎化有加剧的趋势,特别是2090s-SSP126情景下,闽楠潜在适生区丧失面积占现代适生区面积的40.32%,适生区明显北移,生境破碎化最严重,说明气候变暖对闽楠分布构成较大威胁。本研究表明,6种未来温室气体排放情景下,现代闽楠南端潜在适生区均出现不同程度丧失的现象,且闽楠适生区均往北迁移,黔桂两省(区)闽楠丧失区比较大,而黔桂又是闽楠古树分布的集中区,因此,开展闽楠的抗逆性研究已刻不容缓。此外,本研究还表明,虽然闽楠潜在丧失区与新增适生区无显著差异(除2050s-SSP585情景外),但闽楠未来生境将变的更加破碎化,这一结果与浙江楠、丁香叶忍冬(Loniceraoblata)和藤枣(Eleutharrhenamacrocarpa)等濒危植物的研究结论相一致[32,35- 36],说明全球气候变暖将加剧闽楠资源的丧失,亟待从基因组和蛋白组等分子层面提出有效的闽楠种质资源保护方案,以减缓日益变暖的气候状况对闽楠生境的影响。此外,种源传播在不受阻碍的情况下,闽楠有可能向重庆东南部、湖南、湖北、台湾中部及安徽南部等北部潜在边缘分布区扩展,因此潜在适生区可将闽楠作为改造针叶人工纯林重要树种之一。

3.4 闽楠的资源保护

前人研究发现[31,37],闽楠致濒的可能原因有：首先,种子(幼苗)在高湿度生境中受到土壤病原菌的侵染而易变质,而且容易遭受动物捕食,加之闽楠种子休眠特性,从而导致较低的野外闽楠种子发芽率和幼苗存活率低；其次,闽楠是中等耐阴树种,在高郁闭度的天然林环境下,可能成为制约其幼苗生长的主要因素,导致其幼苗的竞争力低于群落其他耐阴树种幼苗；再次,由于闽楠经济价值高,历史上长期遭受大规模滥砍盗伐。野外调查发现,大多数闽楠古树并未分布在已建立的自然保护区内,而是以风水林的形式破碎化地分布在村镇和寺庙旁,加上闽楠种子主要靠重力散布,导致基因流较低[33-34,37],而本研究模拟结果表明,随着气候变暖,未来闽楠分布破碎化现象越来越严重,闽楠资源保护面临严峻挑战。综合闽楠现代适生区图,未来6个温室气体排放情景下适生区图和空间变换格局,发现闽楠潜在适生区的稳定区为闽中北部、桂东北部、赣中南部、湘赣交界、湘南和黔东等山地沟谷阔叶林,尤其闽中北部、湘南和桂东北山区适生度最高,建议将上述地区列入优先保护区域。刘丹[19]利用EST-SSR和g-SSR两种分子标记,均得出闽楠具有较高的遗传多样性,然而其天然分布破碎化程度较高,说明造成其濒危的原因可能是由于滥砍滥伐、生境破坏和其特殊的自身繁育特性导致的[18, 38- 39]。因此,建议应尽快对闽楠天然种群采取迁地保护和就地保护相结合的立体保护模式,如闽楠盗伐现象应当引起重视[40-41],同时应建立闽楠保护小区,并定期巡查。本研究团队已在福建南平来舟林场和顺昌洋口林场建立闽楠种子园及种质资源库,下一步将加快贵州、广西和湖南等地闽楠种质资源的收集与保存工作,以期为更好地保护和利用闽楠种质资源奠定坚实基础。