基于国际贸易和环境气候的生物入侵风险评估综合模型研究初探

段梦格,冼晓青,王书平,钱 程,乔慧捷,李新江

(1.河北大学生命科学学院,保定 071002;2.中国科学院动物研究所,北京 100101;3.中国农业科学院植物保护研究所,北京 100193;4.上海海关动植物与食品检验检疫技术中心,上海 200135)

生物入侵是指生物离开原产地,通过自然或人为因素,到达新的生态环境地区,并在当地形成稳定种群的过程[1-2]。由于在新的环境下缺乏天敌、竞争者等因素,外来生物往往会对入侵地本地物种和生态系统造成破坏,对本地生物多样性、农林牧渔业以及人类健康造成重大经济损失和生态灾难[1,3]。而生态系统健康与人类安全息息相关[4],由入侵物种造成不可逆的生态学效应[5]所带来的生态安全问题会影响人类的健康和社会的可持续发展。然而,生物入侵过程复杂,受到自然环境、地理隔离和人类活动等多种因素的综合影响,这都为量化评估其入侵风险带来阻碍。目前的研究工作大都集中在入侵的某一个阶段,如探讨外来物种随人类活动(贸易、压舱水等)入侵的风险评估、物种在入侵地的潜在分布地预测等,而没有一个方法可以囊括上述入侵过程的各个阶段。因此系统地对入侵物种从原产地-扩散路径-入侵地全链条分析并综合量化评估其入侵危害迫在眉睫。

生物入侵的传播途径主要有3种:有意引入(生物防治、环保性植物、观赏性植物、饲料和宠物等)、无意引入(被感染的货物、包装材料或交通载体)、自然传入(通过风力、水力、鸟类粪便等,利用人类建设的基础设施,扩散到自然状态下不能扩散到的地方)[6-7]。人类的全球化活动,特别是国际贸易被认为是造成有害生物无意引入的主要原因[8-9],是生物入侵的快速途径之一[10]。Maino等在对铃木氏果蝇Drosophilasuzukii的分布预测中,认为人为活动因子对预测结果具有较大的影响[11];人类贸易进入全球贸易阶段以来[12],旅行和贸易以不可思议的速度加快了世界范围内的基因交流[13]。国际贸易不仅可通过直接携带外来生物造成生物入侵,也能通过栖息地的破坏、城市化、大气污染和气候变化等间接影响生物入侵[14-15]。因此,国际贸易是生物入侵的直接与间接驱动力,了解国际贸易在生物入侵中发挥的作用对于保障生物安全和社会经济可持续发展至关重要[14]。

全球入侵物种计划(Global Invasive Species Programme,GISP)的研究表明,预防生物入侵比控制其暴发更重要[9,16]。评估生态风险是预防和控制外来物种入侵的重要方法之一。目前,生态学对外来生物的风险评估方法主要有入侵可能性风险评估、定殖可能性评估,潜在损失风险评估[17-18]。主要根据环境相似性和进化历史两方面进行方法研究。气候环境是影响生物传入后能否定殖的重要因素。气候环境越相似,则两个地区可提供的生态位越相似,从一个地区到另外一个地区的外来生物则更容易定殖,存在较大的入侵风险[19-21]。

入侵地具有满足入侵种生存的空生态位,入侵地和来源地具有相似的气候环境能够满足入侵种的生存,是外来物种能在当地形成稳定种群导致成功入侵的原因之一。如能定量描述物种对未知区域的适应性,即可分析评估物种的入侵风险。生态位模型(ecological niche modeling,ENM)是目前预测物种是否在某个地区具有潜在适生区的重要手段[22],大量的研究工作证明了利用生态位模型对入侵生物进行风险评估的可能。GARP、MaxEnt、CLIMEX、BIOCLIM、ENFA等模型已被用于预测红棕象甲Rhynchophorusferrugineus[23]、异色瓢虫Harmoniaaxyridis[24]、栎列队蛾Thaumetopoeaprocessionea[25]、华山松大小蠹Dendroctonusarmandi[26]等病虫害当前及未来的潜在地理分布。

生物通过某一途径传入入侵地是导致生物入侵成功的另一重要原因[27]。在经济全球化趋势加剧的背景下,跨大洲的国际贸易和旅游以及境内的长距离运输等打破了地理隔离与生态屏障[10],是将物种远距离引入新地区的主要原因。近些年来,随着国际贸易量的增加,因生物入侵引起的恶性事件一直在增加[28]。国际贸易与外来生物入侵的关系一直备受关注。例如,Costello等通过贸易-引入-发现过程的结构模型探究了进口商品的入侵风险与贸易伙伴的关系[21];Hlasny等研究了美国地区的外来昆虫与进口和跨国旅行之间的相关性[29];Westphal等使用回归树在全球尺度下分析了国际贸易对生物入侵的影响[5];Essl等利用外来生物传入记录证实了经济活动是生物入侵的影响因素之一[15];Lantschner等的研究表明全球化是小蠹亚科Scolytinae物种入侵的主要驱动因素[30]等。

由此可见,生物入侵是一个集合了物种自身的环境适应性,扩散能力以及人类活动等多种生物、非生物因素影响的复杂过程,对物种入侵风险的评估需要考虑上述因素综合分析。本研究以重要外来入侵昆虫苹果蠹蛾Cydiapomonella和番茄潜叶蛾Tutaabsoluta(又名Phthorimaeaabsoluta)为例,使用物种的地理分布、目标物种分布地的贸易、入侵地交通运输、大气环境等多种数据,综合分析因物种的自然扩散和人类活动(贸易、交通、旅行等)所带来的入侵过程,并分别对物种的引入风险和定殖可能性进行量化评估。利用物种原产地的贸易数据和交通运输数据计算物种进入口岸后通过运输扩散到其他地区的引入风险(简称TDR);利用贸易数据和地理距离计算物种进入口岸后通过自然扩散到其他地区的引入风险(简称NDR)。将物种的潜在分布作为定殖可能性与前两者结合综合评估国际贸易中的外来生物入侵风险。在本文中,外来物种在国际贸易中到达口岸后通过交通运输导致在中国区域的入侵风险称为TIR,通过自然扩散造成在中国全境区域的入侵风险称为NIR。

1 材料与方法

1.1 数据源

1.1.1贸易数据与检疫数据

从联合国粮食与农业组织(FAO,https:∥www.fao.org/)下载2003年-2015年关于中国进出口相关的国际贸易数据,共22 894条数据,338个商品类型,198个贸易伙伴国。由于一个国家名字可能有多种表达方式,所以需要进行统一合并。根据伙伴国、进出口合并贸易量数据,共有375条数据。其中,中国的进口数据为178条,即与中国进口贸易相关的共有178个国家(图1)。由于本研究以水果蔬菜害虫为例,因此将商品分类,选择水果、蔬菜两个类别的商品。使用标准分数(standard score or z-score)的计算方法将国际贸易进口中的两个类别数据整体进行标准化,以和后续数据保持一致。标准分数的计算方法为:

(1)

图1 2003年-2015年各国对中国的出口贸易量Fig.1 China’s import trade volume from various countries from 2003 to 2015

由于我们只有总的水果和蔬菜的贸易量而没有每个关区具体的进口数据,为了将贸易数据合理地分配到各个海关,我们使用海关的检疫性有害昆虫的检出数据作为依据,即将各海关的昆虫检出比例作为贸易量占比。该数据由合作单位上海海关动植物与食品检验检疫技术中心提供。为了和贸易数据保持一致,我们对检疫数据中的商品进行分类并将各个口岸对应到相关海关,从而整理出各个海关的水果或蔬菜类别中昆虫检出量(分国家或地区)。此外,为了和最新的进境水果指定海关保持一致,根据中华人民共和国海关总署2022年公布的进境水果指定监管场地名单选择水果和蔬菜类别的检疫数据(https:∥github.com/duanmengge/risk-assessment-for-

invasion)。

1.1.2有害生物的地理分布

选择有害生物并根据中国国家有害生物检疫信息平台(http:∥www.pestchina.com/)与欧洲和地中海植物保护组织(European and Mediterranean Plant Protection Organization)建立的EPPO Global Database (https:∥gd.eppo.int)上公布的相关有害生物的地理分布,以国家为单位,选择中国在这些地区的进口水果和蔬菜类贸易量数据用于后续分析。将各地区的进口水果和蔬菜类贸易数据根据对应的昆虫检出数据分配到各个关区(例如,美国水果进境的贸易数据对应美国水果进境的检疫数据)。苹果蠹蛾分布在78个国家或地区,其中具有水果或蔬菜类贸易数据及检疫数据的地区有33个。番茄潜叶蛾分布在99个国家或地区,其中具有水果或蔬菜类贸易数据及检疫数据的地区有35个。苹果蠹蛾与番茄潜叶蛾地理分布地区对中国的出口贸易量及标准化贸易量具体数据见https:∥github.com/duanmengge/risk-assessment-for-invasion。

1.1.3地图和距离

我们搜集了全国42个关区的经纬度。为了量化评估外来物种进入关区后在国内其他地区的扩散概率,从国家基础地理信息中心(http:∥www.ngcc.cn/ngcc/)获取公路与铁路的路网数据,公路路网所用的道路等级为国家高速公路并将其简化处理。利用R包stplanr[31]计算出某个关区到目的地的所有可能的铁路距离和公路距离。通过R包geosphere里的distGeo函数计算得出地理(大圆)距离[32]。交通路网数据用于交通扩散风险分析,地理距离用于自然扩散风险分析。

1.1.4物种发生数据

物种分布数据来源于全球生物多样性信息网络(GBIF,https:∥www.gbif.org/)和相关文献[33-35]。物种发生数据是物种分布点的详细信息,以经纬度为单位。为防止采样偏差经过以下处理:1)去除没有经度或者纬度的数据;2)去除所有重复数据;3)为了匹配建模环境变量的空间分辨率,避免空间自相关性,每个栅格(10′×10′)只保留1条数据。经过处理后,苹果蠹蛾的物种发生数据共有2 177个样本点;番茄潜叶蛾的物种发生数据共有325个样本点(图2)。

图2 苹果蠹蛾(红色)和番茄潜叶蛾(绿色)发生数据Fig.2 Occurrence points of Cydia pomonella (in red) and Tuta absoluta (in green)

1.1.5环境数据

气候环境数据下载于WorldClim数据库(https:∥www.worldclim.org/),共有19个环境气候因子 (bio1～bio19),为了避免多重共线性,结合物种的生物学特性对19个因子进行了Pearson相关性分析,选择环境气候因子用于后续分析。苹果蠹蛾选择bio1(年均温),bio5(最暖月最高温),bio12(年降水量),bio15(降水量的季节性变化系数);番茄潜叶蛾选择bio1(年均温),bio4(气温的季节性变动系数),bio12(年降水量),bio17(最干季降水量)。

1.2 物种的引入风险评估

将筛选出的各个国家对中国水果或蔬菜类贸易量标准化数据依据有害昆虫检出比例分配给各个关区。使用反距离权重法(inverse distance weight,IDW)[36]根据距离对贸易数据进行插值分析。反距离权重法计算公式为:

Dw=e(-β×dg)

(2)

Dp=∑e(-β×dt)

(3)

其中dg为实际各点到各个关区的地理距离,dt为实际各点到各个关区的交通距离,Dw为反距离加权后各点到各个关区的地理距离,Dp为反距离加权后各点到各个关区的交通距离。公式2为计算地理距离的公式。公式3为计算交通运输距离的公式(从一个关区到一个地点不止有一条路,将其反距离加权之后加和在一起为某一地点到某一关区的交通距离),定义β=0.01。

反距离权重法是对采样点进行线性加权来决定输出的栅格值,加权与距离呈负相关,输出点离栅格越远,对输出栅格影响越小。在公路距离内,对其进行反距离加权,并将某个关区到目的地的所有反距离加权路径进行加和,表示某个关区到目的地的公路距离(公式3)。铁路距离内也做了相同的做法。交通距离中公路距离和铁路距离根据权重(Ra,Rb) 进行加和并归一化。本文定义两个的权重相同(1∶1)。公式为:

T=Ra×Da+Rb×Db

(4)

(5)

其中权重Ra=Rb=0.5,Da为通过铁路交通某一点到某一关区的距离,Db为通过公路交通某一点到某一关区的距离。公式5对其合并进行归一化,Dt为归一化后的总体交通数据。

此外,对地理距离同样进行反距离加权计算。

将苹果蠹蛾和番茄潜叶蛾选择的国家的各类型标准化贸易量分配给各个关区后,结合距离数据(交通距离或地理距离),根据如下公式计算得出中国区域内各点的引入风险(公式6)。

(6)

其中,Dij为第i个地点到第j个关区的交通距离(Dtij)或地理距离(Dgij);Cj为第j个关区的标准化贸易量;Ri为地点i物种的TDR(Ri_t)或NDR(Ri_g);n为关区。

1.3 物种的定殖可能性评估

本研究使用生态位模型评估外来生物入侵的定殖可能性。在案例中,我们使用MaxEnt模型对物种的中国区的潜在分布进行预测。在构建MaxEnt模型之前,对其中的特征组合(feature combination,FC)和正则化乘数(regularization multiplier,RM)进行调优。FC主要使用线性特征(L),二次型特征(Q),乘积型特征(P)的不同排列组合(L,Q,P,LQ,LP,QP,LQP),RM选择的区间为1～3,步长为0.5。ENMeval包中的ENMevaluate函数对上述35个组合分别对两个物种进行调优,背景点数设置为10 000,选用70%的训练数据,30%的测试数据。最终采用最小信息准则AICc(the minimum information criterion AICc value,ΔAICc),找到最优参数。优化完成后,采用优化参数使用dismo包中maxent函数进行建模[37-38]。根据kuenm包中的kuenm_proc函数计算最优参数的partial AUC值[39]。并利用ENMeval包中的evaluate函数的计算指标对预测结果进行评估[40]。将建立的模型重新投影到中国区域(10′×10′)进行预测。

1.4 国际贸易中所带来的物种的入侵风险

将物种的引入风险和定殖风险结合在一起评估某个物种由国际贸易带来的入侵风险。公式为:

R=Ri×Rp

(7)

其中,Ri为某个物种的TDR(Ri_t)或NDR(Ri_g),Rp为某个物种的定殖风险,R为物种在国际贸易中的TIR(R_t)或NIR(R_g)。

2 结果与分析

2.1 交通距离和地理距离

铁路和公路距离根据反距离权重、合并、归一化等处理后,共得到42个数据集,分别是中国区域内距离各个关区的交通距离。地理距离通过反距离权重、归一化等处理后同样得到42个数据集,分别是中国区域内距离各个关区的地理距离。

2.2 苹果蠹蛾由国际贸易带来的入侵风险

将筛选后的33个苹果蠹蛾分布地区按照不同的类别(水果、蔬菜)分别根据各关区的昆虫检出量将标准化贸易量分配到各个关区(图3a)。关区标准化后贸易量数据结合交通运输数据和地理距离数据,得出中国区域内各点的引入风险(图3c，3d)。苹果蠹蛾的物种发生数据经过处理之后共有2 177个样本点。使用的4个气候因子都是较为重要的(表1)，且各因子之间的相关性并不高(|r|<0.75)，避免了拟合共线性(图4)。在35组模型调优后选择FC=LQP,RM=1(ΔAICc=0)用于MaxEnt模型,对模型的评估中partial AUC值为1.69,Kappa值为0.54,表现较好,利用此模型对苹果蠹蛾在中国的适生区进行了预测。在MaxEnt模型中4个气候因子的贡献率分别为66.6%(bio1)、30.2%(bio15)、2.6%(bio5)、0.7%(bio12)(表1)。MaxEnt模型模拟了苹果蠹蛾在当前条件下的潜在分布区(图3b),其在全国大部分地区都有一定的定殖可能性,其中云南、四川西部、西藏部分地区、贵州、湖北、陕西、甘肃、新疆等地区的潜在分布可能性较高(Rp>0.3)。当Ri>0.05时,国际贸易货品中所携带的苹果蠹蛾通过交通运输到全国的引入风险(TDR)主要集中在河南、湖北、安徽、湖南、江西等中部地区(图3c),通过自然扩散到全国的引入风险(NDR)在东南部地区更高(图3d)。与TDR相比,NDR在西部地区风险值较高。当R>0.001时,综合考虑交通和自然扩散,苹果蠹蛾在全国大部分地区都有一定的入侵风险(图3e，3f)。风险相对较大的区域主要在中南部地区:云南北部、四川中部、甘肃南部、宁夏部分地区、陕西、山西部分地区、河南西部、重庆部分地区、湖北西部、安徽部分地区、贵州西部(图3e，3f)。东南部地区在国际贸易中的自然入侵风险稍高于交通入侵风险。

图3 以苹果蠹蛾为例的国际贸易所带来的外来入侵风险Fig.3 Alien invasion risk of Cydia pomonella in international trade

图4 苹果蠹蛾入侵风险模型中的4个气候因子的相关性分析Fig.4 Correlation analysis of four bioclimatic variables for ecological niche modeling of invasion risk of Cydia pomonella

表1 气候因子在苹果蠹蛾入侵风险模型构建中的贡献率Table 1 Contribution rate of bioclimatic variables in ecological niche modeling of invasion risk of Cydia pomonella

2.3 番茄潜叶蛾由国际贸易带来的外来入侵风险

将筛选后的番茄潜叶蛾的35个分布地区按照不同的类别(水果、蔬菜)分别根据各关区的昆虫检出量将标准化贸易量分配到各个关区(图5a)。关区标准化后贸易量数据结合交通运输数据和地理距离数据,得出中国区域内各点的引入风险(图5c，5d)。番茄潜叶蛾的物种发生数据经过处理之后共有325个样本点用于进一步分析。使用的4个气候因子都是较为重要的(表2),且各因子之间的相关性并不高(|r|<0.75)，避免了拟合共线性(图6)。在35组模型调优后选择FC=LQ,RM=1(ΔAICc=0)用于MaxEnt模型,对模型的评估中partial AUC值为1.34,Kappa值为0.84,表现较好,利用此参数建模对番茄潜叶蛾在中国的潜在分布区进行了预测。MaxEnt模型中贡献率分别为49.3%(bio1)、44.3%(bio4)、5.5%(bio12)、0.9%(bio17)(表2)。MaxEnt模型模拟了在当前条件下番茄潜叶蛾的潜在分布区(图5c)集中在中国中南部地区,其中云南、四川东部、贵州、重庆、广西、湖南、西藏东南部、甘肃南部、陕西南部、湖北、江苏、浙江、福建等地的潜在分布区面积较大且可能性较高(Rp>0.5)；此外,中南部其他地区如河南、安徽、江西等地部分地区及山东部分沿海地区的潜在可能性同样较高(Rp>0.5)。国际贸易货品中所带来的番茄潜叶蛾通过交通运输到达全国的引入风险(TDR)主要集中在中部地区(图5c),通过自然扩散到达全国的引入风险(TIR)则主要集中在东南部地区(图5d)。综合番茄潜叶蛾由国际贸易带来的引入风险和物种的潜在分布的定殖可能性,将两者归一化之后评估入侵风险(图5e,5f)。在R>0.001时,将TIR与NIR相结合,除了新疆、黑龙江、吉林等地外,中国区内其他各省市均有一定入侵风险(图5e,5f)。番茄潜叶蛾的可能入侵地主要集中在中国东南部地区,福建部分地区的交通入侵风险高于自然入侵风险。

图5 以番茄潜叶蛾为例的国际贸易带来的外来入侵风险Fig.5 Alien invasion risk of the Tuta absoluta in international trade

表2 气候因子在番茄潜叶蛾入侵风险模型构建中的贡献率Table 2 Contribution rate of bioclimatic variables in ecological niche modeling of invasion risk of Tuta absoluta

3 讨论

3.1 外来生物入侵评估方法

从生态学的角度对贸易进行分析,对评估贸易带来的外来入侵风险、制定关于生物入侵风险方面的政策起着关键作用[41]。国际贸易相关的货品经由海关进入并通过物流运往其他地区,入侵生物也随着这个过程扩散到适宜其分布的地区,因此入侵生物在我国境内的入侵风险需要结合环境适宜度、自然扩散与交通运输等多种因素综合评估。本文提出了一种将交通运输距离与自然扩散的空间距离分别对贸易量进行插值权重分析,并以此来评估外来物种引入风险的新方法。本研究在原有的物种潜在分布地预测的基础上,结合了贸易带来的引入风险分析,详细地评估国际贸易带来的某一特定有害生物的入侵风险。本文利用该研究方法,以番茄潜叶蛾和苹果蠹蛾为例,对其入侵风险进行了量化评估。

3.2 苹果蠹蛾由国际贸易带来的入侵风险分析

苹果蠹蛾地理分布的大部分地区均与中国有贸易交流,通过贸易引入的风险较大。且全国绝大部分地区为苹果蠹蛾的适生区。因此一旦引入,就会快速扩张,造成生物入侵。苹果蠹蛾是果树的重要害虫之一,一旦入侵会对国内的农产品种植业造成重大影响。综合TIR与NIR的评估结果,甘肃和陕西是入侵风险相对较高且适生区面积较大的省份,辽宁、河北、山东等省份预测为中高风险区。除了青藏高原、内蒙古和黑龙江北部、广西、广东、海南等部分地区的入侵风险相对较低,其他各地区均有一定的风险。虽然西北地区关区较少,引入风险相对较低,但模型预测表明在西北地区(新疆、甘肃、内蒙古,青海部分地区)有较高的定殖可能性。虽然在国际贸易中所带来的入侵风险不如中部地区,但是一旦入侵成功就会造成不可估量的经济损失。

3.3 番茄潜叶蛾由国际贸易带来的入侵风险分析

相较于苹果蠹蛾,番茄潜叶蛾在世界范围内分布更广(99个国家或地区),通过贸易引入的可能性较大。综合TIR和NIR,番茄潜叶蛾在西南地区的入侵风险较高,其次是中东部地区。云南是圣女果种植总面积和产量最大的地区[33,42],又是番茄潜叶蛾的最适潜在分布区,具有很高的入侵风险;四川盆地及周边地区的入侵风险也较高;作为我国番茄种植面积和产量最高[33],且为番茄潜叶蛾首发地的新疆地区,虽然引入风险不高,但具有一定的定殖风险(Rp>0.1)。同时,考虑到番茄潜叶蛾入侵过程中生态位可能会发生变化,且温室大棚创造的微环境可能为番茄潜叶蛾提供越冬场所等因素,因此不能对新疆地区的番茄潜叶蛾防控放松警惕。

3.4 方法的有效性与局限性

本文通盘考虑了物种的潜在分布、国际贸易、交通运输和人类活动等多种因素,综合各个因素在入侵过程中的权重来评估国际贸易带来的生物入侵风险,并利用两个案例对由于上述因素所导致的入侵风险进行了量化评估,初步量化分析了由国际贸易带来的外来生物入侵风险。人类活动在国际贸易中有着重要作用,是造成国际贸易带来的外来生物入侵难以量化评估的主要原因。本研究尝试用交通运输等数据代替人类活动因素,结合潜在地理分布,量化评估由国际贸易带来的有害生物入侵风险。但入侵过程中的影响因素不止于此,特别是农业生产的特殊性(如温室大棚等创造的人工微环境,跨区域的运输导致的物种扩散不确定性等)导致的入侵物种出现在模型预测区域以外,这些不可控的因素都需要在风险评估过程中加以考虑。此外,由于数据可获得性的问题,本文并未使用对应货品真实的运输量数据,而是使用了道路距离作为替代指标,导致模型性能下降,这些都需要在今后的研究中注意和加以改进。

4 结论

本文提出了一种综合多源数据对外来物种入侵过程进行全链条分析的方法,并以苹果蠹蛾和番茄潜叶蛾为例,对其入侵风险进行了量化评估,讨论了该方法的应用范围和局限性。利用本方法可以综合评判物种扩散过程中各个环节对入侵的影响程度,寻找过程中的关键节点,为精准防控国际贸易中所带来的外来入侵提供了新思路。