基于优化后的MaxEnt模型对珍贵藏药桃儿七潜在适生区分析

赖文峰, 叶兴状, 文国卫, 施晨阳, 张伟皓, 叶利奇, 张国防

(福建农林大学林学院，福建 福州 350002)

联合国政府间气候变化专门委员会(Intergovernmental Panel on Climate Change, IPCC)第5次报告指出，21世纪末全球平均地表温度将比1986—2005年升高0.3～4.8 ℃[1].研究表明，持续性的气候变暖驱使着植物生理生化特性、遗传结构和群落结构的改变[2-4]，甚至可能会导致植物适生范围收缩，生境破碎化程度增加，物种灭绝风险加剧[5].已有研究指出[6]，气候变化正是引起某些植物灭绝的原因.本世纪末，我国多个珍稀植物种群将受到气候变暖的威胁，极大影响了我国生物多样性的安全性[7].因此，为了应对气候变暖危机，制植物物种多样性保护策略至关重要.

气候对植被种群的影响是当前植物多样性保护的研究热点之一.近年来，多数研究者采用物种分布模型(species distribution models, SDMs)来模拟推测物种潜在适生区[8-10]，研究植被种群对气候变化的动态响应.物种分布模型利用样本分布和环境因子数据，基于特定的算法预测物种生态位，最终以存在概率的形式呈现物种出现的可能性[11].现已开发多种物种分布模型，其中最大熵模型(MaxEnt)具有需求样本量少、预测精度高的优点，适用于研究世界级濒危珍贵物种，如世界自然保护联盟(IUCN)红色名录植物预测[12].目前MaxEnt模型已被广泛应用于梓叶槭(Acercatalpifolium)[13]、四合木(Tetraenamongolica)[14]、闽楠(Phoebebournei)[15]、云南红豆杉(Taxusyunnanensis)[16]等珍贵树种.另外，具有完整物种分布和气候因子数据的前提下，MaxEnt模型可应用于史前重要时期末次间冰期(last interglacial, LIG)、末次盛冰期(last glacial maximum, LGM)、全新世中期(mid holocene, MH)孑遗植物分布格局变化、冰期避难所等研究[17]，为孑遗植物谱系地理学提供了一定的辅助研究信息.

桃儿七(Sinopodophyllumhexandrum)为小檗科(Berberidaceae)桃儿七属(Sinopodophyllum)，多年生草本植物，含有鬼臼毒素成分，是传统的珍贵藏药，药用价值较高，也是合成抗癌药物的原料之一[18]，同时，桃儿七起源古老，也是研究鬼臼亚科系统分类学的重要科研材料之一[19].但因其天然繁殖力弱且人为过度采伐而导致桃儿七居群数量迅速减少[20]，现已收录在《国家重点保护野生植物名录》(第二批)，被列为国家Ⅱ级保护植物[21].已有相关研究主要集中于物种所含化学成分分析[22]、生物学特性[23,24]、遗传多样性[25]等方面.目前关于桃儿七MaxEnt模型研究中，主要集中于当代至未来桃儿适生区变化情况[26-28]，未涉及史前地质时期桃儿七分布和运用R软件的ENMeval语言包进行结果优化.

因此，本试验基于R软件ENMeval语言包优化后的MaxEnt模型预测第四纪3个古时期(LIG、LGM、MH)、当代以及未来不同气候情景(RCP2.6、RCP8.5)与不同时期(2050年、2070年)桃儿七地理分布格局，使用贡献率(percent contribution)及置换重要值(permutation importance)、刀切法(jackknife)检验综合评估制约现代地理分布的环境因子，根据响应曲线(response curve)确定桃儿七的适生范围，探讨关键因子如何制约桃儿七分布.结合ArcGIS 10.4制作不同气候情景桃儿七分布图及其未来气候情景下桃儿七适生区变化图，揭示不同气候情景下，该物种适生区的动态变化，为制定野生珍贵藏药桃儿七种质资源保护策略和亲缘地理学研究提供依据.

1 材料与方法

1.1 数据收集

桃儿七分布数据来源于实地踏查、地方植物志和数据库共享平台.数据库共享平台包含Global Biodiversity Information Facility(GBIF, www.gbif.org)、国家标本平台(NSII, http://www.nsii.org.cn/)、中国数字植物标本馆(CVH, http://www.cvh.ac.cn/)等标本收集平台；中国知网(China National Knowledge Infrastructure, CNKI)、web of science(SCI)等各类文献数据库收录的文献记录.一共获取351条桃儿七样本分布记录.将获得的桃儿七分布数据去除人为栽植分布点、重复分布点、模糊记录分布点等原则，然后将初步筛选的分布点进行校正，最终确定209个桃儿七自然分布点并导出文件备用.

1.2 环境变量因子获取与筛选

中国省级行政区划图来源于国家基础地理信息系统数据库(http://www.ngcc.cn/ngcc/).8个气候情景(LIG、LGM、MH、Current、RCP2.6-2050s、RCP2.6-2070s、RCP8.5-2050s、RCP8.5-2070s)中19个气候因子以及当代地形因子(海拔)均来自WorldClim version 1.4数据库(http://worldclim.org)，其中LIG为130 ka Bp的末次间冰期，LGM为22 ka Bp的末次盛冰期，MH为10 ka Bp的全新世中期，Current为当代时期，RCP2.6和RCP8.5分别代表典型浓度目标(representation concentration pathways, RCP)的低(2100年辐射强迫稳定在2.6 W·m-2)与高(2100年辐射强迫稳定在 8.5 W·m-2)的温室气体排放情景.从国家青藏高原科学数据中心(http://westdc.westgis.ac.cn)选择16个土壤因子[基于世界土壤数据库(http://www.fao.org/faostat/en/#data.)的中国土壤数据集(V1.1)].由于数据库中暂无史前地质时期与未来的地形、土壤数据，对于史前时期仅运用气候因子，而对于未来时期依据参考文献[29]中的假设(未来气候情景下地形、土壤数据与当代保持一致).环境因子空间分辨率精度除LGM时期为2.5′外，其余时期所有环境因子精度均为30″.为避免较多环境因子引起模型过度拟合化导致预测结果不精确，利用IBM SPSS Statistics 25软件进行Pearson相关系数分析，保留相关系数|r|<0.8的环境因子，在|r|>0.8的环境因子中仅保留1个最有可能与桃儿七生长关系最为密切的因子，最终筛选出用于预测桃儿七潜在适生区的环境因子(表1).

表1 桃儿七主要环境因子参数

1.3 模型优化与精度评定

近年来MaxEnt模型被广泛应用，Morales et al[30]对于已发表MaxEnt文章进行了综合评价，结果显示部分文章缺少最佳参数选择，导致结果不精确.Warren et al[31]研究表明，过于复杂或简单的模型会降低推断物种栖息地质量的能力.Muscarella et al[32]曾开发ENMeval语言包用于最佳参数选择，并以8种海葡萄属(Coccoloba)植物为例验证该数据包的有效性.R平台ENMeval语言包针对MaxEnt软件中调控倍频(regularization multiplier, RM)和特征组合(feature combination multiplier, FC)参数进行优化[32,33].本研究利用数据包设定RM从0.5～6，每次增加0.5，共12个调控倍频，同时设定10个特征组合，分别为L、QT、H、HP、PT、QH、LQH、LPT、QHP、LQHPT，其中L为线性(linear)，Q为二次型(quadratic)，H为片段化(hinge)，P为乘积型(product)，T为阈值型(threshold).ENMeval语言包将上述120种参数组合进行检验，模拟结果准确性采用受试者工作特征曲线(receiver operating characteristic curve, ROC曲线)与横坐标所围成的区域即受试者工作特征曲线下面积(area under curve,AUC)进行评价.AUC范围为0～1，值越接近1，表示预测效果越好.一般认为，AUC<0.7表示预测效果极差，0.7～0.8表示预测效果比较精确，0.8～0.9表示预测效果很精确，0.9～1.0表示预测效果极其精确.AUCDIFF(the difference between training and testing AUC)代表训练集AUC与测试集AUC之差，or10(10% training omission rate)表示10%训练遗漏率，用于测试模型对物种分布的拟合度，Akaike信息量准则(akaike information criterion correction,AICc)用于评估不同参数组合的拟合度与复杂度，AICc值最小即delta.AICc=0表示该参数组合模拟结果最优，可用于MaxEnt软件建模.将桃儿七样本分布数据和同一时期的环境变量因子输入MaxEnt中，设定75%随机样本数据用于训练集(training data)，25%用于测试集(test data)，设定10次重复，其余均设定默认值.

1.4 适生度等级分区及面积统计

将每个时期MaxEnt模拟结果中的平均值导入ArcGIS 10.4软件中，对导入图层进行重分类.参考已有研究得出l[26]，0～0.2为非适生区、0.2～0.4为低适生区、0.40～0.65为中适生区、0.65～1.00为高适生区.利用ArcGIS软件进行可视化处理并利用面积制表功能，统计各部分适生区面积.

1.5 适生区空间格局变化

由于史前地质时期(LIG、LGM、MH)仅采用气候数据，为避免研究中各时期数据类型及精度不一致，故未绘制史前到当代时期变化图.桃儿七适生区空间格局变化图均是未来各情景与当代相对比.

物种对于气候与环境的变化产生3种生存情景：消失、维持现状、迁移.根据3种可能发生情景、分布区域定义；物种迁入地区为增加适生区；物种维持现状不发生迁移为保留适生区；物种迁出地区代表丧失适生区[34].将ArcGIS中已分类的不同时期适生度等级分区图中物种分布概率P<0.4赋值为0，P>0.4赋值为1，再依据上述原则继续进行重采样赋值和可视化处理，绘制未来时期桃儿七物种中、高度适生区分布格局变化图并计算其面积.

2 结果与分析

2.1 模型优化结果与准确性评价

MaxEnt软件初始默认调控倍率RM=1，特征组合FC为LQHPT，该模型下delta.AICc= 601.3746.经ENMeval 语言包优化后得到RM=3.5，FC为LQH时，AICc最小，delta.AICc=0，则表明其为最优参数组合(表2)，该组合中avg.diff.AUC和avg.test.or10pct值均小于默认模型，说明该参数组合能降低模拟结果的复杂程度和过拟合情况，提高预测的精确性.本研究采用MaxEnt模型预测当代桃儿七潜在适生区，设定对应优化参数，经10次重复训练得到的训练AUC最大值为0.957 8，平均值为0.953 7，最小值为0.947 4，测试AUC最大值为0.961 2，平均值为0.949 3，最小值为0.937 1.所有AUC>0.9，表明模拟结果准确.

表2 MaxEnt模型ENMeval优化评价指标

2.2 影响桃儿七空间分布的环境变量因子及适宜范围

挑选16个环境因子用于建立模型，包含10个气候因子，1个地形因子，5个土壤因子.根据MaxEnt模型得出各个环境因子贡献率、置换重要值(表1).由贡献率可知：海拔(elev, 51.0%)、年降水量(bio12, 26.5%)、最冷季平均气温(bio11, 8.0%)排名前三，贡献率合计达85.5%.置换值中，排前3位的是海拔(elev, 41.4%)、最湿季降水量(bio16, 17.1%)、最干季降水量(bio17, 15.5%)，置换重要值累计74%.

刀切法检验(图1)表明：仅使用单一环境变量时，测试增益(test gain)与AUC数值最高的3个因子为海拔(elev)、最热月最高温(bio5)、等温性(bio3)，在模型的正则化训练收益(regularized training gain)中排名前三的因子是海拔(elev)、等温性(bio3)、最热月最高温(bio5)；依次排除单一环境变量时，模型的正则化训练收益降低较多的是海拔(elev)、等温性(bio3)、最热月最高温(bio5)、等温性(bio3)，测试增益和AUC数值降低较多的是海拔(elev)、最热月最高温(bio5)、等温性(bio3)，说明这些环境变量因子包含更多与桃儿七生长的相关信息.综合分析可知，海拔、年降水、最冷季平均气温是制约当代桃儿七分布的关键因子，最湿季降雨量、最干季降雨量、最热月最高温、等温性也对桃儿七的分布具有一定影响.

bio1:年均温；bio2:月平均昼夜温差；bio3:等温性；bio5:最热月最高温；bio7:温度年较差；bio11:最冷季平均气温；bio12:年降水；bio15:降水量季节性变化；bio16:最湿季降水量；bio17:最干季降水量；elev:海拔；t_cec_soil:土壤的阳离子交换能力；t_oc:表层土壤有机碳含量；t_sand:表层含沙量；t_ph_h2o:表层土壤酸碱度；t_usda_tex_class:表层质地分类

根据响应曲线(图2)确定各影响桃儿七分布关键因子的阈值(一般认为逻辑值≥0.5对应的环境因子范围最适于物种生存[35])：海拔2 000～4 200 m，年降水量500～1 000 mm，最冷季平均气温-8～5 ℃.

图2 关键环境因子响应曲线

2.3 不同时期桃儿七潜在适生区及其面积变化

8个不同气候下，桃儿七整体分布范围大体一致，主要分布于西藏、云南、四川、青海、甘肃、宁夏、陕西等地(图3).在过去3个气候下，从末次间冰期到全新世中期，桃儿七总适生区面积逐渐增大，但高适生区呈收缩趋势，具体表现为总适生区面积由112.32×104km2增加到121.75×104km2，而高适生区面积从20.87×104km2减少到14.85×104km2(表3).当代潜在总适生区面积114.99×104km2，其中低适生区面积48.15×104km2，中适生区面积49.79×104km2，高适生区面积17.05×104km2，分别占全国地表面积的5.02%、5.19%、1.78%.与当代相比，未来RCP2.6情景下，桃儿七总体适生区呈不显著收缩趋势，2050年与2070年总面积分别减少0.14%和0.22%；RCP8.5情景中2050年与2070年总适生区面积均增加，但仅表现为低适生区显著扩张，而中、高适生区范围呈现收缩趋势，总体适生区分别增加0.08%和0.38%.

A:末次间冰期(LIG);B:末次盛冰期(LGM);C:全新世中期(MH);D:当代;E:RCP2.6-2050年;F:RCP2.6-2070年;G:RCP8.5-2050年;H:RCP8.5-2070年.

表3 不同气候情景下桃儿七适生区面积统计1)

2.4 不同气候变化下桃儿七适生区变换格局

由图4和表4可知，未来情景下桃儿七丧失区域主要位于中、高适生区的边缘部分，大致位于青海东部与南部、四川西北部、西藏东部及云南中部与北部，其中RCP2.6情景的2050年与2070年丧失面积均大于RCP8.5同时期，分别为4.78×104km2和4.72×104km2，丧失率为7.15%和7.06%.未来桃儿七增加地区主要位于甘肃中部与东南部、青海东部、四川中部、云南西北部地区，增加最多的两个时期分别为RCP8.5-2070年和RCP2.6-2070年，增加面积为3.02×104km2和2.92×104km2，增加率为4.53%和4.37%.而未来保留区域整体不改变，主要分布于甘肃南部、青海东部、四川中部与西部、西藏东部、云南北部地区，其中RCP8.5情境下的两个时期保留面积最大，分别为63.67×104km2和62.78×104km2，保留率95.25%和93.92%.随着时间推移，未来两种浓度背景下桃儿七中、高适生区面积呈现先减少后增加趋势，但两个等级适生区总面积均小于当代面积，呈现下降趋势.

A:RCP2.6-2050年;B:RCP2.6-2070年;C: RCP8.5-2050年;D:RCP8.5-2070年.

3 讨论

3.1 模型可靠性

MaxEnt模型优化分析表明，濒危植物桃儿七潜在适生区范围主要在陕西南部、宁夏南部、甘肃中南部、青海东部和东南部、四川中西部、云南北部、西藏青藏高原东南部等地区，这与郭彦龙等[26]、张海龙[27]、吕蓉等[28]使用MaxEnt模型预测范围基本一致，且本试验的AUC>0.9，说明本次结果可信度高，对于以后桃儿七种质资源调查与收集具有指导意义.

3.2 关键因子

关于制约桃儿七分布的关键因子研究中，郭彦龙等[26]研究表明，最热季平均气温、年降水量、气温季节性变动系数是影响桃儿七分布的主导因子；张海龙[27]研究表明，海拔、最热季平均温度、极端低温对桃儿七分布起到重要作用；吕蓉等[28]研究表明，等温性、4月份降雨量、海拔为制约桃儿七分布的主要环境因子.而本试验结果表明，影响桃儿七分布的重要环境因子是海拔、年降水、最冷季平均气温.造成各研究结果不同的原因可能是：①采用环境因子类型不同.除气候因素外，地形、土壤等均会影响物种分布[36].因此，采用不同环境数据类型进行预测可能得到不同的结果；②参与建模的环境因子是否进行相关性分析.前人研究过程中未运用相关性分析进行环境因子筛选，变量之间的相关性可能会影响分析结果[26-28]；③保留并参与建模的物种居群样本范围和数量也会导致模拟结果不同[37].

海拔和温度的共同作用对桃儿七生长有着至关重要的影响[38-40]，本试验结果表明，海拔为制约桃儿七分布的最关键因子，其最佳适宜范围为2 000～4 200 m.海拔对桃儿七的影响主要是通过温度而实现[41].高海拔地区低温引起的雪害影响着桃儿七的繁殖策略，其繁殖策略中存在种子数及结实率与海拔呈正相关，单粒种子干重和果实内种子数与海拔呈负相关的规律[42].同时温度控制着桃儿七的生长节律，随着海拔上升，年平均气温和极端低温越低，其物候期随着海拔每升高100 m时间将推后3～5 d[43].年降雨量是降水因子中对桃儿七生长的最重要因子.本试验结果表明，降水因子累计贡献率为30.8%，桃儿七最适于生在长年降水量400～900 mm、降水集中于6—9月的地区[20]，本研究预测年降雨量最适范围(500～1 000 mm)与其大致相符.当代桃儿七适生区主要分布在西南地区，气候大多为雨热同期.研究表明，雨热同期的季风气候影响高寒地区植被的返青期与雨季始期高度重合[44].雨季中桃儿七迅速生长发育，叶片展开，形成大量根系.同时雨季中也具有桃儿七最适生长温度，能够积累越冬和繁殖所需要的营养物质[42].另有研究表明[43]，桃儿七能根据生长环境中水分状况调整繁殖方式，在湿润的地区可以进行营养繁殖，而在干旱的地区则不进行营养繁殖.最冷季平均气温(-8～5 ℃)成为影响桃儿七分布的第3重要因子.极端气候(最冷季平均气温)也是影响植物分布的重要因素之一，由于极端气候的不连续性，常导致物种分布范围发生不同程度的改变[45,46].综合分析，海拔、年降雨量、最冷季平均气温是制约桃儿七分布的关键因子，这与桃儿七适生于高海拔、湿润冷凉的生物学特性相一致.

3.3 气候变化对桃儿七适生区影响

末次盛冰期气候寒冷干燥，该时期极端气温导致部分生物无法适应气候而灭绝，而另一部分生物选择迁移至避难所.生物避难所是具有一定稳定性的生境地理区域，能为冰期生物提供理想的生存繁衍环境，具有丰富遗传多样性的特征[47].横断山脉地区地理位置特殊，位于青藏高原东南边缘，冰川事件使其生境多样化，植物区系丰富，使其有可能成为避难所[48]，例如桃儿七碧塔海自然保护区地处横断山腹地，受冰期影响较小，成为众多古老植物的理想避难所[49].在末次盛冰期，桃儿七高适生区分布相较于其它史前地质时期更为集中，主要位于在东喜玛拉雅地区与横断山区，且多时期下桃儿七位于此处的地理范围大致不变.因此，本文推测东喜马拉雅地区和横断山区存在桃儿七的气候避难所.李永[50]利用扩增片段长度多态性(amplified fragment length polymorphisms,ALPFs)和cpDNA(chloroplast DNA)进行桃儿七亲缘地理研究，结果表明，桃儿七在东喜玛拉雅地区与横断山两个地区均存在避难所，本研究结果与之保持一致.全新世中期较末次盛冰期更为温暖湿润，孑遗生物度过冰期后将从避难所中迁移并重新分布[51].良好的水热条件使得桃儿七更好满足生长发育对热量和水分的要求.气候回暖驱使着桃儿七适生区向外扩张，整体适生区范围更大.未来4个时期桃儿七适生区范围变化不显著，这与桃儿七有较宽的生态幅有关，可以适应多种外界环境[43].未来桃儿七整体适生环境范围基本不变的情况下中、高适生区略有缩减，且适生范围有南部减少而北部与西部增加的趋势，预测未来桃儿七潜在适生区可能会向西部北部等高海拔高纬度地区迁移.郭彦龙等[26]研究表明，桃儿七适生区北部为河西走廊等荒漠区，干旱恶劣的环境条件阻碍了桃儿七向高纬度迁移，故未来桃儿七适生范围极大可能会向西延伸至高海拔地区.

3.4 桃儿七保护策略

通过对比适生区面积变化可知，本世纪末两种气候浓度(RCP2.6、RCP8.5)下桃儿七中、高适生区面积小于当代时期，表明气候变暖会对未来桃儿七居群构成威胁.因此，制定桃儿七保护策略至关重要.未来桃儿七物种可能在西藏东部、四川中部、青海西部、甘肃南部等地增加适生地区，建议将这些地区规划为引种示范区，充分发挥其生态和经济效益.其次，因桃儿七生态幅较广，故可尝试在一些低适生区进行引种试验，扩大种植范围.亟待收集云南、四川和甘肃等桃儿七丧失区的种质资源，以期尽可能保护桃儿七遗传多样性.对甘肃、四川及云南等地的桃儿七高度适生区制定相应保护政策，加大保护宣传力度，尽可能减少人为采摘和放牧活动，保证桃儿七数量不减少和生境不被破坏.除了传统培育保护措施之外，尽早实现人工合成鬼臼毒素代替植物提取以解决桃儿七供需不平衡的现状，对于保护其种质资源也具有重要意义.

4 结论

基于MaxEnt软件模拟分析预测影响当代桃儿七适生区范围的环境变量因子主要为海拔、年降水量、最冷季度平均气温.当代桃儿七潜在适生区主要分布于我国西藏、云南、四川、青海、甘肃、宁夏、陕西等地.末次盛冰期，桃儿七为应对极端气候迁移至横断山脉地区和东喜马拉雅山脉地区进行避难.度过冰期后，桃儿七适生区范围扩大.未来桃儿七适生区将可能向高海拔地区迁移，且会因气候变暖导致适生程度较高的地区丧失.