2019年7月18日, 世界自然保护联盟(IUCN)更新了濒危物种保护“红色名录”, 将超过7000种动、植物列入其中, 并警告称：人类对大自然的破坏造成物种以“前所未有”的速度濒临灭绝。2019年生物多样性和生态系统服务政府间科学与政策平台(IPBES)全球评估报告指出, 当前全球正面临“史无前例”的自然衰退和“加速”的物种灭绝率的局面, 保护和恢复自然需要“变革性改变”[1-2]。全球变暖导致物种分布空间向更北及更高海拔地区迁移[3], 保护区成为物种应对气候变化和栖息地丧失的重要缓冲区域[4], 物种潜在生存空间或保护区选择得到学者持续关注[5-6], 全球变暖已成为国际社会共识。

桫椤(Alsophila spinulosa (Wall. ex Hook.) R. M. Tryon)是隶属于桫椤科(Cyathenaceae)桫椤属(Cyathea Sm.)的一种树形蕨类植物, 由距今1.8亿年前的桫椤类植物演变而来, 是在第四纪冰川劫难中幸存下来的孑遗植物, 对研究古植被演化和蕨类植物系统发育、侏罗纪古环境、地质变迁等具有重要参考价值[7-9], 对维持森林生物多样性亦有重要贡献[10]。桫椤植株多见于中国云南、贵州、四川、重庆、西藏、广东、广西、海南、福建、台湾等省区, 现为国家二级保护植物[11]。桫椤具有先锋树种的特征, 多分布于路边、溪边、沟谷等开阔地带, 正遭受人类毁灭性的干扰, 人工针叶林中适当地梳理林下灌木有利于该种群发展[12-13], 减少过度采伐森林等行为有利于降低桫椤局部灭绝的风险[14]。国内现有研究主要关注典型生长区内桫椤多样性与群落结构特征[8, 15]、生长环境[16]及地理分布[9]等, 关于气候变化引起桫椤潜在适生区分布及其变化规律的相关研究则未见报道。

最大信息熵(Maximum Entropy, MaxEnt)模型是一种基于最大熵理论提出的生态位模型[17-18], 由已知物种分布栅格的环境变量计算约束条件, 在研究区所有栅格内寻找该约束条件下最大熵的可能分布, 具有在较少已知分布点情况下得到较满意预测结果的优点[19], 是目前预测物种潜在分布较好的模型[20], 在预测物种地理分布方面应用广泛[3, 21-25]。本文以中国现有桫椤分布资料, 运用MaxEnt模型预测过去、现在和未来气候条件下桫椤在中国的潜在适生区, 揭示桫椤潜在适生区环境特征、空间分布、变化规律及影响因素, 为保护桫椤应对气候变暖、引种栽培和保护区选址提供参考。

1 材料与方法

1.1 数据来源

整理全球生物多样性资讯机构(Global Biodiversity Information Facility, GBIF, www.gbif.org)、中国数字植物标本馆(Chinese Virtual Herbarium, CVH, www.cvh.ac.cn)、中国知网(www.cnki.net)收录的关于桫椤研究、发现等资料, 保留这些资料中明确标记有桫椤所处地经纬度的信息, 或通过Google earth获取资料中无经纬度但精确到乡(镇)、村地名的坐标；为了避免出现模型过度拟合, 在1 km环境格网数据内只保留一个分布数据, 最终得到现今桫椤分布点189个(图 1)。

图 1 中国桫椤现今植株和第四纪孢粉分布图

Fig. 1 The distribution of A. spinulosa in the present and its spores in the Quaternary

图选项

气候数据来源于世界气候网(WorldClim, www.worldclim.org), 共19个气候因子(表 1), 坐标系为WGS84, 空间分辨率为30 seconds(在赤道地区约为1 km)。包括现今气候(1970—2000年气象观测数据平均值)、末次间冰期(Last Interglacial, LIG, ~120000—140000 cal. yr BP)、末次冰盛期(Last Glacial Maximum, LGM, 21000 cal. yr BP)、中全新世(Mid-Holocene, 6000 cal. yr BP)和2070年模拟气候数据, 模拟气候数据为美国大气研究中心(NCAR)开发的CCSM4模型模拟得到, 该模型对气候变量的模拟和预测准确性都比较高[26], 其中2070年气候数据为CCSM4模型依据典型浓度路径(Representative Concentration Pathway, RCP)2.6模拟的2061—2080年数据产品平均值, RCP 2.6为最乐观温室气体排放情景。地形数据来源于地理空间数据云(Geospatial Data Cloud, www.gscloud.cn)网站, 包括数字高程模型(Digital Elevation Model, DEM)和坡度(Slope), 空间分辨率为90 m, 在ArcGIS(Version 10.2)软件中进行镶嵌, 并用聚合法重采样为1 km空间分辨率, 文中假定末次间冰期以来地形因子不变。地图资料来源于国家基础地理信息中心(http://www.ngcc.cn)。

表 1 环境变量

Table 1 Environmental variables

变量variables

环境变量说明Description of variables

bio_1

年平均温度/℃

bio_2

昼夜温差月均值/℃

bio_3

等温性

bio_4

温度变化方差

bio_5

最热月最高温/℃

bio_6

最冷月最低温/℃

bio_7

年温变化范围/℃

bio_8

最湿季平均温度/℃

bio_9

最干季平均温度/℃

bio_10

最暖季平均温度/℃

bio_11

最冷季平均温度/℃

bio_12

年降水量/mm

bio_13

最湿月降水量/mm

bio_14

最干月降水量/mm

bio_15

降水量变化方差

bio_16

最湿季降水量/mm

bio_17

最干季降水量/mm

bio_18

最暖季降水量/mm

bio_19

最冷季降水量/mm

高程/Altitude

高程/m

坡度/Slope

坡度/°

—

—

表选项

1.2 数据处理

为避免气候因子过多而模型过度拟合, 先用桫椤分布数据(75%作训练集, 25%作测试集)和19个气候因子为输入参数, 运行MaxEnt(version 3.4.0)模型, 用刀切法(Jackknife)检验各气候因子对桫椤潜在适生区的预测贡献；用ArcGIS软件中Band collection statistics工具对各气候因子进行相关性检验分析, 剔除相关系数(R)高于0.8的2个因子中对MaxEnt模型预测结果贡献较小的变量[13], 文中考虑变量的独立性, 将与其他气候因子相关性均较低(R < 0.8)的变量予以保留。最后筛选出昼夜温差月均值(bio_2)、最冷月最低温(bio_6)、年温变化范围(bio_7), 年平均降水量(bio_12)、最干月降水量(bio_14)、降水量变化方差(bio_15)以及高程、坡度共8个环境变量预测桫椤潜在适生区。

采用受试者工作特征曲线(Receiver operating characteristic curve, ROC)下的面积值(Area under curve, AUC)检验MaxEnt模型预测精度, AUC值在0.5—0.6范围即模型不合格；0.6—0.7, 较差；0.7—0.8, 一般；0.8—0.9, 良好；0.9—1.0, 优秀[27]。MaxEnt模型预测桫椤潜在适生概率(P)取值介于0.0—1.0之间, 运用ArcGIS软件将预测结果划分为5个等级：即0.0—0.1为不适生区(Ⅰ)、0.1—0.3为低适生区(Ⅱ)、0.3—0.5为边缘适生区(Ⅲ)、0.5—0.7为适生区(Ⅳ)、高于0.7为高适生区(Ⅴ)[24-25]。

2 结果

2.1 模型预测精度及关键环境变量

训练集、验证集AUC值均高于0.954, 表明模型预测水平优秀。从各变量对桫椤潜在适生区贡献的大小来看(图 2)：年温变化范围、昼夜温差月均值、最冷月最低温的训练增益值约为1.6；年降水量、最干月降水量训练增益值介于1.3—1.4；降水量变化方差、坡度、海拔等因素训练增益值较小。

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

Fig. 2 Jackknife test of importance on environment variables

图选项

各因素对桫椤潜在适生区预测贡献率表明(表 2), 最干月降水量、年降水量、最冷月最低温等对预测桫椤潜在适生区分布起决定作用。桫椤适生概率高于0.5时对环境需求较严格：昼夜温差月均值低于7.6℃、最冷月最低温为3—10℃、年温变化在21.5—27.5℃；年降水量在1200—2100 mm之间、最干月降水量为16—46 mm、降水量变化方差为57—76；海拔160—800 m, 坡度12—28°。

表 2 各变量贡献率

Table 2 Contribution rate of environmental variables

变量Variables

P > 0.5环境变量范围Variable range of P > 0.5

P > 0.1环境变量范围Variable range of P > 0.1

贡献率/%Contribution rate

累积贡献率/%Cumulative contribution rate

最干月降水量Precipitation of driest month/mm

16—46

10—96

36.9

38.1

年降水量Annual precipitation/mm

1200—2100

800—3100

19.8

56.7

最冷月最低温Minimum temperature of coldest month/℃

3—10

-3—15

15.4

72.1

昼夜温差月均值Mean diurnal range

< 7.6

< 9.8

9.8

81.9

坡度Aspect/(°)

12—28

2—48

9.4

91.3

年温变化范围Annual temperature range/℃

21.5—27.5

8—30.5

6.2

97.5

降水量变化方差Precipitation seasonality

57—76

49—96

2.4

99.9

海拔Altitude/m

160—800

< 2600

0.1

100

表选项

2.2 末次间冰期以来气候变化

LIG至2070年, 桫椤潜在适生区年平均温度、最冷月最低温总体增加, 年温变化范围、昼夜温差月均值总体减小；年降水量总体呈增多, 在Mid-Holocene时期达最大, Mid-Holocene至现今减少, 随后增多；降水量方差总体降低, 最干月降水量先增多, 在Mid-Holocene最大, 之后稍有减少(图 3)。总体上看, 随时间迁移, 桫椤潜在适生区气候变暖趋势明显且年温变化范围缩小, 年降水量总体增多且降水量方差减小。

图 3 不同时期潜在适生区环境变量变化

Fig. 3 Change of environmental variables in potential distribution areas of A. spinulosa in different period

LIG：末次间冰期, Last Interglacial；LGM：末次冰盛期, Last Glacial Maximum；Mid-Holocene：中全新世；Present：现今；Future：未来(2070年)

图选项

2.3 桫椤潜在适生区空间分布及变化

2.3.1 桫椤潜在适生区水平位置变化

预测结果(图 4—图 8)表明, 桫椤潜在适生区主要位于中国约32°20′N以南的热带亚热带陆地区域, LIG时期分布在广西、海南、台湾中央山脉东西两侧、广东至福建沿海一带, LIG至LGM时段, 桫椤潜在适生区空间扩展最大, 向北、向西、向东延伸至四川、重庆、贵州、广西、广东、海南、福建、台湾、云南、西藏等省区；自LGM至2070年, 桫椤潜在适生区水平分布变化总体较小。

图 4 末次间冰期桫椤潜在适生区分布

Fig. 4 The potential distribution areas of A. spinulosa in the Last Interglacial

图选项

图 5 末次冰盛期桫椤潜在适生区分布

Fig. 5 The potential distribution areas of A. spinulosa in the Last Glacial Maximum

图选项

图 6 中全新世桫椤潜在适生区分布

Fig. 6 The potential distribution areas of A. spinulosa in the Mid-Holocene

图选项

图 7 现今桫椤潜在适生区分布

Fig. 7 The potential distribution areas of A. spinulosa in the present

图选项

图 8 2070年桫椤潜在适生区分布

Fig. 8 The potential distribution areas of A. spinulosa in 2070

图选项

LIG至LGM, 桫椤潜在适生区中心点由广东中部向西北方向移动至广西北部且距离最远(图 9)；LGM至2070年, 总体由广西北部向东北方向移动至湖南西部, 但距离较小。

图 9 桫椤潜在适生区和环境变量中心点变化

Fig. 9 Change of central of environmental variables and potential distribution areas of A. spinulosa

图选项

2.3.2 桫椤潜在适生区海拔梯度变化

LIG至2070年, 桫椤潜在适生区垂直分布规律明显, 不同时期各级潜在适生区海拔高差大小排序依次均为Ⅱ、Ⅲ、Ⅳ、Ⅴ, 表明桫椤潜在适生区等级愈高, 桫椤分布高差范围愈小(图 10)。Ⅱ、Ⅲ级潜在适生区平均海拔总体先升高, 在现今达最高值后再降低；Ⅳ、Ⅴ级适生区平均海拔先升高、随后呈波浪形升降趋势。

图 10 不同时期桫椤潜在适生区海拔分布箱式图

Fig. 10 Altitudinal boxplot of potential distribution areas of A. spinulosa in different period

Ⅱ：Ⅱ级潜在适生区(0.1 < P < 0.3)；Ⅲ：Ⅲ级潜在适生区(0.3 < P < 0.5)；Ⅳ：Ⅳ级潜在适生区(0.5 < P < 0.7)；Ⅴ：Ⅴ级潜在适生区(P>0.7)

图选项

2.3.3 桫椤潜在适生区面积变化

现今气候条件下桫椤Ⅱ、Ⅲ、Ⅳ、Ⅴ级潜在适生区面积分别为：66.2×104 km2、44.4×104 km2、26.7×104 km2、4.9×104 km2, 总面积为142.2×104 km2。各时期适生区面积变化趋势(图 11)表明, 自LIG至现今, 随气候总体变暖, 桫椤潜在适生区总面积、Ⅲ、Ⅳ、Ⅴ级适生区面积均增加, 在现今气候条件下达到最大；自现今至2070年, 随着气候进一步变暖, 桫椤潜在适生区总面积、Ⅲ、Ⅳ、Ⅴ级适生区面积均减少。

图 11 不同时期桫椤潜在适生区变化

Fig. 11 Change of potential distribution of A. spinulosa in different climates

图选项

Ⅲ、Ⅳ、Ⅴ级潜在适生区面积随气候变暖在不同海拔梯度(间隔100 m)总体呈先增加(在海拔约300—400 m处达最大值), 随后减少的趋势(图 12)。海拔约400—1700 m范围, Ⅱ级潜在适生区面积先增多, Mid-Holocene至现今减少, 随后增至最多；Ⅲ、Ⅳ、Ⅴ级潜在适生区面积亦呈增多, 至现今达最大值后再减少的趋势；海拔约1700 m以上, 主要分布Ⅱ级潜在适生区及少量Ⅲ、Ⅳ级潜在适生区。数据统计表明：77%—98%的Ⅱ—Ⅴ级潜在适生区主要集中在海拔1200 m以下, 81%—97%的Ⅳ、Ⅴ级潜在适生区不超过海拔800 m。总体上看, 自LGM至现今, 桫椤潜在适生区随气候变暖向高海拔地区迁移趋势明显, 自现今至2070年, 随气候进一步变暖, 桫椤潜在适生区在各海拔区间总体呈减少趋势。

图 12 不同时期不同海拔桫椤潜在适生区变化

Fig. 12 Change of potential distribution of A. spinulosa at different altitudes in different period

图选项

2.4 桫椤潜在适生区变化的影响因素

刀切法检验结果表明, 桫椤潜在适生区主要受最干月降水量、年降水量、最冷月最低温等因素影响, 3因素累积贡献为72.1%。选取满足Ⅱ级潜在适生区的最干月降水量(10 mm)和最冷月最低温(-3℃)条件, 获取二者的空间分布交集, 分析其中心点与桫椤潜在适生区中心点的空间变化轨迹(图 9)。LIG至LGM, 最干月降水量(高于10 mm)、最冷月最低温(高于-3℃)中心点由广东中部沿西北方向移动至广西北部, 桫椤潜在适生区中心点亦由广东中部移动至广西北部；LGM至2070年, 最干月降水量(高于10 mm)、最冷月最低温(高于-3℃)中心点广西北部向东北方向移动至湖南西部, 桫椤潜在适生区中心点亦向东北移动至湖南西部, 且Mid-Holocene、现今和2070年桫椤潜在适生区还受降水量方差高于50分布范围进一步限制(图 4—8)。

3 讨论

3.1 桫椤潜在适生区气候特征

现今气候分析表明, 桫椤Ⅳ、Ⅴ级潜在适生区分布于中国热带至亚热带气候区, 季风带来丰沛的降水[28], 全年无霜期多在260 d以上甚至无冬霜, 日平均气温稳定通过≥10℃初终间日数为300 d(川黔渝交界处260 d), 日平均气温稳定通过≥10℃积温达6000℃(川黔渝交界处5500℃)[29]。总体上说, 桫椤潜在适生区光照充足、热量丰富, 这与贵州[8]、福建[16]、海南[9]、广东[30]、四川[15]、广西[31]等省区已有桫椤保护区或研究区年平均气温、极端低温、年降水范围等气候特征、山地地形及海拔范围相似, 且相对湿度较高(约80%—91%)[8, 15, 30], 桫椤常与其他物种构成常绿阔叶林、针叶林、灌木等林相, 森林郁闭度中等偏高[9, 31]。

3.2 桫椤潜在适生区变化规律

由于LIG最干月降水量明显偏低, 该时期桫椤潜在适生区多位于沿海附近, 这与在福建、深圳等地第四纪地层发现桫椤孢粉具有较好的对应关系[32-33]。LIG至LGM, 最干月降水量增多, 最冷月最低温分布区向西、向北扩展且范围最大, 该时段桫椤潜在适生区面积增长最多。LGM至现今, 最冷月最低温增加, 最干月降水量总体增加, 桫椤潜在适生区向北和向高纬度地区迁移、面积增加, 与适生于温暖潮湿环境的豆梨[23]、喜光及温和湿润的鹅掌楸[34]、喜光照的山茱萸[35]等物种同期潜在适生区面积变化具有相似性。现今至2070年, 最冷月最低温持续升高, 而最干月降水量则呈减少趋势, 温和湿润的桫椤生长环境缩小, 导致桫椤潜在适生区向北移动的同时面积减少, 亦与长苞铁杉[36]、水杉[37]等孑遗植物以及珙桐[5]、豆梨[23]等潜在适生区变化具有一致性。

同一海拔梯度内各级潜在适生区面积随气候变化总体呈先增多, 在现今气候条件下面积达到最大随后减少的趋势。已有研究结果表明, 桫椤主要分布于200—1700 m(云南、西藏局部地区达2000 m)低海拔至中高海拔地区郁闭度中等偏高的林中[7, 9, 38], 本文分析结果表明, LIG至2070年, Ⅱ级潜在适生区分布数量在海拔约1600—1700 m内变化明显(图 10、图 12), 这可能反应中国亚热带山地海拔约1600—1700 m范围与落叶阔叶林带、常绿落叶阔叶林带划分最相关[39]。LIG至现今, 最冷月最低温、年均温均呈升高, 最干月降水量、年降水量总体增多(图 3), 促进桫椤潜在适生区在垂直梯度上扩展, 同一海拔梯度内各级潜在适生区面积总体增加。现今至2070年最冷月最低温、年均温均持续升高, 年降水量增多但最干月降水量减少, 可能导致Ⅳ、Ⅴ级潜在适生区减少, 或退化成边缘适生区。

3.3 桫椤潜在适生区变化的影响因素与物种多样性保护

水分与热量条件是对桫椤种群特征影响最大的因子[40], 在水平方向上, 最冷月最低温和最干月降水量影响桫椤潜在适生区分布, 垂直方向上, 海拔对桫椤群落分布起主要作用[41]。

桫椤种群具有一定的自然更新能力[15], 与现存的以裸子植物或被子植物为优势种的生态系统相比, 桫椤的生物量和生产力都比较小, 具有进化上的保守性[42]。入侵物种根系占据土壤空间获取水分及养分, 亦可通过产生对周边桫椤有消极影响的化学物质来抑制桫椤根系生长, 威胁桫椤正常生长发育甚至使桫椤立木度减少[43]。桫椤潜在适生区中心点与环境因子中心点移动趋势表明, 桫椤适应环境变化能力相较于环境变化有一定滞后性。现有桫椤种群中存在以中龄林为主, 幼株数量严重缺失的现象[8], 造成桫椤群落物种结构单一, 种群生境狭窄[41]的现象。桫椤物种多样性保护需在辩证统一地认识物种协同生存、正确认识资源利用与维护、桫椤保护区多元文化与生物多样相互依存的前提下, 适度干预自然生态系统[44], 建议适当扩大生长型桫椤分布区保护范围, 就地对幼株数量严重缺失区域进行幼苗抚育, 以促进桫椤种群恢复、应对全球气候变暖生存空间变窄的问题。

3.4 不确定性

本文基于MaxEnt模型预测桫椤潜在适生区, 训练集、验证集AUC值均高于0.954, 表明模型预测水平优秀。实际生存环境中, 物种相互作用、植被类型、地貌特征、物种自身扩散能力等会对被预测物种潜在适生区产生重要影响[45]。MaxEnt模型预测结果是物种分布的最大可能性, 更侧重于揭示物种潜在地理分布及气候特征, 无法准确表达物种现实分布区域[22]。

从气候数据来看, LIG气候数据虽经过验证, 亦存在一定误差[46], CCSM4气候模拟系统较CCSM3有更高的温度分辨率[47], 但模型误差和内部气候变化有时难以厘清[48]。

现今桫椤在中国分布较广, 现有资料大多记录桫椤所处的县级行政区或保护区, 经纬度则精确到度, 而缺少更精确的坐标信息, 在获取精确到乡(镇)、村位置的桫椤坐标过程中, 虽尽量选择更为湿润的沟谷地带, 但依然存在一定误差, 导致预测结果的不确定性。文中假定末次间冰期以来地形不变, 亦对结果产生一定影响。

4 结论

本文利用现今气候观测数据和过去及未来气候模拟数据, 采用MaxEnt模型分析末次间冰期以来桫椤潜在适生区在中国的分布、影响因素和变化规律, 获得以下主要结论。

(1) 桫椤潜在适生区主要位于中国约32°20′N以南的热带亚热带海拔约1600—1700 m以下山地。LIG时期桫椤主要分布于广西中部、海南、台湾、广东福建沿海一带；LGM、Mid-Holocene、现今以及2070年气候条件下桫椤潜在适生区总体上较为相似, 分布于四川、重庆、贵州、广西、广东、海南、福建、台湾、云南、西藏等省区, 其中2070年潜在适生区在山东省烟台、威海、青岛等地有零星分布。

(2) LIG、LGM时期最冷月最低温、最干月降水量限制桫椤潜在适生区水平分布, Mid-Holocene、现今和2070年, 降水量变化方差进一步影响其适生区空间分布；垂直方向上, 海拔是限制桫椤群落分布的主要因素。

(3) LIG至LGM, 桫椤潜在适生区中心点由广东中部向西北移动至广西北部且距离最远, LGM至2070年, 总体由广西北部向东北移动至湖南西部。

(4) LIG至2070年, 潜在适生区面积、Ⅱ—Ⅴ级潜在适生区面积总体增大, 在现今气候下达到最大, 然后略有减少；Ⅱ、Ⅲ级潜在适生区平均海拔总体升高, 在现今气候下达到最大, 之后降低；Ⅳ、Ⅴ级潜在适生区的平均海拔先升高、随后呈波浪状变化；Ⅱ—Ⅴ级潜在适生区面积在海拔约300—400 m处最大, 约77%—98%的Ⅱ—Ⅴ级潜在适生区分布于海拔1200 m以下, 约81%—97%的Ⅳ、Ⅴ级潜在适生区低于海拔800 m。