全世界约有蝗虫12000余种, 除南极洲外各大洲均有分布。飞蝗(Locutsa migratoria L.)是世界范围内的重要害虫, 已知有1个种10个亚种, 其中西藏飞蝗(Locusta migratoria tibetensis Chen)为青藏高原特有且分布海拔最高的亚种^[1], 主要分布在西藏、四川及青海等省区^[2], 主要以青稞、牧草及杂草为食^[3-5]。据史料记载, 西藏飞蝗曾在青藏高原多次成灾, 蝗灾所到之处庄稼颗粒无收。20世纪90年代以来, 西藏飞蝗在四川甘孜州和阿坝州, 西藏普兰、昌都、阿里及青海玉树等地曾大面积爆发, 尤以2003—2006年该虫发生极为严重, 特别是雅鲁藏布江、雅砻江以及金沙江等横断山脉河谷地带, 蝗灾的爆发对当地农牧业正常生产造成极大破坏, 导致难以估量的经济损失^[6]。

预测物种适生分布是生态学的重要研究领域, 物种分布模型(species distribution models, SDM)逐渐发展并成为研究物种适生性的一个重要工具^[7]。此类模型根据物种已知的“存在”或“不存在”数据, 结合相应的环境变量, 利用一定的数学理论计算物种在目标地区的分布概率, 量化环境变量与物种分布之间的关系, 寻找物种分布的限制因子及生境偏好等。目前, 物种分布模型的常用领域主要为：预测物种潜在分布区、物种分布与气候关系、濒危及珍贵物种适生性预测和古生物地理研究等方面^[8-10], 并起到了良好的效果^[11]。目前常用的物种分布模型中, 最大熵模型MaxEnt(maximum entropy model)的模拟精度要高于其他模型, 加之该软件具有运行时间短、操作简便、运行结果稳定和所需样本量小等特点在业内得到了广泛的认可^[12-13], 已被应用于多种病虫害适生性的模拟预测, 如柑橘黄龙病菌(Candidatus Liberibacter asiaticus)、叶缘焦枯病菌(Xylella fastidiosa)、柑橘木虱(Diaphorina citri Kuwayama)和稻水象甲(Lissorhoptrus oryzophilus Kuschel)等^[14-27], 且模拟的结果与物种实际分布的吻合度高。

目前对西藏飞蝗的研究主要集中在生物学特性、环境胁迫、数量性状及遗传特征等方面^[28], 还未发现对其生态位预测的相关研究报道。本研究基于MaxEnt模型, 研究西藏飞蝗分布与环境变量之间的联系, 揭示影响其分布的主导气候因子, 分析西藏飞蝗在中国的潜在分布区, 旨在为制定合理的防控措施提供重要参考和理论依据。

1 材料与方法 1.1 西藏飞蝗的分布数据与处理

MaxEnt模型模拟需要物种分布数据和环境数据, 分布数据即为“存在”点的经纬度信息。西藏飞蝗的分布数据通过以下方式获得：(1)实地考察, 四川境内的采样点多为实地考察获得, 调查点主要位于西藏飞蝗高发区的甘孜州和阿坝州, 发现西藏飞蝗若虫或成虫即记为“存在”, GPS记录经纬度。(2)查询物种分布数据库。本研究查询的数据库包括“国际农业与生物科学中心(CABI, http://www.cabi.org/)”数据库、“全球物种多样性信息库(GBIF, http://www.gbif.org/)”、中国西南地区动物资源数据库(http://www.swanimal.csdb.cn)、教学标本资源共享平台(http://mnh.scu.edu.cn/)。(3)检索国内外公开发表的西藏飞蝗相关的期刊论文。数据处理：采用缓冲区分析法对获得的分布点进行校对、筛选, 排除空间关联性较大造成的过拟合模拟的影响, 由于环境变量空间分辨率为2.5 arc-minutes(约4.5 km²), 故设置缓冲半径为1.5 km, 当分布点之间距离小于3 km时, 只保留其中一点, 共获得分布点47个, 通过Google earth确定每个分布点的经纬度。参照MaxEnt软件操作手册, 将经纬度坐标输入Excel, 保存格式为*.CSV, 字段包括物种名称、经度和纬度^[29]。由于西藏飞蝗生活在交通不便、气候恶劣的高海拔地区, 所获得的分布点很难完全覆盖其实际分布。研究表明, 在样本量很小的情况下, MaxEnt模型较其他模型预测效果更好、准确性更高, 更适用于小样本数据的预测。

1.2 预测环境变量 1.2.1 环境变量数据来源

本研究所用的环境因子均来源于Worldclim气候数据库, 数据是以全球的气象记录信息为基础数据, 整合插值生成的全球气候栅格数据^[30]。利用ArcGIS软件提取研究区(西藏、四川、云南、贵州和青海)19个生物气候因子(具有较强的生物学意义, 可反映温度与降水的特点及季节性变化特征)、月平均降水量、月平均最高及最低气温等气候指标数据。另外还选取了海拔数据作为地形因子(表 1)。上述数据的空间分辨2.5 arc-minutes(约4.5 km²)。

1.2.2 影响西藏飞蝗潜在分布关键环境变量的筛选

影响物种分布的环境因子多种多样, 包括气候因子、土壤因子、植被分布等, 本文中选取了worldclim提供的24个环境变量构建初始模型。初始模型构建过程中选择MaxEnt软件中的刀切法(jackknife test)检验来测定环境变量对模型预测的贡献大小, 以此剔除对MaxEnt模型预测结果贡献较小的环境变量, 对贡献较大的关键限制因子进行Spearman相关分析, 提高模型模拟的精度, 在此基础上重建西藏飞蝗分布的最大熵模型, 并对模拟结果进行准确性评价^[31]。

刀切法检验使用MaxEnt软件, Spearman相关分析使用Spss分析软件。

1.3 MaxEnt模型的操作与结果评价

MaxEnt模型来源：选用基于最大熵理论的MaxEnt模型软件来预测西藏飞蝗在中国的适生区分布, 该软件可在MaxEnt主页(http://www.cs.princeton.edu/~schapire/maxent/)上免费获取^[32], 当前最新版本为3.3.3 k版。

MaxEnt模型的使用：将上述1.1中47个西藏飞蝗分布点数据和1.2中的环境因子数据添加到软件中, 随机选取75%的西藏飞蝗分布点作为训练集(training data)建立预测模型, 剩余25%的西藏飞蝗分布点作为测试集(test data)验证模型, 选择刀切法(jackknife test)测定各变量权重, 选择创建环境变量响应曲线, 其余参数均选择模型的默认值^[33]。

本研究采用受试者工作特征曲线(Receiver Operating Characteristic Curve, ROC曲线)下面积值, 即AUC值(Area Under Curve, AUC)为判据来评价模型模拟结果。ROC曲线, 又称为感受性曲线, 该曲线根据一系列不同的二分类方式, 分别以假阳性率(1-特异率)和真阳性率(1-遗漏率)为横纵坐标绘制而成, 曲线下面积值即AUC值, 取值范围为[0, 1], AUC值越接近1表明环境变量与分布模型之间相关性越大, 预测结果精度越高^[34]。ROC曲线可使用Spss分析软件获得, 本研究中MaxEnt模型可自动生成ROC曲线, 并给出相应的AUC值。具体评价标准^[35-36]见表 2。

1.4 西藏飞蝗适生等级分区

MaxEnt模型输出结果为西藏飞蝗在世界范围内的存在概率, 数据为ASCⅡ格式, 首先使用ArcGIS的ArcToolbox的格式转换工具, 将该数据转为Raster格式, 使该结果可在ArcGIS中显示, 而后使用“提取分析”功能得到该虫在中国的存在概率分布图^[37-38]。MaxEnt软件模拟输出的结果值在0—1之间, 值越接近1表示物种越可能存在。参考IPCC报告^[39]关于评估可能性的划分方法, 结合西藏飞蝗的实际情况, 利用“Reclassify”功能, 划分分布值等级及相应分布范围, 并使用不同颜色表示, 划分标准为：存在概率＜0.05为不适生区；0.05≤存在概率＜0.33为低适生区；0.33≤存在概率＜0.66为中适生区；存在概率≥0.66为高适生区。

2 研究结果 2.1 主导环境因子筛选及模型检验 2.1.1 影响西藏飞蝗潜在分布主导环境变量筛选

MaxEnt模型预测过程中需使用物种空间分布点数据和环境变量数据, 并不是每个环境变量对物种潜在分布预测都是必需的, 因此需要剔除对MaxEnt模型预测结果贡献较小的环境变量, 筛选得到贡献较大的关键限制因子^[40]。本研究参考雷军成等^[41-42]方法, 利用MaxEnt软件本身提供的模块, 即刀切法测定变量的重要性, 通过比较每个环境变量对模型结果的影响程度来剔除影响较小的变量, 最终获得对模型模拟贡献率较大的变量作为主导变量。本研究中, 依据MaxEnt的运算结果, 不同环境变量对西藏飞蝗潜在分布的累积贡献率见表 3。海拔高度为限制西藏飞蝗分布的主要变量, 单因子贡献率达到了52%, 其次, 8月份平均雨量、1月份平均雨量、等温性、1月份平均温度、12月份平均温度和2月份平均温度对西藏飞蝗潜在分布的累积贡献率为42.6%。

本研究使用Spearman相关系数法计算7个环境因子之间的相关性, 以消除共线性对模型建模过程和结果解释的影响, 环境因子之间的相关系数绝对值大于0.8^[43-44], 则认为两者具有很强的相关性。由表 4可以看出, 7个环境因子两两之间相关系数均未超过0.8, 因此本研究选定上述7种环境变量作为影响西藏飞蝗分布的主导变量, 在此基础上重建西藏飞蝗在中国的分布的最大熵模型, 并进行模拟结果的精度评价, 以获得最优秀的预测模型。

2.1.2 ROC曲线和AUC值对模型准确性检验

图 1是MaxEnt模型模拟输出的ROC曲线图, 图中表明基于所有环境变量和主导环境变量构建的西藏飞蝗在中国的地理分布模型的AUC值分别为0.996和0.993, 依据表 1的AUC值评价标准表, 本次构建模型的预测准确性达到“极好”的标准, 说明此结果可用于研究西藏飞蝗的适生区划。

2.2 西藏飞蝗在中国的地理分布的预测

根据7个主导环境变量与西藏飞蝗的分布数据构建的最大熵模型模拟结果, 按照高适生区、中适生区、低适生区和不适生区4个等级标准划分, 最终获得西藏飞蝗在中国的适生区划图(图 2)。预测结果显示, 西藏飞蝗在当代气候条件下预测的适生区范围为88.79°—113.57°E, 21.52°—38.33°N。其中高适生区主要位于四川和西藏两省区, 包括四川的甘孜州, 西藏的昌都地区、林芝地区、山南地区和拉萨市。中适生区沿高适生区周围分布, 主要集中青藏高原东部地区, 包括四川西部、西藏东部、云南北部等地。

西藏飞蝗当代适生区总面积为170.55×10⁴ km², 占中国国土面积的17.77%。其中, 高、中、低适生区面积分别为22.51×10⁴ km²、63.04×10⁴ km²和85.01×10⁴ km², 分别占总适生区面积的13.2%、36.96%和49.84%。四川和西藏高适生区面积最大, 分别为10.19×10⁴ km²和11.05×10⁴ km²(表 5)。

2.3 西藏飞蝗地理分布与环境变量之间的关系 2.3.1 刀切法检验结果

刀切法检验结果可反映不同环境变量对于分布增益的贡献大小, 该方法分别计算“仅此变量”、“除此变量”和“所有变量”模拟时的训练得分情况, 判定标准为：“仅此变量”时得分较高, 说明该因子具有较高的预测能力, 对物种分布贡献较大；“除此变量”时模型训练得分能力降低较多, 说明该变量具有较多的独特信息, 对物种分布较为重要。图 3为使用MaxEnt软件自带的刀切法分析各环境变量对西藏飞蝗分布重要性的结果图, 由图中可见海拔(Alt)是影响西藏飞蝗分布的关键环境变量, 其训练增益超过2.5；等温性(Bio3)也是较为重要的环境变量, 其单独训练得分超过1.5；另外5种主导变量的重要性排序为2月份平均温度(Tmean2)＞1月份平均温度(Ttmean1)＞12月份平均温度(Tmean12)＞1月份平均雨量(Prec1)＞8月份平均雨量(Prec8), 这5种环境变量的训练得分也均超过1.0。

2.3.2 环境变量对MaxEnt模型预测的影响

图 4是MaxEnt模型绘制的主导环境变量与分布概率之间的反馈曲线, 该图可反映不同阈值下环境变量的取值范围。参照1.4中的分级方法, 本文以0.33为阈来划分适合西藏飞蝗分布的环境变量的范围。结果显示, 海拔的适宜范围为2870—5188 m, 最适高度为3941 m, 海拔2870—3941 m时, 分布概率随海拔升高而增加, 3941—5188 m时, 随海拔升高而降低。等温性的适宜值为38.5—47.2, 最适値为42.9, 在38.5—42.9, 随等温性上升分布概率增大, 在42.9—47.2时, 随等温性升高分布概率降低。1、2、12月份平均温度和8月份平均雨量与分布概率之间的反馈曲线的变化趋势与上述2种变量类似, 只是在变化幅度及取值范围存在差异, 具体见表 6。在适宜范围内, 各主导环境变量的变化对西藏飞蝗的存在概率均有一定的影响, 而取值在适宜范围之外时对该虫分布影响的概率逐渐减小。

3 讨论

研究表明^[45], 环境变量的选择会对生态位模型的预测结果造成一定的影响。应用MaxEnt模型预测物种地理分布的很多研究, 采用世界气候数据库(http://www.worldclim.org)的Bioclim数据集中的19个生物气候因子作为全部或主要环境变量。19个环境变量是以温度和雨量为基础数据, 根据不同需求计算演生而来, 因此这些变量之间存在不可避免的自相关及多重线性重复等问题^[46-47]。研究业已证实, 这些高度相关的变量在模型预测过程中会引入冗余信息, 影响预测结果。在建模过程中, 为避免上述问题, 应首先对环境变量进行相关分析和有效筛选。本研究中, 参考孙敬松等^[48-50]方法, 以各因子对物种分布的贡献率来评价变量的重要性, 有效剔除了贡献率较小的环境变量并对筛选的变量进行了相关性分析, 获得了7种主导环境变量并对模型进行重建, 减少了冗余信息对模拟结果的影响, 提高了预测结果的准确性。

目前对模型精度评价应用最广泛的方法为ROC曲线法(即AUC法)^[51], 由于AUC不受诊断阈值的影响, 并可提供所有阈值范围上的性能评价结果, 因此目前被公认为生态位模型评价领域的最佳评价指标。AUC的取值范围在0.5—1, 值越接近于1模型精度越高。本研究中基于主导环境变量的模拟的训练集AUC平均值为0.993(非常接近1), 预测结果达到“极好”水平, 说明此次模型预测的地理分布与西藏飞蝗实际分布拟合度较高。另外, 本研究利用地理信息系统软件ArcGIS对MaxEnt输出的栅格文件进行后期处理, 使目标物种的分布数据与环境变量数据在栅格单元上相对应, 有效的降低了系统误差, 进一步提高了数据的准确性。

本研究预测结果显示, 西藏飞蝗在中国的适生范围为88.79°—113.57°E, 21.52°—38.33°N, 总面积达170.55×10⁴ km²。其中高适生区主要位于四川和西藏两省区, 中适生区则主要集中青藏高原东部地区, 包括四川西部、西藏东部、云南北部等地。上述区域地理上包含了四川的多数草原(若尔盖草原、红原草原、塔公草原、龙灯草原和毛垭草原等)以及西藏雅鲁藏布江及雅砻江河谷。李庆等^[3]研究发现, 西藏飞蝗在草地环境下产卵块数最多、草地环境更适合其生活, 西藏飞蝗可对草地造成严重危害。因此西藏飞蝗在上述地区的危险系数极高, 有关部门应引起高度重视, 严防西藏飞蝗聚集成灾。封传红等^[52]建立了有效积温与地理位置的模型, 获得了西藏飞蝗的潜在分布面积, 结果表明, 1998年西藏飞蝗分布面积最大, 为14.3×10⁴ km²。本文中高适生区面积为22.51×10⁴ km², 中适生区面积则达63.04×10⁴ km², 与封传红的研究结果相差较大, 这可能是因为使用的模型及变量不同所导致。

物种-环境关系是研究物种生态学需求和空间分布的重要方面, 本研究分析了西藏飞蝗存在概率与主导环境变量之间的关系, 获得了相关的反馈曲线。结果表明, 西藏飞蝗的存在概率随着主导环境变量(海拔、等温性、1月份及8月份平均雨量、1月份2月份及12月份平均温度)的变化而改变。7个主导变量中, 海拔和等温性是影响该虫分布最为重要的环境变量, 表明西藏飞蝗扩散和繁殖过程中, 一方面受当地海拔范围的制约, 另一方面也受当地等温性的影响。本研究中西藏飞蝗适宜的海拔范围为2870—5188m, 最适为3941 m。目前已知西藏飞蝗的分布范围仅限于青藏高原地区, 王思忠等^[53]虽曾在四川雅安(海拔560 m左右)室内成功饲养繁殖, 但死亡率高, 且众多研究证明至今未发现该虫在低海拔地区分布, 这与本实验的结果相吻合。等温性反映了温度变化的迟早及幅度, 如当年气温转变迟缓, 则利于西藏飞蝗的胚胎发育, 西藏飞蝗存在概率较高的等温性指标在42.9左右, 这表明该虫在温差较大情况下存在概率较高, 这与上述的研究结果相近。李庆等研究了西藏飞蝗各发育阶段的耐寒性, 结果表明卵期的过冷却点为-22.02℃, 4龄蝗蝻的过冷却点为-6.46℃, 西藏飞蝗以卵越冬。西藏飞蝗蝗卵越冬期长达200 d, 冬季低温会影响其胚胎发育及孵化率。封传红等^[54]研究发现, 温度为限制西藏飞蝗分布的主要限制因子, 夏季高温不是限制因素, 冬季低温可影响西藏飞蝗的安全越冬, 这些都说明冬季低温对蝗卵的越冬至关重要, 本研究认为西藏飞蝗的存在概率与越冬期12月及次年1、2月平均温度有关, 在较低的温度条件下, 该虫仍有较高的存在概率, 但过低的温度将导致其存在概率下降；8月份平均雨量也是影响西藏飞蝗分布的重要因素, 最适值为98.2 mm。研究发现, 西藏飞蝗与飞蝗其他亚种类似, 主要栖居在较低湿的禾本科杂草地带, 该虫一般在8月开始交配、产卵, 适宜的湿度有利于该虫顺利的完成这些生命活动。此处需要注意, 物种存在概率对环境变量的反馈曲线反映的是单一环境变量的作用, 但昆虫的生命活动受到各种环境变量(包括气候因子、寄主条件、天敌种类、植被条件等)的综合影响, 因此此结果可作为判断西藏飞蝗与环境变量之间关系的参考, 但不能完全准确解释二者之间的关系。

尽管MaxEnt模型预测具有操作简便、样本需求量小和预测精度高等优势, 但同其他生态位预测模型类似不可避免的存在一定的局限性。(1)用MaxEnt模型进行预测时所涉及的环境变量除海拔外, 其余变量均为气候变量, 19个生物气候变量为西藏飞蝗实际分布地的气候极值, 即最大值和最小值, 因而, MaxEnt模型所显示的是物种分布的最大可能性, 无法准备表达物种在现实分布的主要地区。上述预测结果更侧重于理解和展现西藏飞蝗潜在地理分布以及揭示适于该物种分布的气候特征。(2)基础生态位是一种理想生态位, 指某一物种在理想的生存条件而不存在其他物种竞争时所占据的最大生态位, 该理论只需考虑非生物因素的影响^[55]。MaxEnt模型对适生区的预测就是基于物种对基础生态位的需求, 但在实际的生存环境中, 生物因子(如物种之间的相互作用、植被类型、地貌特征、物种自身扩散能力和土壤类型等)同样会对预测物种的潜在分布产生重要影响^[56]。例如, 李庆等^[3]研究了植被及土壤对西藏飞蝗产卵及分布的影响, 结果表明, 植被类型、土壤含盐量、含水量和PH值能显著影响西藏飞蝗的产卵及分布。基于上述原因, 可以推测, MaxEnt模型预测的生态位比西藏飞蝗所占据的实际生态位要宽。对此, 在下一步工作中, 除考虑气候因素影响外, 还应注重考量物种间相互作用等生物因子的可信表达, 以改善模型的预测效果。(3)过去30年来, 随着温室气体的不断排放和保持, 全球气候正异常变暖, 而气候变化会引起物种生长及分布格局的变化^[57-59]。本研究中所使用的环境变量来源于世界气候数据库worldclim, 时间跨度为1950—2000s, 缺失了近期10多年的气候数据。因此, 在以后的工作中, 应充实缺失的数据, 使预测结果更为准备可靠。