生态环境因素决定了物种的分布, 其中气候因素对物种分布的影响占主导地位, 对生物多样性、遗传多样性以及适应性都会产生深刻的影响^[1]。末次间冰期(Last Inter Glacial, LIG)由于其记录承载体较丰富、分辨率较高、相互之间可比性较强, 使其成为理解目前全新世间冰期气候自然变化的最佳选择对象^[2]。研究表明, 未次间冰期(LIG)以来的气候变化对物种的分布格局有较大影响^[3-6]。联合国政府间气候变化专门委员会(Intergovernmental Panel on Climate, IPCC)^[7]在第五次评估报告中预测, 未来全球气候将持续变暖, 21世纪末全球气温将上升0.3—4.5℃左右^[8]未来气候变化对目前物种的分布范围和分布面积将会产生巨大的影响, 甚至会导致某些濒危物种彻底灭绝^[9]。因此, 了解不同气候情景对物种分布格局变化的影响, 研究物种的潜在分布与气候条件的关系, 能够为物种的起源进化提供地理分布证据, 为引种以及风险防控提供理论依据。

芒属植物是一类C4多年生草本植物, 具有高光效、高产量、强抗逆性和生产成本低等优点。是一类优良的纤维素能源植物, 属内存在丰富的遗传变异, 且自然分布广泛、生境类型多样^[10]。双药芒(Miscanthus nudipes)是芒属植物中适应极度寒冷、干旱环境的一个物种, 适应多种生境类型且具有较高水平的遗传多样性, 对气候变化也有着相应的应对策略, 是一种优质的芒属植物遗传资源^[11]。我国边际土地共有2.90亿hm², 其中年平均气温在-2℃以下的物理边际土地总面积为538.22万hm^2[12], 分析双药芒的适生分布区以及适宜生境, 能为有效开发利用这些物理边际土地提供新的思路, 分析不同时期双药芒适宜面积和分布中心的变化, 能为芒属植物的起源进化提供地理依据。

利用生态位模型分析植物的分布范围、传播趋势、分布面积已成为研究植物与环境、植物与气候关系的常用手段^[13-15]。目前, 对物种潜在分布区分析的生态模型已有多种, 如Maxent最大熵模型, DOMAIN模型, CLIMEX模型, BIOCLIM模型等。在众多生态模型中, Maxent熵模型最具代表性, 李巧云等^[16]运用Maxent模型预测了芒属两个广布种(芒M. sinensis和五节芒M. floridulus)在我国的现代潜在分布区域, 廖莎等^[17]运用Maxent模型预测了中国特有芒属植物南荻(M. lutarioriparia)和芒属广布种荻(M. sacchariflora)在我国的现代潜在分布区域, 但上述研究并未分析地理因素对分布的影响, 以及历史气候和未来气候的变化对芒属植物的分布的影响也未作出分析。本文通过文献查询, 标本信息和野外调研收集分布点, 利用Maxent模型模拟双药芒的不同时期潜在分布区。并利用ArcGIS软件进行适宜区与分布中心点划分, 通过Jackknife检验、响应曲线和环境限制因子研究, 分析不同时期双药芒适宜区变化, 限制双药芒分布的气候因子和适宜双药芒生长及分布的气候条件, 以期为双药芒的生态驯化、优质基因利用和重要种质资源保存利用提供理论依据。

1 实验方法 1.1 分布点数据的获得

通过实地调查记录GPS经纬度, 查阅文献, 查询中国植物数字标本馆(http://www.cvh.ac.cn)、全球生物多样性信息网络(http://www.gbif.org)、PPBC中国植物图像库(http://ppbc.iplant.cn)、国家标本平台(http://www.nsii.org.cn)、昆明植物研究所标本数据库(http://kun.kingdonia.org)、中国植物志(http://frps.iplant.cn)等数字标本平台, 得到中国双药芒现代分布点。将获取的分布点经纬度数据存为csv格式, 导入ArcGIS软件, 检查并删除错误分布点, 保留拉丁名正确和有图片经纬度信息的分布点^[18]。

1.2 气候变量的选取

用于模型建立的气候数据含温度相关(Bio1—11)、降水相关(Bio12—19)、平均月气候数据(Tmin 1—12、Tmax 1—12、Tavg 1—12、Prec 1—12、Srad 1—12、Wind 1—12、Vapr 1—12)和海拔数据(表 1)。运用GisMap软件裁剪气候图层, 转换成ASC格式备用, 本研究采用的当前气候变量从世界气候数据库中下载^[19], 海拔数据从Google Earth中下载, 图层空间分辨率为2.5 min(大约4.5 km)。末次间冰期(Last Inter Glacial, 约120—140ka)、末次盛冰期(Last Glacial Maximum, 约22 ka)、未来气候(2070年)情景采用WorldClim数据库中的CCSM 4模型^[20], 选用未来全球低排放情景RCP 2.6, 为了保持时空序列上模型的可比性, 其海拔和太阳辐射因子在未来潜在分布模拟中保持不变^[7]。地图地理坐标系皆采用GCS_WGS_1984投影坐标系。

1.3 预测模型的构建

将双药芒分布点数据和裁剪后的环境因子图层导入Maxent软件, 设置软件随机选取75%的点进行建模, 25%的点验证模型, 设置最大迭代次数为1000, Bootstrap重复计算10次, 绘制响应曲线(response curves)和刀切法(jackknife)分析环境因子与双药芒分布的相互关系^[3]。将构建好的模型导入ArcGIS软件, 根据适生值百分比的大小, 将双药芒潜在分布区依次划分为以下五个等级：非适生区(适生值介于0—20%)、低适生区(适生值介于20%—40%)、一般适生区(适生值介于40%—60%)、中度适生区(适生值介于60%—80%)和高度适生区(适生值介于80%—100%)^[21]。运用受试者工作曲线(Receiver Operating Characteristic, ROC)下面积AUC(Area Under Curve)值对模型的准确性进行评价^[22]。选用Maxent软件自带的刀切法(jackknife)工具来判断单个环境变量对双药芒分布的影响^[23], 并利用响应曲线(Response curve)评价双药芒的适宜气候条件。运用ArcGIS将不同时期适宜生境分布范围缩小为单个中心点, 计算不同时期分布中心的位置与变化, 分析双药芒的迁移趋势与历史迁移轨迹。

2 结果与分析 2.1 双药芒的分布现状

在湖南农业大学易自力教授课题组10余年来对芒属植物进行资源分布调查和鉴定的基础上, 通过史料查询、实地考察和标本数据库检索的手段, 对双药芒的现代分布点进行系统梳理, 检查并删除错误分布点, 仅保留拉丁名正确且有图片和经纬度信息的分布点, 共计获得61个双药芒有效分布点(图 1)。按省级行政区划分, 主要分布在四川省、云南省、贵州省、西藏自治区和重庆市。其中云南省丽江市、西藏自治区林芝市、四川省雅安市、阿坝藏族羌族自治州、甘孜藏族自治州、凉山彝族自治州分布点最集中。

2.2 气候因子的筛选

气候因子是影响物种分布的主要因素, 本研究选取104个含温度相关、降水相关、平均月气候和海拔数据的环境因子变量。由于变量之间的共线性会导致分布预测模型过度拟合, 采用主成分分析和Spearman相关性分析, 选取相关性小于0.8的变量(表 2), 相关性大于0.8的则保留贡献率最大的变量。经过筛选, 选取9个环境因子用于模型预测：5个气候因子(Bio 3, Bio 4, Bio 6, Bio 7, Bio 15)、3个平均月气候因子(Prec 5, Prec 7, Srad 6)和1个地形因子(Alt)。

2.3 ROC曲线对Maxent模型预测精度检测

ROC曲线是通过Maxent软件自带函数, 对训练集进行模拟得到的, 曲线下的面积为AUC值。AUC的值一般在0—1之间, AUC值在0—0.5之间时, 表明模型预测失败；AUC值在0.6—0.7之间, 表明预测效果较差；AUC值0.7—0.8, 表明预测效果一般; AUC值0.8—0.9, 表明预测效果良好; 大于0.9时, 表明模型预测效果优秀, 越接近1, 表明模型预测效果越好^[24]。本研究运行10次模型, 取平均值, 得到ROC曲线, 其训练AUC平均值为0.987, 标准差±0.002(图 2), 不同的模型重复稳定性好。依据评价标准, 该模型的总体预测精度达到了优秀水平, 说明该模型预测的双药芒潜在适生分布结果准确。

2.4 影响双药芒分布的环境变量

根据Maxent模型对双药芒气候变量贡献分析结果(图 3)占主导地位的气候因子有6个, 分别是温度变化方差、最冷月最低温度、年温变化范围、5月降水量、6月太阳辐射和海拔。其中, 温度变化方差、最冷月份最低温和年温变化范围可归结为气温因子, 5月降水量是降水因子, 6月太阳辐射是光照因子, 海拔则是地形因子。另外, 当使用单独变量进行模型预测时, 温度变化方差的正规化训练增益值、测试增益值和AUC值最高, 其测试增益大于1.7, 正规化训练增益大于1.6, 受试者工作特征曲线下的面积大于0.93, 因此温度变化方差被认为是影响双药芒分布的主要气候限制条件。

2.5 双药芒最适气候条件

为了进一步分析环境因子对双药芒分布的影响, 将6个主要环境因子分别导入到Maxent模型中, 建立单因子模型, 同时绘制出单变量响应曲线, 即双药芒地理分布概率与主要环境因子的关系(图 4)。

一般认为, 当分布概率大于0.5时, 其对应的生态因子值适宜物种的生长^[25]。根据环境因子的响应曲线(图 4)可以得到温度变化方差在440—645时适宜双药芒的生长及分布；最冷月份最低温在-8.5—0℃时适宜双药芒的生长及分布；年温变化范围在22—30℃时适宜双药芒的分布及生长；5月降水量在70—115 mm时适宜双药芒的分布及生长；海拔范围在1630—3750 m时适宜双药芒的分布及生长；太阳辐射 < 15800 kJ m^-2 d^-1时适宜双药芒的生长及分布。

2.6 不同气候情景下双药芒潜在适生区的变化

与末次间冰期(图 5A)相比, 双药芒现代适宜生境(图 5C)范围收缩了9427.80 km², 现代适宜生境范围呈现出向西藏西南部、河南省南部和湖北省北部地区扩散, 其适生中心依旧在四川省中南部地区无明显变化。相比之下现代高度适宜面积缩减了4023.65 km²(占当前生境的0.63%), 高度适宜区由西藏东部地区、四川省西南部向中部地区收缩, 但是西藏的灵芝县也有小面积的高度适生区, 说明在演化过程中, 一部分向西藏地区扩散, 另一部分向四川中部地区收缩。相比末次间冰期双药芒的现代中度适宜面积缩减了538.73 km²(占当前生境的0.09%), 虽然面积缩减了少许, 但是其范围却相对扩散, 尤其是甘肃省南部, 在末次间冰期无中度适生区, 但是到现代却出现了小面积的适生区。而云南的西部也有小面积的扩散。相比末次间冰期现代双药芒的一般适宜面积增加了6245.92 km²(占当前生境的0.99%), 范围向西藏的中南部、云南中部地区扩散, 陕西省更是由无到有, 出现少面积的一般适宜面积。现代低度适宜面积缩减了11111.33km²(占当前生境的1.76%), 陕西省的低度适生区进一步扩增到河南省和湖北省。在末次间冰期云南的南部有一般适生区分布, 但是到现代后变为不适生区。

与末次盛冰期(图 5B)相比, 双药芒当前适宜生境(图 5C)范围变化较小, 仅仅扩增了606.07 km²(占当前生境的0.10%), 现代生境范围向西藏自治区西北部、云南省南部、陕西省北部地区扩散, 而甘肃省的南部、西藏自治区的东北方向、贵州省的东部地区和云南省的东部小范围区域出现收缩。相比末次盛冰期, 现代高度适生区收缩了3114.54 km²(占当前生境的0.49%), 高适宜区的分布中心变化不大, 但分布范围有小幅度的收缩；中度适生面积收缩了11094.50 km²(占当前生境的1.75%), 一般适生面积扩增了22929.75 km²(占当前生境的3.63%), 低度适宜面积缩减了8114.64 km²(占当前生境的1.28%)。

与现代适生区(图 5C)相比, 2070年双药芒的适宜生境(图 5D)范围收缩了85439.43 km²(占当前生境的13.51%), 双药芒分布中心南移, 西藏自治区的适生区向东部收缩, 湖北省、河南省和陕西省的适生区消失, 云南省适生区向西南和东南方向扩散, 相对现代适生区, 2070年高度适生区面积扩增了2946.19 km²(占当前生境的0.46%)；中度适生面积缩减了23468.49 km²(占当前生境的3.71%)；一般适生面积缩减了54041.49 km²(占当前生境的8.55%)；低度适生面积扩增了10875.63 km²(占当前生境的1.72%)。双药芒在未来2070年的适生范围预测总趋势为适生范围面积缩减, 说明双药芒的潜在适生范围与全球环境变化具有相关性, 气候变暖导致双药芒的适生区面积逐渐减少。

Maxent模型预测结果显示(图 6)：末次间冰期适宜面积为641612.25 km², 末次盛冰期适宜面积为631578.37 km², 现代适宜面积为632184.45 km², RCP2.6-2070年适宜面积为546745.02 km²。可以看出双药芒的总适生面积在末次间冰期最大, 然后逐渐收缩, 到本世纪70年代, 其适生区将缩小94867.23 km²。

2.7 不同气候情景下双药芒分布中心的变化

双药芒在末次间冰期的分布中心位于四川省凉山彝族自治州木里藏族自治县唐央乡, 坐标为北纬29.06°, 东经100.88°(图 7)。末次盛冰期的分布中心位于四川省凉山彝族自治州木里藏族自治县倮波乡, 坐标为北纬28.62°, 东经101.41°, 双药芒分布中心整体向东南方向迁移。在现代气候条件下的分布中心位于四川省甘孜藏族自治州九龙县乃渠乡, 坐标为北纬28.84°, 东经101.59°, 相较于冰期, 分布中心向东北方向迁移。在RCP2.6—2070气候情景下的分布中心位于四川省甘孜藏族自治州九龙县三岩龙乡, 坐标为北纬28.76°, 东经101.37°, 可以看出双药芒在末次间冰期至末次盛冰期, 中心点向东南方向发生大幅度迁移, 但是末次盛冰期以后, 分布中心变化较小, 坐标变化范围在28.62—28.84°N, 101.37—101.59°E。可能是由于末次间冰期以来气候逐渐稳定, 到末次盛冰期以后气候变化逐渐平稳, 双药芒的分布中心稳定在四川省甘孜藏族自治州九龙县。总体来看未来双药芒分布中心有向冰期分布中心迁移的趋势。

3 讨论

气候模型是根据特定的算法, 通过已有数据预测其物种分布, 不同的算法得到的模型之间会有一定差异, 但是这种方法仍然是评估和预测物种分布格局演化的重要工具^[26], 近年来Maxent模型在植物和动物等领域得到广泛的运用^[27-29]。本文基于Maxent模型采用了19个全年气候因子、84个单月份气候数据和1个海拔数据, 通过主成分分析法筛选出独立贡献率大且彼此间相关性小的气候因子进行建模, 再结合现有的双药芒分布点信息, 预测了双药芒在过去、现代和未来的分布格局, 并讨论影响双药芒分布的气候限制因子问题。研究结果表明, 气候变化会对双药芒的适生地理分布格局产生较大的影响, 这与前人的研究结果一致^[30-32]。

双药芒潜在中国现代适生区域集中分布在西南地区(川渝云贵藏), 这与实际分布点分布范围一致。影响双药芒分布的环境因子主要为气候因子, 温度变化方差是双药芒分布的生境限制因子, 前人在进行模型构建时, 大多只考虑气候因子, 但是在实际考察中发现, 除了气候因子, 地形因子也对物种的分布起作用。因此本研究加入了影响植物生存最重要的地形因子——海拔因子, 进行气候模型构建。根据气候因子的响应曲线可知温度变化方差在440—645, 最冷月份最低温在-8.5—0℃, 年温变化范围在22—30℃, 5月降水量在70—115 mm, 海拔范围在1630—3750 m, 太阳辐射小于15800 kJ m^-2 d^-1为双药芒分布的最适环境因子范围。

第四纪的冰期气候变化对物种历史分布影响重大^[33]。研究结果显示双药芒在末次间冰期至末次盛冰期, 中心点向东南方向发生大幅度迁移, 但是从末次盛冰期至未来, 分布中心趋于稳定。总体来看未来双药芒分布中心有向冰期分布中心迁移的趋势。适宜区域在末次间冰期、末次盛冰期得到大范围的扩散, 这可能是由于末次间冰期温度逐渐回暖, 温度波动范围变小, 并逐步趋于稳定, 导致双药芒的生存能力和适应性大幅度提高, 最终致使其分布范围大面积扩散。这一现象与李文庆等在末次间冰期情景下研究四子柳(Salix tetrasperma Roxb.)的分布变化结果一致^[7]。就本文而言, 从末次间冰期至近现代双药芒的适生面积一直呈现逐渐缩减趋势。根据气候模型的预测, 预计到2070年, 低排放情景下双药芒的潜在适生面积将缩减94867.23 km²(占当前生境的15.00%)。这种植物分布区域随气候变暖而出现递减的现象与张晓玲等人关于茶适宜区分布的研究结果一致^[34]。

双药芒是芒属植物中一个能适宜极端低温和高海拔的种, 由于芒属植物具有自交不亲和, 且种间能够杂交, 产生基因交流, 是改良芒属植物抗寒种质的重要遗传资源, 但是由于在人类活动、农业活动和工业活动的影响下, 导致芒属植物分布范围逐渐缩减, 我们希望能够通过探明双药芒的潜在分布和环境限制因子, 用于珍贵种质的收集保存, 为以后芒属植物在寒冷等极端边际土地的引种与开发提供遗传资源与理论基础。