当前(1970—2000)与未来气候变化情境下中国醉马芨芨草潜在分布预测

黄睿杰， 张春艳， 温雨婷， 吴晨晨， 路 浩， 赵宝玉

(西北农林科技大学动物医学院， 陕西 杨凌 712100)

全球变暖对温度和降水格局的改变影响着植物的分布格局，因为它们的分布很大程度上取决于非生物因素，如降水、温度、土壤和海拔等[1]。因此，气候变化对植物分布的影响是近年来研究所关注的焦点。对气候变化如何影响植物分布的研究有助于农业生产和生物多样性保护，并促进生态系统的可持续性。

物种分布模型(Species distribution models，SDMs)可以将物种发生点数据与相关环境变量相结合，预测物种适宜生存的环境[2]。最大熵(Max entropy model，MaxEnt)模型是利用最大熵原理预测物种生境的SDMs，其假设物种发生点由物种分布概率确定，而此概率是由环境因子对发生点的约束性计算出来的[3]。过去的研究已证实MaxEnt相比其他SDMs具有更强的预测能力和准确性[3-4]。由于其具有方便易用、对样本数量需求小和预测表现优越等优点，MaxEnt已广泛应用于生物多样性保护、外来入侵物种防治、病虫害防治和作物栽培等方面。但是目前在毒害草防控领域的运用还较少。

醉马芨芨草(Achnatheruminebrians)，俗称醉马草、醉针草、醉针茅等(图1)，属于禾本科芨芨草属多年丛生草本植物，是分布于我国西北、华北天然草地的主要毒害草之一[5]。由于长期的过度放牧和干旱气候影响，醉马芨芨草在退化草地竞争中占据优势，种群数量和生长面积不断增加，降低草原生产力，放牧动物误食或因饥饿被迫采食醉马芨芨草易引起急性中毒，从而严重影响着草地畜牧业发展和草地生态安全[6]。然而醉马芨芨草在药用、造纸、饲草化利用、机场驱鸟、生态等方面又存在其利用价值[6-7]。目前关于醉马芨芨草的研究多为对其化学成分、生物特性、内生真菌、毒理学等方面研究[8]，而关于其不同气候背景下适生区分布格局的研究尚未报道。本研究将基于优化的MaxEnt模型对醉马芨芨草当前和未来(2050s和2070s)4种气候情景下适生区分布进行模拟和预测，并探究影响其分布的关键环境变量。研究将为醉马芨芨草的防治和资源利用、草地生态学研究和草地放牧管理提供理论依据和参考。

图1 醉马芨芨草单株及其生境Fig.1 Individual plant and habitat of A. inebrians

1 材料与方法

1.1 醉马芨芨草发生点数据

研究所用醉马芨芨草发生点数据共515个，来源于：(1)本课题组2014—2021年在中国西北地区的野外调查，使用手持式GPS(GARMIN GPSMAP 621 sc)记录植物群落野外生长的经纬度；(2)中国数字植物标本馆(https://www.cvh.ac.cn/)和全球生物多样性信息机构(Global Biodiversity Information Facility，https://www.gbif.org/)的网络数据库。对于一些没有地理坐标但有准确的地理描述的标本，我们使用百度地图坐标拾取系统(https://api.map.baidu.com/lbsapi/getpoint/index.html)查询坐标。

MaxEnt模型的一个假设是物种的发生点在所有区域都是系统或随机抽样的[3]。但是在实际调查中，往往会受到各种原因导致采样偏差，从而使模型过拟合影响预测能力[9-11]。因此，在删除错误点后，采用ENMTools软件(https://github.com/danlwarren/ENMTools)对发生点数据进行筛选[12-13]，最终每个5 km×5 km的栅格内只含1个发生点数据。

1.2 环境变量数据

研究考虑了60个可能影响醉马芨芨草分布的环境变量，包括19个生物气候变量、36个气候变量、5个土壤变量以及海拔(表1)。生物气候变量虽然由气候变量衍生，但是二者并不相同[14]。当前(1970—2000)与未来(2050s与2070s)的生物气候、气候变量以及海拔数据来源于世界气候网站(Worldclim,https://www.worldclim.org/)[15]，而土壤数据来自于世界土壤数据库v1.2(https://daac.ornl.gov/)[16]。政府间气候变化专门委员会(Intergovernmental Panelon Climate Change，IPCC)综合考虑人口、经济、技术进步和资源利用等因素，提出了共享社会经济路径(Shared socioeconomic pathways，SSPs)来代表全球可持续发展路径、中等可持续发展路径、地方可持续发展路径、不平衡发展路径和常规发展路径[17]。研究选择4种SSPs(SSP126，SSP245，SSP370和SSP585)进行预测，它们分别代表碳排放由低到高的未来气候情景。在SSP126情景中，全球CO2排放量大幅减少，到本世纪末全球气温升高约1.8℃；在SSP245情景中，CO2排放量在本世纪中叶开始下降之前徘徊在当前水平，到本世纪末气温将上升2.7℃；在SSP370和SSP585情景中，CO2排放量将大幅度增加，到本世纪末平均气温分别上升3.6℃和4.4℃[17]。研究使用北京气候中心气候系统模式(BCC-CSM2-MR)作为全球气候模式(Global climate model，GCM)，因为已有研究证实该模式非常适合模拟中国的气候变化[18-19]。

由于环境变量之间的多重共线性会导致模型过拟合，影响模型的输出结果，所以使用ENMTools对环境变量进行相关性分析[13，20]。如果一组环境变量的相关性系数绝对值大于0.7，则根据刀切法(Jackknife)检验结果除去贡献较低的那一个[21-22]。

表1 研究使用的环境变量Table 1 List of environmental variables used in this study

1.3 模型优化及评价

研究使用MaxEnt 3.4.4软件预测醉马芨芨草的适生区分布[3]，随机选取75%的发生点数据用于训练模型，而剩下25%的数据用于检验。采用以往研究中常用的10倍交叉验证(Ten-fold cross-validation)确保模型的准确性[11]。

之前的研究表明特征类(Feature class,FC)和正则化乘数(Regularization multiplier,RM)的设置对MaxEnt模型预测的稳健性影响较大，且默认设置的参数通常不是最优的[23-24]。因此，研究采用R 3.6.3软件(https://www.r-project.org/)和Kuenm包(https://github.com/marlonecobos/kuenm)优化MaxEnt模型。Kumen包通过计算候选模型的偏ROC(Partial ROC,pROC)、遗漏率(Omission rates,OR)和AICc值综合考虑模型的统计显著性、预测能力和复杂性[26]。最终选择具有统计学显著性(pROC=0)，小于遗漏率阈值(OR<5%)，且AICc值最小的模型参数[26]。

AUC即受试者工作曲线(ROC)下的面积，是一种不依赖阈值的模型表现评估，常用于评价MaxEnt模型[11]。AUC值越接近1，则表示模型越可靠。通常0.6

1.4 适生区等级划分与变化比较

研究使用最大训练灵敏度加特异性逻辑(Maximum training sensitivity plus specificity logistic,MTSS)阈值对物种的“存在-不存在”进行划分。研究证实该方法在模型输出基于仅存在(presence-only)数据时，受“出现点-背景点”比例以及物种流行度影响较小，可降低流行度低的物种的遗漏误差和流行度高的物种的错判误差[27]，且具有客观性、等效性、判别力强等特点[11]。根据模型输出的逻辑值(logistic value)，即适生指数，将醉马芨芨草的生境划分为：非适生区(0～MTSS)、低适生区(MTSS～0.3)、中适生区(0.3～0.5)、高适生区(0.5～1)。

使用ArcGIS 10.7软件与SDMTools 2.5 (http://www.sdmtoolbox.org/)[28]来可视化当前和未来气候条件下醉马芨芨草适生区分布的差异。即根据MTSS阈值制作“存在-不存在”的二进制图层，并通过未来各气候情景适生区预测的二进制图层中减去当前适生区模拟的二进制图层观察变化。使用ArcGis 10.7进行地理数据分析，并使用Microsoft Excel和GraphPad Prism 8软件统计数据与绘制统计图。

2 结果与分析

2.1 发生点与环境变量筛选

经过筛选，最终使用229个醉马芨芨草发生点数据以及8个环境变量(bio2，bio15，bio19，elev，t_ece，t_oc，t_ph，tamx12)用于建模。这些环境变量之间的相关性系数绝对值均小于0.7(图2)。

图2 用于建模的环境变量相关性Fig.2 Correlation of environmental variables used to modeling

2.2 模型优化与评价

使用Kuenm包一共创建了1 160个候选模型，所有候选模型都具有统计学意义(pROC=0)，其中139个候选模型符合5%遗漏率标准(图3)。最终选择符合5%遗漏率并且AICc值最小的模型。最优模型pROC=0、遗漏率=3.51%、AICc=5 335.42，其参数为：RM=3.3，FC=TH。然后用该RM和FC参数运行MaxEnt，并进行10次交叉验证，平均AUC为0.955，标准偏差为0.011 7。这表明模型预测结果具有很好的准确性和可信度，可用于模拟醉马芨芨草的适生区分布预测。

2.3 影响醉马芨芨草适生区分布的关键环境变量

利用刀切法检验环境变量对醉马芨芨草适生区分布预测收益的重要性(图4)。仅使用12月平均最高气温(tmax12)时，所获得的增益值最高，其次是海拔(elev)、最冷季降水量(bio19)和降水变异系数(bio15)。这表明这些环境变量中含有很多有效信息。除去tmax12后，增益值减少最多，其次是bio19，elev和bio15。这说明这些环境变量中含有许多其他变量中不可替代的重要信息。同时，上述4个环境变量对模型的累积贡献率为96.8%(表2)。因此，tamx12,elev,bio19和bio15是影响醉马芨芨草分布预测的关键环境变量。

图3 所有候选模型的遗漏率与AICc值自然对数Fig.3 Omission rate VS the natural logarithm of AICc of all candidate models

图5为模型输出逻辑值与4个关键环境变量的响应曲线。4个关键环境变量与逻辑值之间都为单峰响应关系。tmax12，elev，bio19和bio15分别在—1.83℃，3 082.69 m，4.47 mm和97.65时，模型输出逻辑值达到峰值。根据MTSS来划分环境变量的适生阈值(表2)。tamx12，elev,bio19和bio15的适生阈值分别为—7.5～7.9℃，1 053.8～4 590.9 m，1.6～24.6 mm和25.5～111.5。其他环境变量的适生阈值如下：bio2(>10.09℃)，t_ece(>0 ds(m-1)，t_oc(>0%)，t_ph(>5.06)。

图4 刀切法检验环境变量对分布预测收益的重要性Fig.4 Effects of environmental variables on the gain of distribution prediction using Jackknife test

图5 模型输出逻辑值与关键环境变量的响应曲线Fig.5 Response curves between logistics output and key environmental variables

表2 关键环境变量的贡献率与适生阈值Table 2 Percent contribution and adaptive threshold of key environmental variables

2.4 当前醉马芨芨草适生区分布

在当前气候条件下，醉马芨芨草适生区分布模拟如图6所示。结果表明：醉马芨芨草总适生区面积约为2.33×106km2，约占中国大陆面积的24.24%，主要分布于新疆、西藏、四川、青海、甘肃、内蒙古、宁夏、陕西和山西等省区。高、中、低适生区面积分别为5.48×105km2，9.70×105km2和7.87×105km2，分别占中国大陆面积约5.71%，10.10%和8.42%。

2.5 未来醉马芨芨草适生区分布

根据未来气候数据预测了4种共享社会经济途径情景下2050s和2070s醉马芨芨草适生区分布(图7)，并对不同时期和气候情景各等级适生区占中国大陆面积比例与总适生区平均海拔进行统计分析(图8)。图8E显示，在SSP245，SSP370和SSP585情景下，适生区面积逐渐增加，其中SSP370情景增加最明显；在SSP126情景下，2050s增加而2070s减少，但是2070s总适生区面积相对于当前是增加的。图8 A,B,C和D提示，在4种未来气候情景下，高适生区面积变化不大，中适生区面积有小幅度减小，低适生区面积有较大幅度的增加。在4种未来气候情景下，总适生区平均海拔皆高于当前，至少增加118.75 m(图8 F)。研究将当前与未来不同气候条件下醉马芨芨草适生区空间分布进行对比(图9)。结果显示，在未来不同情景下山西地区适生区有很明显的缩小，其次是新疆西南部；在新疆西北部、西藏中部、青海中部、内蒙古东北部以及黑龙江和吉林的部分地区有明显的扩张。总体上在东南和西南方向缩减，向西北和东北方向扩张。

图6 当前气候背景下醉马芨芨草适生区分布模拟Fig.6 Simulation of distribution of A. inebrians suitable area under the current climate注：审图号GS(2022)1873，下同Note:On drawing No. GS(2022)1873,the same as below

图7 未来(2050s,2070s)4种共享社会经济途径气候背景下醉马芨芨草适生区分布预测Fig.7 Prediction of suitable area of A. inebrians under four shared socioeconomic pathways scenarios in the future (2050s,2070s)

图8 未来气候背景下各等级适生区面积占比(A,B,C和D)；各时期总适生区面积占比(E)；各时期总适生区平均海拔变化(F)Fig.8 Proportion of suitable areas at different levels under future climate background (A,B,C and D)；proportion of total suitable area in each period (E);average elevation change of total suitable area in different periods (F)

图9 当前与未来不同气候条件下醉马芨芨草适生区变化对比Fig.9 Changes in the area of suitable habitat of A. inebrians under different future climate conditions relative to that under current climate conditions

3 讨论

本研究利用优化后的MaxEnt模型和土壤、地形、气候和生物气候等环境变量，模拟和预测了在当前和未来四种气候情景下醉马芨芨草适生区的分布。最终模型的遗漏率、AICc和AUC分别为3.51%，5 335.42和0.955，环境变量的响应曲线较为平滑。这表明我们的模型输出结果是稳健可靠的。

温度、降水、海拔等环境因素影响植物的生长发育，同时也制约着植物物种的分布[29]。刀切法检验结果表明，影响醉马芨芨草分布的关键环境变量是12月平均最高气温(tmax12)、海拔(elev)、最冷季降水量(bio19)和降水变异系数(bio15)。这些环境变量对最终模型的累积贡献达到了96.8%。利用MTSS对这些环境变量的适宜性阈值进行了划分。各关键变量的适宜范围是：tmax12(—7.5～7.9℃)，elev(1 053.8～4 590.9 m)，bio19(1.6～24.6 mm)和bio15(25.5～111.5)。这与之前报道醉马芨芨草主要生长在高山及亚高山草原，海拔1 700～4 200 m的山坡、草地、田边等，具有抗寒、耐旱等生物特性基本相符[6,30]。此外，陈雅琦等[31]证实醉马芨芨草中的内生真菌使其具有良好的耐盐碱性。岳永寰等[32]发现醉马芨芨草种子在pH 6～10范围中有良好的发芽率(>90%)，提示其对极端环境有极强耐受性。在本研究中，表层土壤盐度(t_ece)适生阈值>0 ds·m-1，在盐度持续升高的情况下依然保持很高的适生指数(0.66)；表层土壤PH(t_ph)适生阈值>5.06，pH 6～10范围时有较高的适生指数(0.38～0.65)。即醉马芨芨草有很强的的盐碱耐受性，与前面的研究一致。

在目前的气候条件下，醉马芨芨草适生区主要分布在新疆、西藏、四川、青海、甘肃、内蒙古、宁夏、陕西和山西等省区，约占中国大陆面积24.24%。其中高适生区主要分布在新疆、青海、甘肃、宁夏、内蒙古、西藏以及四川西部。模拟结果与已知醉马芨芨草地理分布基本一致[6]。在未来4种气候变化情景下(2050s，2070s)，醉马芨芨草适宜生境的分布面积有所增加。这可能与其强大的环境适应能力有关。在SSP245，SSP370，SSP585情境下，醉马芨芨草适生区呈增加趋势，SSP370最为明显，SSP245和SSP585次之。在SSP126情景下，醉马芨芨草适生区先增加后减少，但仍然比当前适生区面积要多。这可能是因为适当增加温度可以促进植物的扩张，而持续升温则会带来一定负面影响[33]。预测未来醉马芨芨草适生区平均海拔将增加，并且向西北和东北方向扩张。这与之前的研究观点类似。例如，Lenoir等[34]发现随着全球变暖气候发生变化，大多数植物的地理分布将会向高纬度和高海拔地区转移。李晓辰[35]预测大白刺(Nitrariaroborowskii)、霸王(Zygophyllumxanthoxylon)、沙拐枣(Calligonummongolicum)等多数荒漠植物的适生区在未来有向北增加趋势，且分布重心也将向高纬度地区迁移。段义忠等[36]预测四合木(Tetraenamongolica)未来适生区呈现破碎化分布，并有向东北、西北高纬度地区迁移的趋势。

但值得注意的是，本研究存在一定局限性。因为生物因素、人为因素、植被覆盖度等同样会影响植物生境适宜性与传播[37-39]。例如，家畜可促进醉马芨芨草的传播，尤其是羊更容易增加其种子的扩散强度[40]。此外，风媒介和水媒介也是醉马芨芨草种子扩散的途径[41]。人为的过度放牧使草地退化，使醉马芨芨草在竞争中占有较强优势，从而加速其扩增[42]。因此，在实际应用中应综合考虑模拟结果的不确定性，才能更有效地预测未来气候变化条件下醉马芨芨草适生区分布格局，从而对这种毒害草的防控与资源化利用做出合理的决策。

4 结论

综上所述，利用优化后的MaxEnt模型能较好的预测醉马芨芨草分布。12月平均最高气温(tmax12)、海拔(elev)、最冷季降水量(bio19)和降水量变异系数(bio15)是主导其分布的关键环境变量。随着未来气候变化，醉马芨芨草适生区面积将增加，并向高纬度和高海拔地区迁移。本研究为醉马芨芨草的防治和资源利用、草地生态学研究和草地放牧管理提供理论依据和参考。