基于MaxEnt模型预测未来气候变化情境下红树秋茄(Kandelia obovata)在中国潜在适生区的变化

应邦肯,田 阔,郭浩宇,杨晓龙,李伟业,李 启,骆宇晨,张秀梅,*

1 浙江海洋大学水产学院,舟山 316022 2 浙江省舟山市水产研究所,舟山 316000

了解物种的地理分布与环境因素之间的关系是当前生态学领域的重要科学问题[1-2]。政府间气候变化专门委员会(IPCC)于2014年发布的第五份评估报告表明,未来全球将持续变暖,到本世纪末(2100年)地球的最低气温将比1986-2005年的年平均气温升高0.3-4.8℃[2]。剧烈的气候变化,将会导致许多物种的地理分布发生变化[3]。因此,评估气候变化对物种的分布影响,有助于管理者应对与物种分布范围变化所产生的相关挑战[4]。

物种分布模型(Species Distribution Models,SDMs)是根据物种的实际分布来估计地理空间中物种分布的一种方法,是预测物种分布与环境间关系的常用手段。基于不同的算法规则及预测目的,SDMs衍生出多元化的预测模型,如最大熵模型(Maximum Entropy Model,MaxEnt)、生物种群增长模型(CLIMEX)、基于生物气候数据的生物气候和域模型、生态位因子分析模型(Ecological Niche Factor Analysis,ENFA)和遗传算法模型(Genetic Algorithm for Ruleset Prediction,GARP)等[5-6]。其中由Phillips于2004年提出的MaxEnt模型,应用最为广泛[7-9]。MaxEnt具有优异的预测精度,尤其是在缺乏物种分布数据的情况下,可以通过主要生态环境因素的筛选定量描述物种的潜在栖息地,以现实物种栖息地分布的模拟。因此,在预测入侵物种和气候变化下物种的潜在分布方面有较为广泛的应用[10-13]。MaxEnt的原理是找出最大熵的概率分布,所生成的栖息地适生指数即为物种发生概率,概率的精度往往受模型参数的限制;利用不同模型预测物种分布范围和变化,是目前该领域的一个研究热点[14-16]。

红树林是自然分布于热带潮间带海陆交汇区域的植物群落,是重要的海岸湿地资源;秋茄是红树中最耐寒的品种之一,也是目前向北移栽的主要品种[17]。红树秋茄自然分布的北界在福建福鼎(27°20′N),人工引种的北界为温州乐清(28°13′N)。但随着全球变暖加剧以及我国生态湿地保护修复力度增加,红树能否大规模向北移植,正日益成为红树林生态系统保护领域关注的焦点[18-21]。然而,其栖息地的分布以及该物种在适应气候变化下的可持续性分布北界仍不清晰,且缺乏理论依据支持,这也导致了红树向北移植能否成功存在较大争议。

SDMs的应用始于陆地物种,但它们在海洋物种中的应用频率较低,尤其是在海陆过渡区域,由于环境条件较为复杂多变,影响了其对海洋物种的预测精度[22-23]。本研究基于秋茄在中国及周边国家和地区当前的分布数据,结合ArcGIS和MaxEnt模型,筛选了14个陆地生物气候变量和2个海洋表层环境因子,通过调整海洋与陆地变量的数据分辨率,使其一致,再将两组数据进行重合堆叠,取其有效数据的重叠部分即为海陆交汇区域,模拟了秋茄当前潜在分布区域的气候变化,并预测了2050年和2100年,秋茄种群在我国北界的可能适生区域,以为未来在该区域利用红树秋茄进行生态系统修复和增加碳汇储备提供参考依据及理论基础。

1 材料与方法

1.1 分布位置数据来源与处理

秋茄的分布数据主要来源于全球生物多样性信息网(https://www.gbif.org(访问日期,2022年03月29日))、中国植物标本馆(http://www.cvh.ac.cn(访问日期,2022年03月30日))中已有的秋茄分布数据,并查询了关键字索引“秋茄”的文献256篇,形成了秋茄的分布数据库。其次,为防止同一范围内出现重复数据,影响模型预测精度。我们基于世界气候数据库(www.worldclim.org)的最高分辨率(2.5arcmin,约为5km),对单个最高分辨内的重复数据进行了删除,保证每5km范围内只存在一个分布数据。最终,共获取了141条秋茄的分布数据[24]。

1.2 环境变量及处理

1.2.1陆地环境变量及处理

本研究使用WGS84坐标系和2.5弧分分辨率,从世界气候数据库(www.worldclim.org)中选择了19个生物气候变量(Bio1-Bio19),作为陆地环境变量进行建模。为避免模型过度拟合,对共线性较高的两个变量(Pearson检验,相关性≥0.8),只保留其中具有生物学意义,且在初始模型中贡献最大的变量(表1)。

表1 筛选后影响秋茄生境分布的环境因子Table 1 Environmental variables affecting habitat distribution of Kandelia obovata after screening

根据CCSM4模型确定未来的气候数据(2050年代和2100年代)[25-26]。该模型包含IPCC第五次排放报告中的四种排放情景。根据其预测到2100年,全球温升的变化,我们分别选择了CO2低浓度排放路径(RCP2.6,温升<2℃);中等浓度排放路径,与现阶段排放量不变(RCP4.5,4℃);高浓度排放路径,化石燃料密集型(RCP8.5,6.9℃)三种情景[27-29]。在以上三种排放情景中,分别选用2041-2050和2091-2100时间段的预测结果,并使用相同的处理方式对未来环境因子进行预处理。

1.2.2海洋环境变量及处理

由于秋茄生长于陆海过渡生境的特殊性,本研究还结合了海洋表层环境变量对秋茄种群分布的影响进行了探讨。从OceanColor、Bio-ORACLE、Halpern与MARSOEC数据库中筛选了海洋环境数据,其中全球海洋生物扩散模型环境数据库(Bio-ORICLE)具有较高分辨率(5arcmin,约为9.2km)和较丰富的数据种类[30-31]。并且,该数据库对所有环境数据统一了坐标系,并适用于通用的SDMs软件文件格式。因此,本次研究历史海洋环境因子选取全球海洋生物扩散模型环境数据库(Bio-ORICLE,http://bio-oracle.org)中基于2000-2014期间月平均的表层海洋环境数据编制成的栅格格式图层,作为本次研究的海洋环境历史数据与陆地环境因子进行组合。

将海洋环境因子导入ArcGIS 10.8(坐标系WGS1984),导出修正为2.5弧分的栅格文件,保证与陆地环境变量的数据分辨率和数据范围一致。将修正后的海洋环境因子导入ArcGIS 10.8,使用ArcToolBox中的转换工具将所有海洋环境变量的栅格数据转为ASC格式数据文件,随后对ASC源文件对数据范围进行修正规范为坐标系默认范围。

1.3 现有分布位置的修正

将筛选处理后的环境因子与秋茄分布位置数据导入ArcGIS 10.8(坐标系WGS1984),陆地与海洋分布的环境数据重合堆叠的部分即为潮间带区域。为保证秋茄的分布位置处于环境因子认定的潮间带范围内,对秋茄的分布坐标进行了人工修正。同时,为保证数据的准确性,将经纬度精确到小数点后3位。

1.4 MaxEnt模型构建与评价

MaxEnt模型处理问题的方法是将物种分布的预测转化为一个概率模型,使物种分布预测的随机性成为一个概率分布,并找到最优分布概率[32]。熵值的计算结果随着各分布数据相关环境因子的输入和迭代次数的增加而增大。最后得到熵最大的状态,即最接近真实事物的状态。数学上,给定一个随机变量ε,它有n个不同的可能结果X1,X2,…Xn发生的概率是p1p2…,那么ε的熵Pn可以用表示:

本研究采用受试者工作特征曲线法(Receiver Operating Characteristic,ROC曲线)以评估模型精度。ROC曲线下的面积为AUC值(Area Under the receiver operating characteristic Curve),AUC值越大,表示环境变量与预测物种地理分布模型间的相关性越大,越容易分辨该物种有无分布,AUC取值范围在0-1之间,越接近于1,预测效果也就越好[33]。

1.5 MaxEnt模型参数设置及修正

修改基础设置,选定随机取样,随机测试比例设定为25%,重复10次试验降低试验偶然性引起的误差。随后修改高级设置,选择编写地图数据(Write plot data)选定选择输出格式为Logistic,选定运行数据输出内容绘制预测图片(Make pictures of prediction)、创建响应曲线(Create response curves)、使用刀切法测定变量的重要性(Do jackknife to measure variable import),其余参数保持默认[34]。随后输入已修正的csv.格式的分布位置数据与已筛选并修正后的陆地、海洋环境因子数据,设定输出文件夹后运行模型进行试验。

1.6 秋茄适生区域等级划分

根据秋茄现有分布状况,并参照黄晓君等的适生性指数Fitness index(FI)划分方法。本研究采用了自然断点分类法,将秋茄的分布区划分为4级;高适生性(FI≥0.456)、中适生性(0.456>FI≥0.207)、低适生性(0.207>FI>0.051)和不适生性(0.051≥FI> 0),绘制出秋茄潜在适生区分布图(图1)[35]。

图1 秋茄分布样点与潜在适生区Fig.1 The occurrence sites and potential suitable habitats of Kandelia obovata

1.7 利用回归模型对模拟结果进行验证

在用MaxEnt模型得出现有秋茄适生区分布后,根据贡献度和训练增益中筛选出,影响秋茄分布的主子。将ArcGIS中秋茄现今分布点位(n=141)上的主要环境因子和适生性指数提出,利用广义加性模型(R语言 mgcv包和nlme包)对主要环境因子与适生性指数进行回归验证,以评估主要环境因子的筛选是否准确。

2 结果

2.1 MaxEnt模型评估

利用MaxEnt软件模拟基于141个秋茄现分布点和16个环境变量的当前值。10次试验平均AUC值为0.990(±0.001),预测范围覆盖了全部实际分布,表明MaxEnt模型所得到的数据是准确可靠的(图2)。此外,对未来(2040s和2090s)模拟的AUC值均大于0.989,表明MaxEnt模型预测性能良好[36-37]。

图2 基于MaxEnt模型预测秋茄潜在分布的受试者工作特征(ROC)曲线Fig.2 The receiver operating characteristic curve (ROC) predicting the current potential distribution of Kandelia obovata based on MaxEnt model图中AUC表示受试者工作特征曲线下面积,Area under the receiver operating characteristic curve

2.2 秋茄适生区分布的主导因子

根据Jackknife测试法对影响秋茄现今潜在分布的环境因子进行分析。对秋茄潜在分布位置影响程度最深的环境因子(贡献率>10%)按照贡献率排序分别为:海洋表层平均水温(temperature)(39.9%)、等温性(bio_03)(25.9%)、最暖季度降水量(bio_18)(12.5%),累计贡献率达78.3%,表明以上三个环境因子是影响秋茄潜在分布区的主导环境因子(表2)。基于刀切法分析的不同气候因子的训练增益,结果表明,训练增益得分前4的气候因子依次为海洋表层平均水温(temperature)、年平均气温(bio_01)、最热月份最高温度(bio_05)、最暖季度降水量(bio_18)(图3)。本研究综合贡献率排序与刀切法训练增益排序结果,对以上环境因子取并集作为主要影响秋茄分布的目标环境因子进行讨论。

图3 刀切法检验训练数据中各变量的重要性Fig.3 Variables importance of jackknife test in training data 图中All variable 为全部环境变量;bio_01:年平均气温 Annual Mean Temperature; bio_02:平均气温日较差 Mean Diurnal Range; bio_03:等温性 Isothermality; bio_05:最热月份最高温度 Max Temperature of Warmest Month; bio_07:气温年较差 Temperature Annual Range; bio_08:最湿季度平均温度 Mean Temperature of Wettest Quarter; bio_12:年降水量 Annual Precipitation; bio_13:最湿月份降水量 Precipitation of Wettest Month; bio_14:最干月份降水量 Precipitation of Driest Month; bio_15:降水量季节性变化 Precipitation Seasonality; bio_16:最湿季度降水量 Precipitation of Wettest Quarter; bio_17:最干季度降水量 Precipitation of Driest Quarter; bio_18:最暖季度降水量 Precipitation of Warmest Quarter; bio_19:最冷季度降水量 Precipitation of Coldest Quarter; salinity:海洋表层平均盐度 Average salinity of the ocean surface; temperatrue:海洋表层平均水温 Average temperature of the ocean surface

此外,主要环境因子所提供的最大生存概率,其中bio_01(年平均气温)的最低,仅为0.449。其余为最暖季度降水量bio_18(0.866)>海洋表层平均水温temperature(0.839)>最热月份最高温度bio_05(0.703)>等温性bio_03(0.516)(图4)。

图4 主要环境因子对秋茄适生概率的响应Fig.4 Response curve of environmental variables

2.3 秋茄分布区的气候特征

适生性指数≥0.456的区域为模型得到的秋茄的高适生区,与现实分布区吻合度高。因此,本研究利用存在概率大于0.456时对应的气候因子范围,表征秋茄分布区的该气候因子特征[38]。秋茄分布区的气候特点为:等温性(bio_03)为23.43-33.99,最热月份最高温度(bio_05)>31.7℃,最暖季度降水量(bio_18)>740.61mm,海洋表层平均水温(temperature)>24.9℃(图4)。

2.4 秋茄的潜在分布与变化

本次研究使用MaxEnt模型,根据RCP2.6、RCP4.5和RCP8.5排放情景所对应的环境因子,分别对秋茄本世纪中叶和世纪末中国大陆分布北界区域进行了模拟(图5)。

图5 不同排放情景下秋茄适生分布区预测Fig.5 Response curve of environmental variables on the northern boundary

当前秋茄在我国大陆沿岸的高适生区如图1所示,北界在浙江省温州市瓯江入海口附近(27°98′N),基本与我国目前秋茄人工引种分布北界重合,本研究将高适生区认定为秋茄潜在自然分布区[38]。可能引种的区域(低适生区)北界位于江苏省盐城市大丰区三龙镇沿海(33°29′N)。

模型结果表明,在RCP2.6、RCP4.5、RCP8.5气候变化情景下,秋茄的潜在分布区发生集中变化的区域在浙江与江苏一带,以秋茄自然分布区北界变化最为显著,基于低浓度排放(RCP2.6)情景下,2050年与2100年的秋茄高适生区相对现今秋茄高适生区均有所北移,但低适生区有南撤现象。而中高浓度排放情景(RCP4.5、RCP8.5)下,秋茄高适生区与低适生区北界均发生北移。

3 讨论

通过MaxEnt模型了解秋茄的具体分布与潜在分布区对于秋茄的北移引种具有重要的参考意义。本研究预测了在不同碳排放情景下本世纪中叶和世纪末秋茄的北界分布。研究结果表明,高适生区范围与目前已出现大范围人工种植与文献报道的自然分布基本吻合,而低适生区范围也有秋茄耐寒试验的研究[39]。因此,本研究对高、中、低适生区分别定义为,已长期存在的(≥5年),较大范围的秋茄自然分布区与人工种植林区、秋茄可种植试验区、秋茄生长极限试验区,建议以此为标准根据秋茄分布边缘北移变化制定相关移植策略。

3.1 MaxEnt模型的可靠性

模型的精度取决于样本的覆盖度、区域和样本量的多少,而AUC值是模型的最佳衡量指标。本研究收集了141分布数据,尽可能覆盖了秋茄的分布范围[12]。由于秋茄属于典型的潮间带红树植物,为保证模型的准确性,本次模拟结合了海洋表层环境变量。但目前MaxEnt模型在物种分布预测方面仅用于单纯的陆地环境变量或者海洋环境变量,并未有陆地环境变量与海洋环境变量结合讨论的先例,从本次试验最终得到的分布模拟结果来看,所得到的高适生区域与收集到的自然分布点位达到高度一致,并呈现与海岸带重合的条带状结果,这为后期潮间带物种的分布预测提供了良好案例。本次研究所进行的十次试验平均AUC值为0.990(±0.001),并将所得结果中高适生区与收集到的秋茄分布位置进行比对后,发现基本一致,说明使用主导环境因子利用MaxEnt模型来预测秋茄的潜在生境分布情况的性能出色,预测结果精度高,有效避免了过拟合现象,可信度高。

根据贡献度和训练增益中筛选出的主要影响秋茄分布的5个环境因子,基于141个秋茄分布点的适生性指数,用广义加性模型进行了结果验证。5个环境因子中,年平均气温(bio_01),等温性(bio_03),最热月份最高温度(bio_5),海洋表层平均水温(temperature)与适生性指数都有极显著关系(P<0.001)。最暖季度降水量(bio_18)有极显著影响(P<0.01)。验证模型结果显示整体方差解释率高达94.7%(R2=0.915),表明这些变量能够在秋茄的分布概率上解释90%以上的方差,模型拟合结果优良(表3)。

表3 环境因子与适生性指数的显著性分析Table 3 Significance analysis of environmental factors and fitness rate

3.2 适生性指数对环境变量的响应

根据模型结果可知,温度(年平均气温、最热月份最高温度、等温性)、降水(最暖季度降水量)和海洋表层平均水温是影响中国秋茄分布的主要环境因素。目前我国秋茄的自然分布北界位于福建福鼎,属于典型的亚热带区域。相关研究表明红树林的光合作用最适生叶温为28-32℃,张乔林等学者的研究试验表明,温度决定了秋茄的纬向分布,包括气温与水温[40-43]。陈鹭真等发现秋茄的北界可至暖温带,去除气温影响,水温对秋茄分布的影响更为显著,冬季适生水温为10-20℃,夏季适生水温为15-25℃[44-46]。在仅考虑单一环境因子影响下,海洋表层平均水温可提供的秋茄最大生存概率为0.839,而年平均气温能够提供给秋茄的最大适生分布概率仅为0.449,这表明海洋表层平均水温对秋茄的分布有着更加显著的影响。等温性是反映温度变化的迟早和幅度的一个指标,这关系到植物的生长、发育,影响植物的温度敏感度和有效积温。温度敏感度是指温度每升高1℃,植物物候期变化的天数。温度敏感度越高的物种,在群落中的覆盖度和生物量会占据优势。积温是指超过一定起点温度(例如 温度>5℃)的逐日平均温度累积和。温度升高使植物每日积温速率增快,使得植物生长发育的有效积温需求能更快被满足。因此,当等温性指标为23.43-33.99,能大概率的促进秋茄的展叶、发芽和种群发育。

此外,秋茄人工林北界与自然北界均位于季风性气候区域,该地区受台风影响较大,会导致降水突然性增多,降水的异常会导致秋茄单位叶面积的净同化速率减少,使得红树生长缓慢。同时,也会降低当地海水盐度,但秋茄对盐度的耐受性较高,这可能是导致盐度在模拟结果中贡献度较低的原因[45-47]。根据已有的研究结果,漳江口区域在5-8月雨量充沛,最暖季度降水量约900mm,有利于红树林生长发育。而本研究发现秋茄高适生区的条件在最暖季度降水量(bio_18)>740.61mm,这与漳江口区域的研究结果相近。秋茄更宜生长在最热季度期间降水相对较为充沛的地区,浙江区域降水的季节性变化与不稳定性较华南区域更为明显。因此,这可能会导致秋茄在浙江地区难以形成广泛的连续自然分布,移栽时,应对当地的降水量进行充分的考虑。

在回归模型中,edf为有效自由度,当edf=1时,环境因子与适生性指数呈线性相关,edf>1为非线性相关。Ref.df和F是在 ANOVA 检验中使用的检验统计量,用于检验平滑的整体显著性,F值越大,则P值越小,方程越显著,拟合程度越好(表3)。由此可知,主要环境因子与适生性指数都是非线性显著相关(P≤0.001),其中最热月份最高温度 (bio_05)和海洋表层平均水温(temperature)与适生性指数的非线性相关最为复杂,其次是最暖季度降水量(bio_18)和年平均气温(bio_01),最后是等温性(bio_03)。由此,我们可以推测秋茄的分布可能极易受到极端气候因素的影响。

3.3 秋茄的地理分布及变化

根据模型结果(图5),在不同的碳排放情景下,秋茄在中国大陆北界变化差异较大。到2050年,低浓度排放情景下(RCP2.6),自然分布的北界会移至温州苍南(27°50′N)到乐清(27°59′N)附近。可移植的适生区会至江苏省南通市如东县海域(32°26′N)附近,舟山群岛具有分布的可能性。RCP4.5条件下,秋茄自然分布的北界可达浙江省玉环南部与乐清雁荡山附近(28°16′N),至江苏省盐城市射阳县附近(33°41′N)可进行秋茄移栽。高浓度排放下(RCP8.5),预测秋茄自然分布北界与可引种区域,基本与RCP4.5情景下一致。

到2100年,低浓度排放背景(RCP2.6)下,秋茄的低适生区与2050年相比,有所南撤且长江口北岸低适生区会全部消失。RCP4.5情景下,秋茄的自然分布北界可达浙江省瑞安附近(27°46′N),低适生区北界与2050年基本持平,但在江苏省连云港市秦山岛会出现零星分布。RCP8.5情景下,我国大陆秋茄自然分布北界可达浙江省温岭市(28°18′N),并在瑞安市到温岭市一带滩涂湿地区域出现大范围自然分布区域。低适生区可达山东半岛南岸(36°52′N)且在连云港以南范围内连续性出现。

根据姚遥等研究,RCP4.5排放条件与目前全球变暖趋势更为接近,故更能代表未来实际气候环境变化趋势[48]。在中等浓度排放下,2050年高适生区北移相对较为明显,可达到浙江省玉环南部与乐清雁荡山附近海域(26°16′N)。在2100年高适生区北界有所南撤,这可能与秋茄在最热季度需要较多的降水量所影响[49]。RCP4.5背景下,低适生区实现了对舟山群岛的全覆盖,甚至有中适生区的存在,这揭示了未来在舟山区域内进行较大范围的秋茄引种试验将成为可能。目前,秋茄人工引种在浙江舟山岱山(30°18′N)以及江苏南通(32°15′N)进行了试验,并开展了越冬抗寒驯化实验,在舟山岱山县已有移栽案例[39]。此外,该区域的秋茄引种除需要考虑极端天气变化外,因降水导致的低盐度也应受到重视[20]。基于模型的预测,未来在中国大陆可种植秋茄的最北界,可能会达到江苏省盐城市射阳县附近。因此,这进一步揭示秋茄分布有北移的倾向。

MaxEnt模型操作简单,样本需求量小,预测精度高,这对我们理解秋茄的分布和群落的扩散迁移趋势有一定参考价值,但在模型运算中仍有较多局限性。(1)本模型所涉及的变量均为气候变量,而影响秋茄种群的分布除了气候因素外,湿地底质类型、潮汐作用、海平面高程也都是影响其种群分布的重要因素。因此,上述预测结果侧重于展现秋茄种群未来可能的一种分布趋势,以及影响该种群分布的重要气候特征。(2)除环境因素外,物种个体的生长特性、种群扩散、迁移能力以及物种间的相互作用也会影响秋茄种群的分布。(3)基于抗寒特性的新品种选育和移植,也会增加秋茄未来分布的不确定性。

综上所述,(1)研究通过MaxEnt模型分析表明温度、海洋表层水温和降水是主导秋茄种群分布的主要气候因子,同时明确了秋茄适生区的各主要气候因子特征。等温性(bio_03)为23.43-33.99,最热月份最高温度(bio_05)>31.7℃,最暖季度降水量(bio_18)>740.61mm,海洋表层平均水温(temperature)>24.9℃。(2)根据RCP2.6、RCP4.5和RCP8.5排放情景下所对应的环境因子预测,秋茄种群的地理分布较现今有北移趋势,北至长江口附近,甚至到达江苏沿岸。这为上述区域利用红树秋茄进行生态系统修复和增加碳汇储备提供了参考依据。