基于MaxEnt模型青藏高原老芒麦适生区模拟预测

韩梦丽， 白史且， 孙盛楠， 鄢家俊， 张传杰， 张昌兵， 张靖雪， 游明鸿， 李达旭， 严学兵*

(1.扬州大学动物科学与技术学院,江苏 扬州 225000；2.四川草原科学研究院，四川 成都 610000)

老芒麦(ElymussibiricusL.)别名西伯利亚野麦草，是禾本科(Gramineae)小麦族(Triticeae Dumort)披碱草属(Elymus)多年生疏丛型中旱生植物，是披碱草属的模式植物[1],它是分布于北半球寒温带的一种比较古老的野生草种[2-3]。老芒麦抗寒性强，有的老芒麦种质资源在零下四十度的低温和海拔4 000 m左右的高寒地区仍旧能够安全越冬[4]。我国是披碱草属的主要分布区之一，其主要分布于西北、华北、东北以及西南的四川和西藏[5]，尤其是在青藏高原地区。该属植物不但是人工草地建植及高寒草甸、沙化地等脆弱生态区生态恢复的重要物种，也能为小麦(Triticumaestivum)、大麦(Hordeumvulgare)等麦类作物的分子育种提供重要的基因资源[5-6]。

从IPCC近130年以来的报告显示来看，全球地表平均温度已上升了0.85℃，未来地表温度还将呈现持续上升的趋势[7]，全球气候的变化不仅会改变降水格局，同时也会影响物种的分布格局[8]。现代物种的分布在很大程度上受第四纪冰期和间冰期交替出现的影响[9]。具有“世界屋脊”之称的青藏高原是世界上海拔最高、面积最大和地质发育最年轻的高原，青藏高原地区的高寒植物对气候变化的响应十分敏感，它们的进化往往和气候变化、人类活动、生命的进化等密切相关[10-11]，而且由于青藏高原特殊的地质、地理分布状况，一直是开展生态学的天然理想实验室[12]。因此通过预测气候变化下物种潜在的分布格局，再结合其它方面(如谱系地理学、遗传多样性等)的研究，阐明青藏高原老芒麦的进化和起源，提出青藏高原老芒麦种质资源的保护和利用策略，还可为青藏高原高寒植物的保护提供宝贵经验。

目前，国内外常用生态位模型主要有基于生物气候数据的Bioclim和Domain模型、生态位因素分析模型(Ecological niche factor analysis,ENFA)、遗传算法模型(Genetic algorithm，GARP)及最大熵值模型(Maximum Entropy,MaxEnt)等[13]。最大熵模型(MaxEnt)作为机器学习算法模型的一种，是目前应用最广泛、预测效果较好的物种分布模型之一[14-15]。MaxEnt模型相比较其他模型，操作更加简单、预测结果更加准确，已经广泛应用于气候变化下物种的潜在分布的预测。根据毛志远等[16]的研究来看，运用MaxEnt模型和GARP模型对石蒜(Lycorisradiata)适生区预测时，Maxent软件的预测结果更为准确；陈积山等[17]基于MaxEnt模型开展了羊草(Leymuschinensis) 物种生境适宜性研究；武自念[18]采用MaxEnt预测气候变化背景下扁蓿豆(Medicagoruthenica)的适宜性生境分布区域及影响其分布的主导气候因子;郭斌等[19]研究气候变化背景下川西北高原多年生垂穗披碱草(Elymusnutans)种植适生区分布。另外，在分布点数据有限的情况下，MaxEnt模型比其它模型预测结果更好、准确性更高，更加适合小样本数据的预测[20]。根据陈新美等[21]对4个物种34个不同样本量对MaxEnt模型预测物种分布精度和稳定性的影响研究，发现样本量的大小对MaxEnt模型预测物种空间分布的精度影响不大。本研究运用MaxEnt模型和ArcGis软件对老芒麦在末次间冰期(LIG)、末次冰期冰盛期(LGM)、当前时期(Current)和2080年(Future)4个时期在青藏高原的分布格局进行模拟预测，结合分布格局的变化和影响老芒麦在青藏高原分布的环境因子进行分析，推测产生这种分布格局的成因，以期为青藏高原种质资源的收集、保存和利用提供有效的理论基础。

1 材料与方法

1.1 数据来源

研究材料选择在青藏高原广泛分布的披碱草属模式植物—老芒麦，分布数据分别来自野外采样和数据库，包含青海、四川、甘肃、西藏四个地区的55个采样点(表1)，基本涵盖了老芒麦的有效分布范围，另外从全球多样性信息网GBIF (http://www.gbif.org/) 获得了33个样点信息，将这88个地理分布数据的经纬度以及物种名称输进Excel中，并转换成csv格式进行保存，以便接下来的数据分析。

环境数据来源于Worldclim(http://www.worldclim.org)的20个环境变量(表2)。这些环境因子包括末次间冰期(Last Interglacial，120—140千年)、末次冰期冰盛期(Last Glacial Maximum，21千年)、当前(Current，1950—2000年)及2080年(Future)4个时期[8]。

1.2 模型软件使用方法

本试验主要运用了MaxEnt version 3.4.1，ArcGis10.2软件。

1.2.1环境变量数据的处理 为了提高模型预测的准确度，避免模型的过度拟合，需要对环境数据进行预处理。分析20个环境变量之间的相关性和贡献率，避免自相关，去掉相关性高且对模型预测结果贡献率小的变量[22]。首先使用20个环境变量初步建立模型，根据MaxEnt预测结果中的刀切法(Jackknife，MaxEnt自带的模块，其结果可以反映不同的环境变量对预测分布的贡献大小。该方法计算了“仅此变量”、“除此变量”和“所有变量”3个状态进行模拟时的训练得分，当“仅此变量”得分较高时，说明这个因子对物种分布的贡献较大，预测能力较高，当“除此变量”得分降低较多时，说明该变量拥有较多的特别的信息，对物种的分布较为重要[23])。检验20个环境变量对预测结果的贡献率，再通过ArcGis中的波段集统计(Band collection statistics)进行相关性检验，最后从相关性系数绝对值大于0.8的2个环境变量之间筛选出贡献率大的环境因子用作再次建模[24]。

表1 老芒麦88个居群分布点的地理位置

表2 用于Maxent的环境变量

1.2.2MaxEnt模型的使用 将88个csv格式的物种分布数据和筛选过后的环境变量(经过ArcGis转换为asc格式的文件)导入MaxEnt模型中的Sample和Environmental layer模块中。设置好参数[15]:勾选create response curves，do jackknife to measure variable importance等选项，Replicated run type设置为Bootstrap；测试集设定为25%，训练集为75%，重复次数为10次，点击Run。

1.2.3MaxEnt模型精度 ROC曲线被公认为是评判模型精确度的最优指标[25]。以假阳性率(1-特异度) 为横坐标，真阳性率(灵敏度)为纵坐标绘制而成的曲线即为ROC曲线，ROC曲线下的面积为AUC值。当ROC曲线的斜率为45°的对角线时，也就是AUC值为0.5，表示该模型为随机模型，结果的精确度为50%；AUC值大于0.5时，表示模型的输出结果要好于随机模型；AUC值小于0.5时，表示模型的输出结果不如随机模型。一般认为，AUC值为0.5～0.7时模型结果一般；AUC值为0.7～0.9时模拟结果很好；AUC大于0.9时模拟的结果最好最精确[26]。

1.2.4ArcGis中绘制适生分布图 将MaxEnt预测结果中的asc格式的文件通过ArcMap的conversion工具转换为栅格格式，点击工具箱中的重分类工具对老芒麦的适生区进行划分。ArcGis软件中等级划分的方法很多，根据崔相艳等[27]的研究发现：Standard deviation对高适生区划分的范围过窄，Geometrical interval以及Quantile对高适生区划分的范围过宽，Manual默认为10个等级,划分过细，Equal interval采用等间距划分，过于人为。而自然间断法(Natural breaks)分类可以很好地“物以类聚”，类别之间的差异明显，而类内部的差异是很小的，每一类之间都有一个明显的断裂之处，考虑了数据分布的统计特征比较自然，近年来被广泛使用。综上，本研究运用Natural breaks对物种的适生区等级进行划分，分为非适生区、低适生区、中适生区和高适生区4个等级，插入合适的图例、指北针和比例尺等，得到老芒麦4个时期在青藏高原的适生分布图。采用ArcGis软件工具箱的Spatial Analysis Tools的Zonal Statistic的分区几何统计计算适生区的面积[28]。

2 结果与分析

2.1 关键环境变量

基于当前时期的分布，利用刀切法分析环境变量对老芒麦适生区贡献大小的研究结果如图1所示。对老芒麦分布贡献较为重要的环境变量有：年均温(Bio1)、最冷季均温(Bio11)、最暖季均温(bio10)、最冷月最低温(Bio6)、海拔(Elev)、最干季均温(Bio9)和年降水量(Bio12)等；再结合相关性分析结果以及环境变量贡献率，最终筛选出8个关键环境变量，分别为：年均温(Bio1)、昼夜温差月均值(Bio2)、等温性(Bio3)、年温差(Bio7)、最冷季均温(Bio11)、年降水量(Bio12)、降水量变异系数(Bio15)和海拔(Elev)。最终使用这8个环境变量用于本研究最终的建模和后续分析，其余3个时期关键环境变量的筛选同当前时期的操作方法一致。

2.2 MaxEnt模型精确度的验证

基于当前时期的结果，通过MaxEnt模型输出的ROC曲线表明(图2)，MaxEnt模型的AUC值为0.944，根据AUC评价指标，MaxEnt模型的预测结果准确度极好，说明该模型的预测结果较为可靠。

图1 刀切法检验的环境变量

图2 MaxEnt模型模拟的ROC曲线

2.3 老芒麦在青藏高原的潜在适生区

基于MaxEnt模型的预测结果，得到老芒麦在青藏高原4个时期的适生区分布(图3)和适生面积(表3)。由图3和表3可知，在当前时期老芒麦在青藏高原主要分布在祁连山、巴颜喀拉山脉以及横断山脉等地区。青藏高原的总面积约为250万km2，从4个时期物种分布的变化看，末次间冰期(图3B)到末次冰期盛冰期(图3A)，总适生面积减少，但是中适生区和高适生区面积却有所增加，分别增加1.2%和1.5%，适生区域逐渐向青藏高原的东部边缘低纬度地区扩散，藏北高原腹地的低适生区面积发生退缩，而喜马拉雅区域西部则出现老芒麦的踪迹，说明这里出现了适宜老芒麦生存的生态环境。比较末次冰期(图3B)和当前时期(图3C)的物种分布，发现老芒麦整体分布呈现退缩趋势，尤其是高适生区，羌塘以南藏南高原适生区分布减少，雅鲁藏布江中高适生区减少。比较当前时期(图3C)和2080(图3D)年，发现高适生区面积几乎没有变化，仅低适生区和中适生区面积分别增加0.6%和1.6%，整体的分布向青藏高原东部边缘低纬度地区移动。

表3 四个时期不同适生区的面积

图3 基于MaxEnt预测的末次间冰期(LIG)、末次冰期盛冰期(LGM)、当前时期(Current)和2080年老芒麦在青藏高原的分布

图4 环境变量重要性刀切法检验

2.4 老芒麦的适生区分布与环境变量关系

利用刀切法检测当前时期8个关键环境变量对物种潜在分布的贡献大小，结果如图4所示。年均温(Bio1)、最冷季平均温(Bio11)是影响老芒麦适生区分布的关键变量，其次分别是年降水量(Bio12)、海拔(Elev)、降水量变异系数(Bio15)、年温的变化范围(Bio7)、等温性(Bio3)，而昼夜温差月均值(Bio2)贡献率不明显。结合贡献率分析，年降水量(Bio12)、年均温(Bio1)、海拔(Elev)昼夜温差月均温(Bio2)、的贡献率分别为38.1%，27.5%，15.3%和7.3%，累计贡献率达到了95.3%。综上分析，影响老芒麦现代分布的主要环境因子为年均温(bio1)、最冷季平均温(Bio11)、年降水量(Bio12)和海拔(Elev)。

利用MaxEnt模型绘制出主导老芒麦分布概率的四个环境因子的响应曲线(图5)，由图可知，老芒麦适生的年均温(Bio1)范围在—16.8～24.1℃，最适值为3.5℃。最冷季平均温(Bio11)的范围是—28.39～18.16℃，最适值为—3.2℃。年降水量(Bio12)的范围是20～3 033 mm，最适值为596.2 mm。海拔(Elev)的适值范围为105～6 406 m，最适值为3 240 m。

图5 老芒麦分布对4个环境因子的响应曲线

3 讨论

许多学者认为，现代物种的分布很大程度上是受第四世纪末末次冰期盛冰期的影响[29]。青藏高原由于其地势复杂，形成了许多山脉河谷，是许多植物的生物避难所。从本研究的预测结果上看，在冰期，老芒麦的分布区域发生扩张，老芒麦的冰期避难所之一很可能在青藏高原的喜马拉雅区域的西部邻近昆仑山的地方。冰期结束后，适生区面积进一步减少。未来60年内老芒麦分布整体向东部区域扩张，分布面积进一步减少，但和当前时期相比，变化并不明显，呈现稳定的状态。

随着全球变暖的加剧，青藏高原的生态系统受到显著影响，老芒麦在青藏高原分布也将受到更大的影响，生物避难所极可能发生转移。胡忠俊等[12]对青藏高原紫花针茅物种分布格局模拟发现羌塘以南藏北核心区域及喜马拉雅山西部阿里地区是紫花针茅的两个可能避难所。刘文胜等[18]发现羌塘高原和可可西里地区是现代青藏苔草分布最适宜的区域。李宁宁等[31]研究了青藏高原云杉在气候变化情景下未来的分布，发现在未来气候变暖的背景下有利于紫果云杉的扩张。本研究中老芒麦在4个时期的分布格局的变化也有相似的规律。从冰期到现在再到未来，老芒麦的分布发生先收缩再扩张的变化，在未来气候变暖的气候背景下，反而有利于老芒麦的生存，但高适生区和中适生区变化不大，比较稳定，整体分布向青藏高原东南部的横断山脉移动。青藏高原东部边缘地区山体高大、地形复杂，是典型的山地气候，这一地形优势有利于老芒麦的分布。本文对未来气候变化下预测老芒麦在青藏高原的分布的研究将有利于将来为老芒麦的保护提供一定的理论基础。

据MaxEnt分析结果显示，温度(年均温和最冷季平均温)、年降水量和海拔是影响老芒麦当前分布的主要环境变量。由图4可知，年均温的和最冷季平均温的正规化训练增益(Regularized trainning gain)最大，其次才是年降水量，由此可推论，在当前气候条件下温度要比降水对老芒麦的影响更大，年均温是影响老芒麦分布的最关键环境变量。曹文侠等[30]的研究表明披碱草属不同材料表型性状与年均温具有明显相关性，年均温对披碱草属不同材料表型性状的遗传变异有重要的影响。根据《中国气候变化蓝皮书》[31]，自2016年以来，青藏高原地区降水量持续异常偏多，积雪覆盖区的积雪覆盖率是2002年以来最高，祁娟[32]对披碱草属植物野生种质资源生态适应性研究发现降水量对披碱草属的遗传多样性有显著影响。而从本研究结果中可以得知，年降水量为596.2 mm时，最适合老芒麦的生长，这为披碱草属植物的生态遗传学研究提供了很好的环境变量。Yan等[33]发现披碱草属植物种群遗传多样性随着海拔呈现先升高后降低的趋势，最大值出现在3 300 m左右。由图5可知，适合老芒麦生长的海拔范围在105～ 6 406 m，老芒麦的最适值为3 240 m，这为高山植物在最适海拔区域具有更高的遗传多样性这一假说提供了证据。这些研究发现很好的证明了本文所筛选的关键环境变量(年均温、最冷季平均温、年降水量和海拔)具有现实意义。另外，我们可以根据老芒麦的扩张路线以及避难所的移动，推测青藏高原东缘地区将会是老芒麦分布的中心，这有助于我们更加快捷有效的对老芒麦种质资源提出利用和保护策略以及理解青藏高原植被的起源和进化问题。

4 结论

本研究表明，影响老芒麦分布的主导环境变量是年均温(Bio1)、年降水量(Bio12)、最冷季平均温(Bio11)以及海拔(Elev)。气候的改变引起青藏高原降水格局和生物小气候的改变，对老芒麦的潜在分布及变化趋势产生了一定的影响，生物避难所会发生迁移。2080年老芒麦的适生区发生扩张，适生区有向青藏高原东南部边缘横断山脉移动的趋势。本研究可以为青藏高原种质资源的收集、保存和利用提供有效的理论基础，还可为青藏高原高寒植物的保护提供宝贵经验。