入侵植物是生物多样性丧失的主要原因之一，入侵植物相较于本土植物在捕获阳光、争夺土壤水分及养分等方面优势显著，具有更强的竞争能力^[1-2]。如今全球几乎所有区域都受到入侵植物的影响，入侵植物给各地生态系统的平衡与稳定带来极大冲击^[3]。中国是一个生态系统高度多样化的国家，是植物入侵较严重的国家之一^[4]。入侵杂草众多，其中以菊科最多，来源地多数为美洲，肿柄菊Tithonia diversifolia A. Gray和金腰箭Synedrella nodiflora (L.) Gaertn.位列其中^[5-6]。

肿柄菊和金腰箭均为菊科入侵植物，原产于中南美洲，现广泛分布于我国云南、华南等地 ^[7-8]。肿柄菊入侵后形成密实的大面积单优势种群落^[7]，其植株构件间的生长关系随着群落环境变化而改变，有利于其入侵不同环境，增强其适应环境的能力^[9]。金腰箭入侵后通过高种子产量和快速营养生长扩张种群^[10]。二者常见于农田、路边、林地及河岸等生境，通过生态位竞争严重危害农作物生产，并威胁区域生物多样性^{[4, 11]}。

鉴于上述肿柄菊和金腰箭的入侵特性及扩散威胁，准确预测其潜在适生区成为生态防控的关键。最大熵模型(Maximum entropy models，MaxEnt)是一种物种分布预测的常用软件，只需要物种分布数据和环境数据，就可以统计分析出物种的当前分布和未来分布，且只需要少量样本就可以达到较高的精确度^[12]，该模型被广泛应用于植物的各种分布区域的模拟，是外来物种潜在适生区预测的首选模型^[13]。基于上述模型特性，本研究通过MaxEnt模型和ArcGIS 软件等对影响入侵植物肿柄菊和金腰箭分布的主导气候因子进行分析，使用ROC曲线法评价模型预测的准确性，并对其当前和未来不同气候条件下(至2040年)在中国的潜在适生分布进行研究，以期为肿柄菊和金腰箭这2种杂草在中国的生态治理提供可靠的数据支撑。

1.   材料与方法

1.1   物种分布预测方法

物种分布预测方法参考当前主流方法^[14-15]进行，主要包括数据的查找与处理、模型的优化及构建、可视化地图的制作等。首先进行物种分布点收集及处理；继而查找对物种分布产生影响作用的各类环境因子，并运用SPSS软件进行相关性分析，同时使用ENMeval对模型进行优化，随后构建MaxEnt模型，最终借助ArcGIS软件绘制可视化图。

1.2   采用的软件和程序

MaxEnt软件3.4.1版、ArcGIS软件10.4.1版、SPSS Statistics软件20.0版、R语言3.6.3版、ENMeval程序包以及Microsoft Excel 2016。

1.3   环境变量因子获取

从世界气候数据网站(https://www.worldclim.org/)下载2.1版、分辨率为5 min(10 km²)的1970—2000年的19个生物气候变量(bio1~bio19)以及高程海拔(elev)、坡度(slope)和坡向(aspect)数据，并于此网站下载国家气候中心系统模式(BCC-CSM2-MR)的4个气候情景模式 (SSP1-2.6、SSP2-4.5、SSP3-7.0和SSP5-8.5)下2021—2040年的世界范围内19个生物气候变量数据(bio1~bio19)，从联合国粮农组织网站的世界土壤数据库下载17个土壤属性变量。一共得到温度、降水量、海拔和土壤特征的环境变量数据39个(表1)。

表  1  环境变量因子

Table  1.  Environmental variable factors

变量 Variable	描述 Description		变量 Variable	描述 Description
bio01	年平均温度		elev	高程海拔
bio02	昼夜温差月均值		slope	坡度
bio03	等温性		aspect	坡向
bio04	温度季节变化标准差		awc_class	土壤有效含水量
bio05	最暖月最高温		t_bs	表层土壤容重
bio06	最冷月最低温		t_CaCO3	表层碳酸钙含量
bio07	年均温度变化范围		t_CaSO4	表层硫酸钙含量
bio08	最湿季度平均温度		t_cec_soil	表层土壤阳离子交换能力
bio09	最干季度平均温度		t_clay	表层黏土含量
bio10	最暖季度平均温度		t_ece	表层盐分含量
bio11	最冷季度平均温度		t_esp	表层碱分含量
bio12	年均降水量		t_gravel	表层土壤砾石含量
bio13	最湿月降水量		t_oc	表层有机碳含量
bio14	最干月降水量		t_pH_H2O	表层酸碱度
bio15	季节降水量变异系数		t_sand	表层底层沙含量
bio16	最湿季度降水量		t_silt	表层粉沙粒含量
bio17	最干季度降水量		t_teb	表层交换性盐基
bio18	最暖季度降水量		t_usda_tex_class	表层Usda土壤质地分类
bio19	最冷季度降水量

下载: 导出CSV 

| 显示表格

使用ArcGIS读取过拟合分析后的物种分布数据及环境数据，将提取的环境变量和物种分布数据导入MaxEnt模型运行10次，得到环境数据的平均贡献率，并对提取的环境变量通过SPSS进行Pearson相关性分析。计算变量间的相关性矩阵，剔除平均贡献率为0的环境变量，并将相关系数R>0.8的2个环境变量中贡献率较低的1组舍去。最终，得到与肿柄菊相关的环境变量因子25个：bio02、bio04、bio07、bio10、bio11、bio12、bio14、bio15、bio18、bio19、awc_class、drainage、t_bs、t_CaCO₃、t_CaSO₄、t_cec_soil、t_clay、t_ece、t_esp、t_gravel、t_oc、t_pH_H₂O、t_sand、t_silt、t_teb；得到与金腰箭相关的变量因子26个：bio02、bio04、bio05、bio06、bio07、bio08、bio12、bio13、bio15、bio17、bio18、bio19、awc_class、drainage、t_bs、t_CaCO₃、t_cec_soil、t_clay、t_ece、t_esp、t_gravel、t_oc、t_pH_H₂O、t_sand、t_silt、t_teb。

1.4   物种分布数据拟合分析

通过在中国数字植物标本馆、全球生物多样性信息网络、中国科学院植物研究所植物标本馆等数字标本库检索目标植物学名，得到肿柄菊和金腰箭分布点信息，将相关数据存为CSV格式。将去掉无效数据后的物种分布数据通过ENMeval软件及环境数据进行过拟合分析，对分布数据的经纬度根据栅格大小去掉同一栅格内重复数据，得到肿柄菊自然分布点3859个、金腰箭自然分布点2402个。

1.5   标准地图数据获取

从标准地图服务系统(http://bzdt.ch.mnr.gov.cn/)下载GS(2022)4312号标准地图。

1.6   物种环境变量因子的分布数据模型优化

分别将得到的肿柄菊和金腰箭的环境变量因子使用R语言ENMeval进行分布数据模型优化，将分布数据随机分成2部分，其中，75%用于训练集，25%用于测试集， 对MaxEnt软件中的调控倍频(Regularization multiplier, RM)和特征组合(Feature combination, FC)参数进行优化，RM设置40个调控倍频，从0.1到4.0，每次增加0.1，FC一共有5个特征：Linear(L)、Quadratic(Q)、Hinge(H)、Product(P)、Threshold(T)，设定29个特征组合(H、L、LH、LP、LPH、LPT、LPTH、LQ、LQH、LQP、LQPH、LQPT、LQPTH、LQT、LQTH、LT、P、PH、PT、PTH、Q、QH、QP、QPH、QPT、QT、QTH、T和TH)，R语言程序包对上述1160种参数组合进行测试检验，用最小赤池信息量准则的AICc模型评估ROC曲线下面积(Area under curve，AUC)，用AUC评估模拟结果的精确度，使用AICc值为“0”的参数组合用于MaxEnt模型建模。

1.7   环境变量因子和分布数据的模型运算

分别将肿柄菊和金腰箭相关的环境变量因子和处理后的分布数据输入MaxEnt模型，其中75%的分布数据为训练集，25%作为测试集，RM与FC参数选择模型优化后的参数，选择制作响应曲线，分析单个环境变量因子和分布概率之间的关系，使用刀切法(Jackknife)检验环境变量的重要程度。以AUC来评价模型的精确度，AUC越大，预测结果的精确度越高，模型重复运行10次，取均值作为模型运算结果。

1.8   物种分布图的绘制

分别将肿柄菊和金腰箭的MaxEnt模型优化后的数据和中国地图数据读入ArcGIS软件，进行掩膜提取，得到对应物种在我国的物种分布图，并利用自然间断点对物种分布图进行重分类，得到非适生区、低适生区、中适生区和高适生区4个等级。

2.   结果与分析

2.1   肿柄菊和金腰箭的MaxEnt生态位模型预测结果

基于25个环境变量因子和肿柄菊3859个分布点数据，MaxEnt软件调控倍率RM=1.2，特征组合FC为LH。经MaxEnt 10次重复训练得到的AUC为0.875(图1a)，预测结果可靠。

图  1  肿柄菊(a)和金腰箭(b)生态位模型预测结果

Figure  1.  The prediction results of niche model of Tithonia diversifolia A. Gray (a) and Synedrella nodiflora (L.) Gaertn (b)

下载: 全尺寸图片 幻灯片

基于26个环境变量因子和2402个金腰箭分布点数据，MaxEnt软件调控倍率RM=3.5，特征组合FC为LPTH，经Maxent 10次重复训练得到的AUC为0.903(图1b)，预测结果可靠。

2.2   影响肿柄菊和金腰箭分布的环境变量因子

肿柄菊使用25个环境因子建立模型，贡献率前10的环境变量因子和贡献率分别为温度季节变化标准差(bio04，37.2%)、最冷季节平均温度(bio11，33.6%)、最暖季节降水量(bio18，6.9%)、高程海拔(elev，4.2%)、年均温度变化范围(bio07，4.2%)、坡度(slope，3.8%)、表层交换性盐基(t_teb，2.5%)、最干月降水量(bio14，1.3%)、昼夜温差月均值(bio02，1.3%)和土壤有效含水量(awc_class，0.9%)，前10的环境变量因子贡献率之和为95.9%。贡献率越高表明植物分布受此环境变量影响的概率越高。

金腰箭使用26个相关环境因子建立模型，贡献率前10的环境变量因子和贡献率分别为最湿月降水量(bio13，49.6%)、年均温度变化范围(bio07，16.5%)、年均降水量(bio12，10.6%)、温度季节变化标准差(bio04，3.8%)、表层交换性盐基(t_teb，3%)、表层土壤烁石含量(t_gravel，2.9%)、昼夜温差月均值(bio02，2.6%)、高程海拔(elev，2.6%)、最湿季度平均温度(bio08，2.3%)和季节降水量变异系数(bio15，1.5%)。前10的环境变量因子贡献率之和为95.4%。

刀切法检验环境变量的结果如图2所示。综合MaxEnt模型和Jackknife检验结果可知，温度季节变化标准差(bio04)、最冷季节平均温度(bio11)是影响肿柄菊分布的主要环境变量，2个环境变量因子贡献率之和为70.8%；而年均温度变化范围(bio07)、年均降水量(bio12)、温度季节变化标准差(bio13)是影响金腰箭分布的主要环境变量，3个环境变量因子贡献率之和达76.7%。

图  2  环境变量重要性的刀切法检验结果

Figure  2.  Results of Jackknife method for the importance of the environmental variables

下载: 全尺寸图片 幻灯片

2.3   影响肿柄菊和金腰箭分布的重要环境因子的相关阈值

存在概率≥0.5对应的环境因子范围是物种适宜分布的阈值。根据响应曲线(图3)确定影响肿柄菊和金腰箭分布的贡献率高的环境因子的阈值。温度季节变化标准差(bio04)处于−200%~300.02%(图3a)以及最冷季节平均温度(bio11)处于13.38~22.66 ℃(图3b)时肿柄菊存在概率≥0.5。当年均平均温度变化范围(bio07)处于3.69~18.06 ℃(图3c)、年均降水量(bio12)处于1542.86~4563.11mm(图3d)以及最湿月降雨量(bio13)处于300.05~2439.07mm(图3e)时金腰箭存在概率≥0.5。

图  3  环境因子响应曲线

Figure  3.  Environmental factor response curve

下载: 全尺寸图片 幻灯片

2.4   肿柄菊和金腰箭在当前气候下的潜在适生分布区

将肿柄菊和金腰箭在MaxEnt模型模拟优化后的数据和中国地图数据分别导入ArcGIS，得到肿柄菊和金腰箭在当前气候下的潜在分布区(图4)。肿柄菊当前总适生面积约152.72万km²，约占国土面积的15.91%，高适生区、中适生区和低适生区面积分别为31.83、42.46和78.43万km²，肿柄菊的适生范围(图4a)主要在西藏南部、云南、四川、重庆以及华南地区，其中，高适生区和中适生区零散分布于位于我国西藏南部、云南、广西、广东、福建、海南和台湾。金腰箭当前总适生区面积约250.68万km²，约占国土面积的26.11%，高适生区、中适生区和低适生区面积分别为33.82、56.49和160.37万km²，金腰箭的适生范围(图4b)主要在西藏南部、云南、四川、重庆以及华南地区，其中，高适生区和中适生区零散分布于位于我国西藏南部、云南、广西、广东、福建、海南和台湾。

图  4  肿柄菊(a)和金腰箭(b)在中国的当前适生区

Figure  4.  Current suitable areas of Tithonia diversifolia A. Gray (a) and Synedrella nodiflora (L.) Gaertn (b) in China

下载: 全尺寸图片 幻灯片

2.5   肿柄菊和金腰箭在未来气候模型下的潜在适生分布区

表2为在不同时期气候情景下肿柄菊和金腰箭在中国的潜在适生区面积。

表  2  不同时期气候情景下肿柄菊和金腰箭在中国的潜在适生区面积

Table  2.  Potential suitable areas of Tithonia diversifolia A. Gray and Synedrella nodiflora (L.) Gaertn in China under different climate scenarios in different periods ×10⁴ km²

物种 Species	时期 Tperiod	气候情景 Climate scenario	低适生区 Low suitable areas	中适生区 Middle suitable areas	高适生区 High suitable areas	总适生区 Total suitable area
肿柄菊 Tithonia diversifolia A. Gray	当前 Current period	-	78.43	42.46	31.83	152.72
	2040s	SSP1-2.6	145.34	51.93	38.62	235.89
	2040s	SSP2-4.5	147.99	52.50	35.91	236.40
	2040s	SSP3-7.0	142.08	51.60	36.98	230.67
	2040s	SSP5-8.5	142.18	50.33	37.23	229.74
金腰箭 Synedrella nodiflora L. Gaertn.	当前 Current period	-	160.37	56.49	33.82	250.68
	2040s	SSP1-2.6	146.40	123.85	93.78	364.03
	2040s	SSP2-4.5	73.23	35.04	19.84	128.11
	2040s	SSP3-7.0	145.16	54.02	54.99	254.17
	2040s	SSP5-8.5	178.97	39.65	26.29	244.92

下载: 导出CSV 

| 显示表格

SSP1-2.6情景下，2040年肿柄菊总适生区面积约235.89万km²，较当前分布面积增加83.18万km²，主要为低适生区面积增加66.91万km²。金腰箭在此气候情境下2040年总适生区面积约364.03万km²，较当前金腰箭分布面积增加约113.35万km²，主要为高适生区增加59.96万km²和中适生区面积增加67.36万km²。

SSP2-4.5情景下，2040年肿柄菊总适生区面积约236.40万km²，较当前分布面积增加83.68万km²，主要为低适生区面积增加69.56万km²。金腰箭在此气候情境下2040年总适生区面积约128.11万km²，较当前分布面积减少122.57万km²，主要为低适生区面积减少87.14万km²。

SSP3-7.0情景下，2040年肿柄菊总适生区面积约230.67万km²，较当前分布面积增加77.95万km²，主要为低适生区面积增加63.65万km²。金腰箭在此气候情境下2040年总适生区面积约254.17万km²，较当前分布面积相对变化不大。

SSP5-8.5情景下，2040年肿柄菊总适生区面积约229.74万km²，较当前分布面积增加77.02万km²，主要为低适生区面积增加63.75万km²。金腰箭在此气候情境下2040年总适生区面积约244.92万km²，较当前分布面积相对变化不大。

图5为4种气候情景模式下肿柄菊和金腰箭在中国的适生区，肿柄菊在2040年，SSP1-2.6情景下(图5a)、SSP2-4.5情景下(图5b)、SSP3-7.0情景下(图5c)、SSP5-8.5情景下(图5d)下总适生区明显增加，主要为长江以南部分非适生区变为低适生区，广东省、广西壮族自治区和福建省部分中适生区变为高适生区，低适生区变为中适生区。

图  5  4种气候情景模式下肿柄菊(a~d)和金腰箭(e~h)在中国的适生区

Figure  5.  Suitable areas of Tithonia diversifolia A. Gray (a-d) and Synedrella nodiflora (L.) Gaertn (e-h) in China under four climate scenarios

下载: 全尺寸图片 幻灯片

金腰箭在SSP1-2.6情景下2040年(图5e)总适生区明显增加，其中川渝地区部分低适生区变为高适生区，华南部分中适生区变为高适生区，渤海和黄海沿岸区域部分非适生区变为低适生区。金腰箭在SSP2-4.5情景下2040年(图5f)总适生区明显减少，其中主要为华南地区部分低适生区变为非适生区，广西南部及藏南地区高适生区减少为中适生区。金腰箭在SSP3-7.0情景下2040年(图5g)总适生区变化不大，低适生区面积减小，同时川渝地区和华南地区部分低适生区变为中适生区和高适生区。金腰箭在SSP5-8.5情景下2040年(图5h)总适生区变化不大，低适生区面积增大，藏南及广西壮族自治区南部的中高适生区面积减少。

3.   结论与讨论

MaxEnt模型的精度取决于样本的覆盖度、区域和样本量的多少 ^[16]。本研究通过国内外植物网站分别收集肿柄菊、金腰箭大量分布数据，经过处理得到肿柄菊分布点3859个，金腰箭分布点2402个，分布点数据数量充足。AUC 是受试者工作特征曲线(ROC)下面积，AUC是衡量模型结果的指标，AUC的取值范围是0~1，数值越大表明环境变量与预测物种地理分布之间的相关性越大，模型判断力越强，模型预测效果越准确^[17]，本研究结果显示肿柄菊的AUC为0.875，金腰箭的AUC为0.903，模型预测结果均较为可靠。

植物的分布是光照、气温、降水量、土壤及地形等多种因素综合作用的结果，但在一定尺度上，气候因子是制约其适生区分布的主要因素^[18-20]。环境因子的数量、种类和精度对模型预测的准确度至关重要^[21-22]。本研究选取了19个气候因子、3个地形因子和17个土壤因子，共计39个环境因子用于模型预测，预测环境变量范围广泛。环境因子的分辨率为5 min (10 km²)，精度较高。环境变量间存在着较高的空间相关性^[17]，为提高模型的精确性，避免引入的无用信息影响模型预测结果，针对不同植物要对上述环境因子进行筛选。最终肿柄菊使用25个环境因子建立模型，贡献率前10的环境变量因子之和为95.9%。其中，影响肿柄菊分布的主要环境因子为温度季节变化标准差和最冷季节平均温度，贡献率之和达70.8%，对其分布起着关键的决定性作用。当温度季节变化标准差为−200%~300.02%以及最冷季节平均温度为13.38~22.66 ℃时，肿柄菊存在概率≥0.5。采用26个相关环境因子建立金腰箭的模型，贡献率前10的环境变量因子之和为95.4%。其中，影响金腰箭分布的主要环境因子为平均温度变化范围、年均降水量和最湿月降雨量，贡献率之和达76.7%，对金腰箭的分布起着关键的决定性作用。当环境处于平均温度变化范围3.69~18.06 ℃、年均降水量达到1542.86~4563.11 mm以及最湿月降雨量为300.05~2439.07 mm时，金腰箭存在概率≥0.5。这些气候条件分别为肿柄菊和金腰箭的生存提供了适宜的生态位。

当前肿柄菊已在我国南部的热带、亚热带一些省份普遍分布^[23]，本研究结果显示肿柄菊在中国当前的高适生区和中适生区分布于我国西藏南部、云南、广西、广东、福建、海南、台湾，总适生面积约152.72万km²。对于未来肿柄菊的分布，不同气候情景下肿柄菊的潜在适生分布范围相差不大。肿柄菊未来预测分布相较当前入侵分布总适生区明显增加，主要为中国长江以南部分，非适生区变为低适生区，低适生区约以倍数增加，并且其未来中高适生区有向广东、广西和福建区域扩张的趋势，这些地区要加强防范。

现阶段金腰箭的高适生区和中适生区零散分布于我国西藏南部、云南、广西、广东、福建、海南、台湾，总适应区面积为250.68万km²，这与《中国植物志》^[24]记载基本一致。对于未来金腰箭的分布，不同气候情景下金腰箭均有逐渐向我国东南部地区扩张的趋势，这些地区要加强对其监控及防治。