植被对气候变化的响应一直是全球气候变化和生物地理学研究的热点问题^[1]。气候是影响区域乃至全球范围内物种和植被分布最重要的环境因子, 气候变化是影响物种分布的决定性因素^[2-4]。20世纪以来, 受人类干扰以及自然因素的影响, 全球气候变化以变暖为主要特征^[5]。随着全球气候变暖, 对物种生活史、物种分布、群落组成、植被格局以及生态系统的结构与功能产生了严重的威胁^{[3, 6-7]}。所以, 研究物种的潜在分布区及对未来气候变化的响应, 对生态系统的可持续发展及制定生物多样性保护策略具有重要的作用^{[6, 8-10]}。

物种分布模型(Species Distribution Models, SDMs)是基于物种已知的分布数据和相关环境因子, 对物种分布及其生态需求进行模拟的一种方法^[11], 现已被广泛应用于气候变化对物种潜在地理分布影响的研究^[12-14]。物种分布模型的发展, 始于BIOCLIM模型^[15]的开发和应用, 随后涌现出HABITAT、DOMAIN、生态位因子分析模型(ENFA)、马氏距离(MD)、边界函数方法(BF)、最大熵模(Maxent)、广义线性模型(GLM)、广义加法模型(GAM)、分类与回归树模型(CART)、推动回归树模型(BRT)、多元适应性回归样条(MARS)等基于统计的和基于规则集的遗传算法(GARP)、人工神经网络(ANN)等模型^[16-19]。其中, Maxent模型具有稳定性和扩展性。基于Maxent最大模型和ArcGIS对物种的潜在地理分布进行预测, 从而使研究者和决策者能提出相应的策略来降低气候变化对生物多样性造成的影响^[20]。

桫椤(Alsophila spinulosa), 是桫椤科、桫椤属蕨类植物, 有“蕨类植物之王”赞誉, 有“活化石”之称。桫椤曾是地球上最繁盛的植物, 与恐龙同期。由于地质变迁和气候变化, 特别是第四纪冰期的影响, 桫椤大都罹难, 只有在极少数的“避难所”才能寻其踪迹^[21]。如今, 残存的桫椤受森林破坏、人为砍伐或自然枯死的影响, 其种群自然繁衍越来越困难, 分布区也大幅度收缩, 已处于濒危状态。在1999年, 我国将桫椤科所有物种列为国家二级重点保护植物^[22]。然而, 目前对桫椤的研究主要集中在生态位^[23-24]、群落特征^[25-27]和毛竹(Phyllostachys edulis)干扰^[28]等方面, 少见桫椤对气候变化响应的报道。本研究结合泸州市、赤水市内自然保护地内大量样点信息和环境因子, 基于GIS运用最大熵模型(Maxent)模型预测桫椤的适宜分布区, 为赤水河地区野生桫椤种质资源的保护和在适宜分布区进行人工繁育提供理论依据。

1 研究区概况

研究区域地处四川省泸州市和贵州省赤水市, 地处四川盆地盆边南部中山植被地区。属中亚热带湿润季风气候区, 雨量充沛, 立体气候和地区差异显著的气候特点。区域内主要有长江、赤水河和习水河三大水系, 溪沟支流众多。年平均气温为18.1℃, 年均降雨量1000 mm, 年均相对湿度82%, 无霜期340-350 d。土壤类型主要为黄壤、山地黄壤、山地黄棕壤。局部河谷焚风影响较显著, 气候干热、降水较少。其主要植被类型有亚热带低山和中山偏湿性常绿阔叶林, 亚热带低山常绿针叶林和亚热带丘陵低山竹林^[29]。

2 研究方法 2.1 样地设置

参照巴拿马Barro Colorado Island(BCI)样地建设技术和经验^[30], 以《生物多样性观测技术导则陆生维管植物(HJ710.1-2014)》的样地建设规范, 结合林分和地形等因素进行样地设立, 共设置20个20 m×20 m的样方, 对样方内所有胸径大于等于1 cm(DBH≥1 cm)的乔木和灌木进行记录登记(不考虑芭蕉、竹类为草本植物和桫椤为蕨类植物), 以胸径为依据对每木进行观测, 记录物种名称、株高、胸径、冠幅等, 并记录样地的郁闭度、海拔、坡度、坡位等环境因子^[26]。

2.2 数据来源及处理 2.2.1 桫椤分布点数据及预处理

本研究主要通过样线踏查、文献资料获取桫椤分布点数据。共在研究区域内布设了110条调查样线, 总长度约210 km, 海拔跨度135-1257 m, 基本覆盖了桫椤在赤水河流域分布的主要海拔段。共得到80个桫椤分布点位, 20个样地信息, 最后把20个样地信息转换成ASCII格式待用。

2.2.2 环境变量数据及预处理

根据桫椤的生境, 选择21个气候变量和1个海拔变量共计22个环境变量^[31-32]作为模型的备选变量(表 2)。其中, 21个气候变量从WorldClim数据库(http://www.worldclimorg)下载获得；海拔变量由中国科学院计算机网络信息中心地理空间数据云平台下载30 m分辨率的DEM计算所得^[33]。在ArcGIS 10.6中, 将22个环境变量(表 2)图层的栅格大小、坐标系统一, 并将图层边界统一, 再将环境变量图层转化为Maxent模型软件需要的ASCII格式。

2.2 模型过程

将桫椤分布点数据和筛选得到的22个环境变量数据导入Maxent 3.4.1中, 随机选取75%的分布点作为训练数据集用于建立模型, 剩余25%的分布点作为检验数据集, 其余设置为默认, 用于模型验证。为保证模型结果的稳定性, 进行10次自举法重复。以软件内建的变量贡献分析、响应曲线和jackknife检验分析模型中各环境变量的相对重要性及其对桫椤生境适宜性的影响。以受试者工作特征曲线(ROC)下的面积(Area Under Curve, AUC)对模型结果的优劣进行评价, 其评判标准为：AUC在0.5-0.6为Fail；0.6-0.7为Poor；0.7-0.8为Fair；0.8-0.9为Good；0.9-1为Excellent^[34]。

2.3 生境适宜性划分

根据10次重复后Maxent输出的平均生境适宜度指数(Habitat suitability index, HSI)对桫椤的生境适宜性进行评价^[35]。在ArcGIS 10.6中对预测分布图以自然间断点分级法(Jenks)进行重分类, 将生境分为不适生区、低适生区、中适生区和高适生区4类。

3 结果与分析 3.1 模型精度评价

根据平均ROC曲线评价结果(图 3)显示：桫椤生境预测模型训练数据AUC值与验证数据的AUC值均为0.999, 明显大于随机分布的AUC值(0.500), 由此, 表明Maxent模型的预测结果为“Excellent”, 预测结果并非随机, 拟合程度相对较高, 结果是可信的, 可以用于桫椤的生境评价。

3.2 桫椤分布与环境因子响应分析

物种的分布区主要取决于其主导环境因子, 在Maxent模型刀切法(Jackknife)分析各个生物环境因子对分布增益的贡献结果中表明(表 3), 主要环境因子及其贡献率为：Bio_2(18.21%)、Bio_4(12.69%)、Bio_14(6.58%)、Bio_15(6.87%)、Bio_18(39.65%)、Alt_2.5(6.82%), 以上6个因子总的贡献率达90.81%。其中, 又以Bio_18、Bio_2、Bio_4和Bio_15(6.87%)为最主要的影响因子, 贡献率达77.42%。

通过对4个最主要的影响因子的响应曲线分析可以发现桫椤对不同环境因子的偏好性(图 4)。桫椤偏好于温度日平均范围(Bio_2)偏低的生境, 生存概率>0.5的适应范围约为5.8-6.8℃, 其最适生境Bio_2值为6.5℃左右, 昼夜温差太高或太低均不适宜桫椤分布；温度季节性(Bio_4)在750左右对其生境影响最为显著, 生存概率>0.5的适应范围约为700-830；桫椤对最暖季的降水量(Bio_18)有明显偏好, 500 mm时达到峰值, 为桫椤最适分布条件, 其生存概率>0.5的适应范围约为400-600 mm；降水量季节性变异性系数(Bio_15)在65左右对其生境影响最为显著, 生存概率>0.5的适应范围约为60-70。

3.3 生境适宜性评价

生境适宜性是指生境能支持一个特定物种的潜在能力, 对物种的生存和发展有着极为重要的意义。根据Maxent模型预测结果(图 5), 桫椤在赤水河地区的适生区的分布范围在空间上存在明显连续性, 从高适生区到不适生区成辐射状。模型预测高适生区主要集中在四川(泸州、宜宾、自贡)、重庆、贵州(遵义), 这与实地调查的结果基本一致, 也与文献资料已知桫椤种群在西南地区的集中分布区相似。

3.4 桫椤高适生区土地利用类型的时空变化

从1980-2018年高适生区土地利用类型^[36]的变化可以看出, 近38年以来土地利用类型面积变化最大的为林地、草地和建筑用地。林地的正增长与草地的负增长的主要原因为天保工程、退耕还林工程的稳步推进^[37]。建筑用地的增加主要是城镇化的推进使城市规模不断扩大。

根据桫椤在赤水市、纳溪区、合江县、宜宾市的生境现状来看, 林地的增加主要是竹林面积的增加。竹类是需水性、再生性很强的植物, 生态位比桫椤更宽, 适生性更强^[24]。桫椤的适生生境正在被侵占、萎缩。

4 讨论

目前对动植物生境进行适宜性评价的模型主要包括3种：生态位模型、机理模型和回归模型, 相比后两者, 生态位模型会更注重考虑物种出现位点和环境变量因子, 在一定程度上生态位模型的可信度更高^[38-39]。而Maxent模型是生态位模型中的一种典型代表。

4.1 影响桫椤生境适宜性的环境因子

根据张华等^[40]、陈爱莉等^[41]和王绮等^[42]分别对胡杨、紫楠和毛榛研究发现, 植物对降水、温度的变化响应特别敏感。桫椤喜生存在温暖、潮湿、荫蔽、水分充足、土层肥厚和排水良好的环境中。基于Maxent模型的桫椤生境质量评价表明, 在赤水河地区桫椤的生境质量非常高, 且具有较好的连通性。桫椤的生境质量受多个环境因子的共同影响, 其最适生境主要分布于最暖季的降水量约500 mm, 温度日平均范围约6.5℃, 海拔800 m左右, 最干月份降水量约20 mm的区域内。结果表明：虽然, 桫椤为高大的树状蕨类, 但其孢子萌发、生长过程对降水、温度较为敏感, 这与蕨类植物喜阴湿温暖的环境相符。

4.2 桫椤实际分布区与预测分布区生境适宜性评价

植物志记载桫椤生于山地溪傍或疏林中, 在国内当前主要分布在贵州、云南、四川、重庆和广西等地^[43-44]。本研究显示桫椤高适生区集中在四川(泸州、宜宾、自贡)、重庆、贵州(遵义)等区域, 本区域主要为丹霞地貌, 特殊的地址构造为桫椤的繁育提供必备条件^[45], 高适生区面积高达49842 km²。但生境破碎化, 人为干扰, 人工植被^[46-47]的影响往往会对桫椤种群造成不利影响。据成晓霞^[28]的研究, 随着毛竹干扰强度的增加, 群落中桫椤和其他树种的立木度逐渐降低, 桫椤群落的总立木数(含毛竹)减少, 物种组成趋于简单化, 幼株缺乏, 种群衰退趋势明显。在毛竹群落中桫椤根系生物量密度均远低于毛竹, 且二者呈负相关^[48]。

4.3 桫椤高适生区土地利用类型的时空变化与生境现状

赤水河地区的丹霞地貌生境为桫椤的种群延续提供了天然的条件, 这一区域零星分布众多桫椤种群。

通过调查研究桫椤的生境现状发现, 竹林生境中的桫椤多生于林缘, 郁闭度较小的阔叶林中的桫椤可与阔叶树种共存。本次调查研究在合江县佛宝自然保护区内的溪沟边红椿林林缘发现两片较大桫椤种群, 可见桫椤喜生水源充足、郁闭度小的林分中。随着天保工程、退耕还林的稳步推进, 长江上游的植被覆盖率得到明显提升^[49-50]。赤水河地区的造林树种多以毛竹为主, 这对桫椤来讲可能导致其种群衰退, 如若不加以干预, 长久来讲可能导致毛竹林中的桫椤绝灭。

5 结论和建议

根据22个历史环境变量数据预测显示, 桫椤在现代气候条件下的潜在地理分布主要集中位于云、贵、川。限制桫椤生长的重要环境因子变量分别为：温度日平均范围、温度季节性、最暖季节降水量和降水季节性变异性系数, 适宜桫椤生长的环境因子变量范围为温度日平均范围为5.8-6.8℃、温度季节性为700-830、最暖季降水量为400-600 mm和降水量季节性变异性系数为60-70。赤水河区域桫椤的适生区范围广泛, 连通性较好, 能够为桫椤提供良好的避难所。但本区域内的阔叶林逐渐被大面积的毛竹林取代, 而桫椤为树状蕨类植物, 植物体的生长发育维持需要消耗大量的养分, 负责供给养分的地下部分却多为柔软纤细, 入土不深的不定根^[51]。在根系发达的针阔树种、毛竹的地下竞争中处于劣势^[52], 种群更新低、衰退趋势明显^[28]。加之桫椤分布的海拔较低, 人为活动可达性较高；成年株限制了孢子繁育的繁育特点影响^[53]及未来气候变化趋势的不确定性。所以, 虽然本区域内桫椤种群数量相对较多, 但其生存现状不容乐观。因此, 应在孢子成熟期人工收集桫椤孢子, 在适生区内水热条件充足的区域开展人工繁育；在受毛竹中度及重度干扰及针阔混交林郁闭度较高的桫椤林中适度间伐、疏伐, 避免桫椤的适生区被进一步侵占；根据实际加强普法教育, 识别、保护桫椤, 并在桫椤集中分布区设立警示标志和围栏或建立保护小区等保护措施。