喜旱莲子草气候生态位的特征及在入侵过程中的变化

田德锋，李耕耘，张文驹

(复旦大学 生命科学学院 生物多样性与生态工程教育部重点实验室， 上海 200433)

喜旱莲子草(Alternantheraphiloxeroides(Mart.) Griseb)是最恶性的入侵植物之一，原产于南美洲，现侵入北美洲、亚洲和大洋洲等的许多国家和地区，对农业(如牧场和农田)和自然水生生态系统造成严重威胁[1]。研究者已经从不同的角度对其进行了深入的研究[2-9]，发现喜旱莲子草在原产地以有性繁殖为主，多种昆虫兼食或专食该植物[2]，而在入侵区域，喜旱莲子草以无性繁殖为主，几乎没有天敌限制；前期的研究还发现喜旱莲子草在原产地种群的遗传多样性丰富，在入侵地种群的遗传多样性则严重降低，中国境内的种群甚至可能是同一株系的克隆后代[5]。到目前为止，人们对这一物种的气候生态位，特别是入侵种群的气候生态位的变化却鲜有研究，而这些特征对于预测该物种的入侵并进行管理是极为关键的信息。在全球气候变暖的背景下，探究喜旱莲子草的气候生态位及入侵过程中可能的变化显得尤为重要，特别是需要从全球尺度来解析其入侵种群气候生态位的变化及与入侵地理范围之间的联系。

已经有多种方法利用地理空间(Geographic space, G-space)和环境空间(Environmental space, E-space)模拟和测试物种生态位变化[10-11]，这些方法已经为评估凤眼蓝(Eichhorniacrassipes)、紫荆泽兰(Ageratinaadenophora)等恶性入侵种的入侵格局和生态特征提供了检测工具[12-13]。地理空间模拟是通过提取观测记录(经纬度信息)对应的环境信息来构建物种分布模型(Species Distribution Model， SDM)[14]，然后将该模型投影到地理空间(G-space)中构建潜在分布适生区。与SDM在地理空间中执行相关模型算法不同，E-space中常用的模型由Broennimann等[1]利用主成分分析(Principal Component Analysis, PCA)方法构建，该方法的优势是可以在相同气候空间中描述多个引入种群的现实气候生态位，这有助于比较不同入侵种群与原产地种群的生态特征。环境PCA的另一优势在于它可以对观测记录进行核平滑函数处理，从而降低了低密度潜在采样偏差带来的风险。

本研究在全球范围内广泛收集喜旱莲子草的观测记录，并利用上述两种生态位模型来解析该物种在原产地(南美洲)及各主要入侵区域(包括亚洲、北美洲和大洋洲)的气候生态位。我们关注的问题是喜旱莲子草在不同区域的入侵种群所具有的气候生态位与原产地种群的气候生态位相比，发生了什么样的变化，及这种变化与入侵地理区域之间的联系。这些分析将有助于在宏观尺度上理解喜旱莲子草入侵种群在地理和生态空间变化的特征，从而为该物种的入侵机制和生态管理提供建议。

1 数据与方法

1.1 分布数据的采集和处理

喜旱莲子草的分布数据包含文献记载和开放数据库共享的物种分布数据记录[15]。开放数据库包括： 全球生物多样性信息网络(GBIF, https：∥www.gbif.org/)[16]，澳大利亚生物多样性信息网络(Atlas of Living Australia, https：∥www.ala.org.au/)[17]，物种链接网络(specieslinks, http：∥www.splink.org.br)，美国国家综合数字化生物采集网络(iDigBio, https：∥www.idigbio.org/)，美国纽约植物标本馆(The New York Botanical garden Virtual herbarium, https：∥www.nybg.org/),中国国家标本资源共享平台(NSII, http：∥nsii.org.cn)，中国教学标本资源共享平台(http：∥mnh.scu.edu.cn)和中国生物多样性与生态安全大数据平台(bioone, http：∥www.bio-one.org.cn/)。建模用的分布数据去除了GBIF收录的iNaturalist(https：∥www.inaturalist.org/)网站提供的信息(存在很多鉴定错误)，同时仅保留1970—2015年的喜旱莲子草分布信息以和建模变量相匹配。为降低空间聚类和空间采样偏差可能造成的影响，使用R包Spthin划定4 km范围作为分布点间的稀疏范围。地理背景图层来自自然地球网站(Natural Earth, https：∥www.naturalearthdata.com/)所提供的栅格地球底图。经数据矫正及清理后，各大洲实际参与建模的分布数据为南美洲(340个)、亚洲(590个)、北美洲(1 927个),大洋洲(287个)。如图1(a)所示，这些种群分别缩写为SA(南美原产地种群)、AS(亚洲入侵种群)、NA(北美入侵种群)和OC(大洋洲入侵种群)。欧洲区域种群因分布点过少和明显采样偏差的原因不参与建模统计。

图1 喜旱莲子草在全球范围内的地理分布图及各区域种群的气候生态位重叠示意图Fig.1 The geographical distribution of Alternanthera philoxeroides in the world and the overlap map among the climatic niches in different continents(a)为喜旱莲子草全球分布图，黑色和灰色点表示实际观测分布点。(b)和(c)分别表示喜旱莲子草不同区域种群环境信息的标准主成分分析和基于环境背景信息重建的现实气候生态位。(b)图中不同环境变量的中文释义见方法中“环境变量筛选”部分。(c)图中虚线范围内表示每个区域种群70%的观测密度气候(去除30%边缘气候)，外围实框线表示包含观测信息在内的100%环境背景信息。在图1中，共同使用(c)中的图例，绿色、红色、蓝色、黄色依次标识原产地(南美洲，SA)、北美洲(NA)、亚洲(AS)和大洋洲(OC)，其中(a)中的面表示地理分布区域，(b)中的点表示种群，(c)中的线包络表示生态空间。

1.2 环境变量筛选

研究最初选取27种宏观环境变量来表征喜旱莲子草的气候生态位特征，其中包括19个生物气候变量[18]以及其他8个涉及温度、辐射和水文的宏观变量[19-20]，这些变量常被认为是影响动植物分布的重要因素。参照前人的建模经验筛选环境变量[21-24]，这些经验包括优先选择影响物种分布的近端作用变量，去除存在强共线性关系的变量(|r|>0.8),去除生物气候变量中部分温度和降水结合的变量等。总共筛选到8个环境变量用于表征喜旱莲子草的气候生态位，包括BIO2(日温度平均范围)，BIO10(最暖季节的月平均温度)、BIO11(最冷季节的月平均温度)、BIO12(年总降水量)、BIO16(最潮湿季节的总降水量)、BIO17(最干燥季节的总降水量)、WET(年降水天数)、GHI(总水平辐射)。上述8个气候变量常被用于模拟入侵植物的地理分布，它们能够表征喜旱莲子草的生理和生存所必需的环境条件。

1.3 环境空间(E-space)中比较生态位变化

首先，统计建模环境变量在不同大洲间的差异，使用非参数Kruskal-Wallis检验组间变量差异(FDR矫正)，以及利用Games-Howell测试事后检验组内变量差异。其次，基于原产地和3个入侵地分布记录和随机背景点的环境信息构建物种已实现的气候生态位。环境生态位的评估框架由Broennimann和Warren等人开发[10-11]，利用前述筛选的环境变量来构建两个主轴为代表的二维主成分环境背景空间(PCA-env)。

基于上述构建的环境背景空间，研究进一步测试喜旱莲子草入侵种群与原产地种群(SA)的生态位重叠、生态位保守性和生态位动态指数。进化支之间的生态位重叠采用Schoener的D值[25]量化，D值的度量范围为0(无重叠)到1(完全重叠)。生态位保守性检验包含生态位等同性和相似性测试两部分，前者是评估两个生态位是否可互换，仅考虑实际发生的气候信息，后者则是检验两个分类群或进化分支的环境背景空间是否相似。两种测试都是比较随机分布与实际观测的生态位重叠值(D值)来评估两进化支的生态位是否比随机的显著不等同或不相似。等同性测试的随机分布是随机置换两进化分支的观测信息，相似性测试的随机分布分别随机构建两进化分支的环境背景，测试使用的随机分布重复100次以确保零假设可以被高置信度拒绝。为检验测试结果是否对不同环境背景都具有相同的敏感性[26]，研究设定3种环境背景来反映喜旱莲子草在自然状态下最大可扩散范围，包括设定0.25°(约25 km)、0.125°(约15 km)的点缓冲区以及去除5%边缘气候的0.25°缓冲区。另外，建模过程中发现使用ecospat包默认的矫正密度函数会导致对应环境信息的密度分布出现边缘偏移，以往入侵种生态位的研究也有指出该问题[12]，因此我们选择使用原始物种分布密度来反映其占据的生态空间。此外，研究还使用100%的环境空间计算了不同入侵种群与原产地种群间的生态位动态指数(稳定、扩展和未填充)，所有指数的值均以百分比形式表示。其中生态位稳定组分反映的是进化分支间共享的气候条件，生态位扩展组分是指入侵种群所占据新的气候条件或气候组合，生态位未填充组分则是入侵范围可利用但尚未被占用的气候条件，以上3种指数的计算参见Guisan等[21]的详细论述。

最后，尽管R软件包ecospat提供了上述比较物种环境生态位的方法(如生态位重叠指数D值)，但这种以温度和降水等相关变量来构建环境空间的方法通常未能考量重叠关系与具体变量之间的联系。因此，研究还通过比较不同种群在主成分环境轴上得分的差异(Kruskal-Wallis测试)，并结合主成分轴对应的主要贡献变量来进一步揭示出进化支间生态特征差异与共性[27]。上述环境空间的生态位分析工作全部基于R软件完成[28]，其中气候环境空间的构建参考使用ecospat包的原始案例[29]和后来由Silva等[30]修改的脚本。

1.4 地理空间(G-space)中比较适生区变化

使用5种不同的方法模拟和预测物种的潜在适生分布区。这些方法包括一种回归方法(广义线性模型，Generalized Linear Model, GLM)，4种机器学习方法(人工神经网络(Artificial Neural Networks, ANN)、随机森林(Random Forest, RF)、广义增强回归建模(Generalized Boosted regression Modeling, GBM)、最大熵模型(MAXENT))。SDM建模的分布数据使用来自各大洲的观测记录，环境数据同样使用1.2节中筛选的8个宏观环境变量，利用biomod2包的BIOMOD_tuning函数分别对上述多个模型调参。调参后的模型各自重复运行5次，使用真实技能统计值(True Skill Statistic, TSS)、KAPPA系数和接受者操作特性曲线下面积(Area Under Curve, AUC)来评估建模的拟合效果[31]，并计算每次迭代过程中环境变量的响应重要性。去除每次迭代运行中单模型评估TSS值低于0.6的模型，将剩余模型构建均值整合模型(Mean Ensemble Model)，使用整合模型可以在一定程度上减少单一模型运行的识别风险。物种分布数据随机分割70%作为训练集，剩下的30%作为测试集,模型输出易于解释0～1连续概率分布用于表征喜旱莲子草的潜在适生区。研究还使用多元表面相似性(Multivariate Environmental Similarity Surfaces, MESS)分析来严格地比较入侵区域与原产地种群的气候相似性，其工作方法是将原产地种群训练环境信息映射到入侵大洲[32]，并保留大洲环境变量相似性值大于零的交集区域来构建与原产地种群气候相似的区域。

2 结 果

2.1 环境空间的生态位变化

参与建模的气候变量在4个大洲的区域种群组间水平上均存在显著差异，组内比较中除大洋洲入侵种群与南美洲种群的日温差(BIO2)无显著差异外，其他所有环境变量均与原产地种群具有显著差异(表1)。其中温度相关变量(BIO2、BIO10、BIO11和GHI)在亚洲、大洋洲和南美洲的边缘种群间呈现出相同的最小值排序格局，即AS

表1 喜旱莲子草原产地种群与入侵地种群气候变量差异检验表

表2 喜旱莲子草原产地种群与入侵地种群气候变量描述性统计

表3 整合模型中建模气候变量的重要性评估表

在0.25°环境背景条件下，标准主成分分析的前两轴占总气候变化的61.32%(图1(b))，其中第一轴主成分主要由BIO16、BIO12、WET等水分变量解释，第二轴主成分由GHI、BIO11、BIO10等温度变量解释。气候生态位主成分分析中原产地种群与亚洲、北美洲和大洋洲种群的气候生态位重叠程度依次降低(表4)，0.125°的环境背景条件的测试结果与0.25°的明显差异。亚洲种群占据着最广泛的气候空间，不仅包含着原产地种群几乎所有的气候范围(SA-AS unfilling： 2.5%),还具有一定程度的生态扩展(Expansion： 29.1%，图1(c))。与亚洲种群相比，北美洲和大洋洲种群未填充比例各自占据原产地种群的61.49%和45.79%。北美洲与大洋洲入侵种群的差异之处在于前者的生态位稳定比例非常高(SA-NA stability： 97.05%)，而后者则存在较高生态位扩展比例(SA-OC expansion： 42.40%)。

表4 喜旱莲子草原产地种群与入侵种群的气候生态比较

生态位重叠值和生态位保守性测试在不同环境集中均相对稳定，不同入侵种群的稳定和扩展比例受不同环境空间背景条件的影响略有差异。入侵种群在0.25°和0.125°环境背景条件下其生态位动态指数均相对稳定，去除5%的边缘气候会导致入侵种群的未填充比例上升和稳定性比例下降。其中，去除边缘气候会导致大洋洲种群的稳定比例降幅达23.65%，这反应出其边缘气候主要与原产地种群重叠的部分相关。另外，生态位保守性测试在所有气候条件下均显著地支持入侵种群与原产地种群生态位不等同，但没有证据表明两者的气候背景显著相似或不相似，因此生态位保守性检验不支持喜旱莲子草在入侵期间发生显著生态位分化[11]。标准主成分轴的测试结果还揭示出不同入侵种群占据的生态空间与原产地种群的生态空间的差异(表5)。与原产地种群相比，亚洲种群的水分变量(BIO16、BIO12、WET)相对保守而生长温度变量(GHI、BIO11)差异显著，北美洲种群的生长温度变量相对保守而水分变量(BIO16、BIO12、WET)差异显著。而原产地种群与大洋洲种群相比，其两个主成分轴同时包含温度和水分变量，并且均存在显著差异，这可能与大洋洲种群的扩展和未填充比例都相当高有关。

表5 整合模型中建模气候变量的重要性评估表

2.2 地理空间的预测分布

使用整合模型构建的SDM预测模型均具有较高的模型性能和预测能力(表6中最小TSS>0.7,最小AUC>0.9)，预测区域与喜旱莲子草实际入侵区域相吻合，其中亚洲和北美洲种群的地理北界均约为36°N。多元表面相似性投影(MESS projection)可以反映原产地种群与入侵区域的环境相似性，3个入侵种群都具有较差的预测映射能力，其中中国境内的实际入侵区域与环境相似性投影差异尤为显著(图2(b)、(e))，这是因为中国境内的总水平辐射(GHI)和最冷季节的平均温度(BIO11)两个温度变量均显著低于原产地种群而到导致该区域的相似性值低于零(图3，见第142页)。大洋洲区域的MESS投影与实际入侵范围也明显不一致，MESS投影区域显示澳大利亚的北部比当前入侵的西南区域存在更适合入侵的气候条件(图2(d)、(g))。另外，结果还发现尽管原产地种群的投影区域与北美洲种群实际入侵区域的内部存在小范围差异，但两者的边界范围基本一致(图2(c)、(f))。研究还进一步比较同属北半球的北美洲和亚洲在关键环境变量(BIO11、BIO12)梯度上的投影差异，如图4(见第142页)所示，亚洲和北美洲的入侵种群均被限制在北纬35°以南的范围内，其中亚洲种群能够在冬季月平均温度低于5 ℃，且湿润季节的降水量高于400 mm的区域广泛分布，而美洲种群北界被较为严格地限制在冬季月平均温度高于5 ℃，且湿润季节降雨量(BIO16)小于400 mm的环境梯度边界范围内。

表6 喜旱莲子草原产地与入侵地种群的整合模型建模优度评估表

图2 喜旱莲子草原产地和入侵地种群的潜在适生区模拟及原产地种群在入侵区域的气候相似性映射Fig.2 Simulation of potential suitable areas of populations innative and invasive areas of Alternanthera philoxeroides and climatic similarity mapping of native populations in invasive areas(a)～(d)分别为使用喜旱莲子草种区域种群(以大洲为单位)分布数据构建整合模型在南美洲(SA)、亚洲(AS)、北美洲(NA)和大洋洲(OC)的平均适生区预测投影；(e)～(g)分别表示亚洲、北美洲和大洋洲与南美洲(SA)原产地种群气候相似的区域。图中连续性适生等级的颜色梯度所代表的概率值如图例所示。

图3 原产地种群气候环境信息与亚洲区域的多元表面相似性(MESS)投影Fig.3 The climate and environmental information of the population of origin (SA) and the multi-surface similarity (MESS) projection of the Asian region每张小图中右下角的小标为使用的环境变量，环境变量的中文释义等同表1注。图中灰色区域表示与原产地种群气候相似的区域，标RMESS的图为前8张图中灰色区域的交集。中国境内的GHI和BIO11与原产地种群气候差异显著，这导致下标RMESS图中中国境内没有与原产地种群相似的气候投影。

图4 北美洲和亚洲入侵种群在BIO11和BIO16两个环境梯度图中的实际地理分布投影Fig.4 The actual geographic distribution of invasive populations in North America (NA) and Asia (AS) in two environmental gradient maps of BIO11 and BIO16 projection(a)，(c)和(b)，(d)分别表示亚洲区域和北美洲区域的环境梯度图。

3 讨 论

本研究比较了入侵物种喜旱莲子草在原产地(南美洲)和3大入侵区域(亚洲、北美洲、大洋洲)的气候生态位特征，与原产地种群相比，入侵地种群的气候生态位都发生了不同程度的变化(图1)。生态位等价性测试和环境变量间的统计检验结果皆显示原产地种群与3大入侵地种群实际发生的气候生态位皆存在显著差异，但生态位保守性测试表明入侵地种群还未发生显著生态位分化(表4)。该结果同样在生态位重叠(D值)和生态位动态指数中得到印证。如北美和大洋洲种群与原产地种群的生态位重叠值都很低(D： 0.2～0.3)，但北美和大洋洲种群的生态位构成存在明显不同，前者的气候生态位近乎全部来自于原产地种群(stability： 97.05%)，而后者则存在高比例的生态位扩展组分(expansion： 42.40%)。原产地种群的生态位是该物种与当地环境(包括非生物与生物)长期相互作用的结果，而入侵种群实际占据的生态位则与入侵历史、入侵区域的环境过滤、生物相互作用、遗传瓶颈及适应性进化等多种因素相关[21]，为此需要一一进行讨论。

环境过滤可能是造成北美入侵格局的主要成因。如图2(c)和图2(f)所示，北美入侵种群的实际分布北界与其原产地种群的气候投影基本一致，这暗示着北美洲的宏观环境过滤可能显著的限制其北向扩散。研究进一步发现，与亚种入侵种群相比，虽北美洲种群地理北界也为0°～5°(BIO11)的低温梯度，但亚洲种群在此低温梯度内分布更为广泛。造成上述差异的原因可能是北美洲种群在低温区域的夏季总降水量(BIO16)低于400 mm，而亚洲种群则高于500 mm，也即低温和夏季降水共同限制了北美洲种群的北向扩散(图4)。

遗传与进化的作用突出体现在入侵亚洲的种群分布上。如前言所述，前期的工作表明喜旱莲子草原产地种群拥有丰富的遗传类型，而入侵亚洲的种群曾经历过严重的遗传瓶颈，遗传多样性非常低，入侵中国的种群甚至可能是同一株系的克隆后代[4-6]。虽然遗传多样性很低，但本研究发现入侵亚洲的种群目前甚至占据着比原产地种群更为广泛的地理空间和生态空间(图1，图2)。如图2(b)和图2(e)的差异和结果所述，亚洲种群能够占据着比原产地更冷的低温环境，从而扩大其气候生态位。近期Luo等[9]的研究表明入侵中国的喜旱莲子草对温度梯度具有很强的表型可塑性，这可能是其能够在亚洲占据更冷生态空间的主要原因。我们推测这一超强的表型可塑性很可能是在入侵过程中获得的，因为已经有研究表明一些入侵种可以通过快速的适应性进化适应新的环境[33-34]，而且喜旱莲子草的繁殖方式确实发生从有性到无性的进化。有性生殖能力的丧失(也意味着有性生殖可塑性的丧失)有可能会通过补偿作用增加无性繁殖的可塑性，这使其能够适应更广泛的温度变化，特别是具有繁殖能力的地下肉质根可以帮助北方陆生种群忍耐更低的冬季低温。

大洋洲的入侵格局可能更加复杂，生态位动态指数显示其未填充和扩展比例都相当高(表4)，目前难以解释这种扩展组分的主导成因。气候相似性分析为揭示复杂格局提供了部分证据，如图2(d)和图2(g)所示，澳大利亚东北部与原产地种群具有相似的气候条件，但实际入侵种群却主要集中在西南区域。进一步观测卫星图发现澳大利亚南西南部有许多重要港口、城市，而回归线以南的东部海岸地带则有较多的森林，且人口密度较低，这意味着更少的人类干扰，所以区域尺度的生物竞争可能是入侵澳大利亚的喜旱莲子草种群不能充分利用根据原产地种群的环境空间所预测的适生区域，从而其生态位拥有较高比例的未填充组分。

入侵物种在入侵过程中生态位的保守性或漂移(niche shift)一直充满争议[35-37]，但这些争议甚少综合考量不同入侵区域对其入侵种群生态位变化的影响。就喜旱莲子草的研究结果而言，其入侵种群的生态位的变化与不同入侵区域的环境过滤、入侵种的遗传特质和生物相互作用等多种因素有关，并且不同入侵区域的地理格局成因也可能各不相同。考虑到入侵大洲的宏观环境过滤(或者说环境可用性)对入侵种群的入侵范围起到主要限制作用，并且其他不同影响因素间具有可交互性，因此这种驱动因素复杂性可能会产生“噪音”，从而很容易掩盖诸如生殖系统和遗传潜力等单一因素的微妙影响[38]，这可能是导致大洲间入侵格局各异的重要原因，也可能是有关入侵物种生态位漂移研究会得出相反结果的原因之一。上述研究结果还提示： 在为入侵物种提供管理和预警政策时，应该综合考量入侵种在不同入侵区域的气候生态位变化以及地理入侵格局，这将有助于更深入挖掘入侵区域的入侵成因和预测当前及未来气候变化下的入侵趋势。