气候变化、人类活动的干扰及物种自身扩散能力等对物种分布格局具有重要的影响^[1]。在全球气候变暖的大背景下, 气候威胁日趋严重, 预测气候变化引起的物种分布格局变化、物种损失及周转情况已成为生态学的重要研究领域。相关研究发现, 在生物多样性保护、维护生态系统可持续发展等方面具有重要理论价值^[2—3], 也有助于开发合理有效的保护策略应对气候变化带来的潜在风险^[4—5]。东北地区是中国森林分布的主要地区之一, 植物物种丰富度较高, 对未来气候变化的响应十分敏感^[6]。有研究表明, 由于气候变暖和人为活动的影响, 东北地区出现了森林生态系统退化、湿地萎缩、冻结层退缩、沙化土地扩大等一系列生态环境问题^[7], 但是森林关键树种分布格局对气候变化响应, 尚缺乏多树种、多气候模式未来响应研究。

目前, 生态位模型已广泛应用于濒危动植物保护、物种入侵及全球变化对物种分布和多样性格局影响等相关领域^[8—10]。生态位模型主要在模拟物种的实际分布拟合的基础上, 对其在未来气候变化下的潜在分布进行预测。常用的生态位模型主要有Bioclim、Domain、Garp、Climex、MaxEnt等^[11]。其中MaxEnt模型以最大熵理论为基础, 根据物种实际分布点和环境数据确定的函数关系, 作为约束条件寻找熵最大的概率分布, 从而预测物种的潜在分布区^[12]。与其他模型相比, MaxEnt模型仅需要物种出现点数据, 且对样本数量要求较低^[13]。近年来, 许多学者应用MaxEnt模型利用WorldClim数据库v1.4中的生物气候因子, 对未来不同气候情景下的物种适宜分布区、物种丰富度分布格局等进行模拟预测^{[3, 14]}, 使得探究物种地理分布格局对社会经济变化和政策干预下的响应成为可能^[15]。

联合国政府间气候委员会IPCC对未来气候场景的更新, 要求更加深入的模拟分析。2020年WorldClim数据库将全球气候和天气数据从v1.4和v2.0更新至v2.1^[16], 把历史气候数据从1960—1990年代更新到了1970—2000年代。在v1.4中, 未来气候包括4个具有代表性的CO₂浓度路径(RCP), 即rcp26(低排放路径)、rcp45(中排放路径)、rcp60(中高排放路径)和rcp85(高排放路径)^[17], 是IPCC全球耦合模式比较计划第五阶段(CMIP5)内容。在v2.1中, 作为IPCC-CMIP6的内容, 未来气候设置了四种共享的社会经济路径(SSPs), 即SSP126(Sustainability, 可持续路径)、SSP245(Middle of the Road, 中间路径)、SSP370(Regional Rivalry, 区域竞争路径)和SSP585(Fossil-fueled Development, 化石燃料为主发展路径)。与CMIP5相比, SSPs能更好地反映社会经济发展与气候情景间的关系^[18], 囊括了更加具体的未来气候周期, 模拟结果更接近于实际观测^[19]。

基于上述, 本研究以东北地区为研究对象, 收集12个建群树种的分布数据及可能影响其生存分布的环境因素, 应用MaxEnt模型对它们的潜在适生区及主导环境因子进行模拟预测。拟解决以下几个科学问题：1)影响东北地区树种的主导环境因子有哪些？2)东北地区不同树种分布对不同气候情景的响应如何？3)东北地区树种受未来气候变化的威胁状况及其物种损失和周转状况？是否受物种可迁移性的影响？相关研究以期为树种的有效保护提供重要的科学依据, 为东北地区植物多样性保护, 进而更好发挥其生态功能提供重要的理论参考。

1 材料与方法 1.1 物种分布数据获取及处理

本研究选取东北地区12个建群树种^[20](表 1), 其中, 针叶树种包括红松(Pinus koraiensis)、红皮云杉(Picea koraiensis)、兴安落叶松(Larix gmelinii)、樟子松(Pinus sylvestris var. mongolica), 阔叶树种包括山杨(Populus davidiana)、春榆(Ulmus davidiana var. japonica)、白桦(Betula platyphylla)、黄檗(Phellodendron amurense)、水曲柳(Fraxinus mandshurica)、胡桃楸(Juglans mandshurica)、蒙古栎(Quercus mongolica)、辽东桤木(Alnus sibirica)。树种的地理分布数据来源于以下3个途径：1)科技部科技基础性工作专项重点项目和国家自然基金重点项目构建数据库中的野外实测数据；2)iPlant.cn植物智平台(http://www.iplant.cn/)植物图片及标本记录数据, 对于缺少精确经纬度信息的数据, 通过百度地图定位到县或乡镇来确定其地理分布点信息；3)全球生物多样性信息平台(GBIF数据库, https://www.gbif.org/), 剔除研究区域外及重复的分布数据。为避免物种采样点不均衡, 过于集中造成数据冗余, 从而影响模型的预测精度, 利用ArcGIS软件将获得数据分布范围转换为与气候数据分辨率(2.5arc-min)一致的等面积5km×5km的网格单元, 每个网格内仅保留一个数据点, 导出至Excel, 保存成CSV格式, 得到树种当前分布数据。

1.2 环境变量数据获取及处理

本研究共选取23个影响物种分布的环境变量, 包括19个生物气候因子、土地利用类型栅格数据LUCC(2007年的MODIS土地利用产品)及地形因子数据包括海拔Dem、坡度Slope、坡向Aspect。土地利用类型数据及地形因子数据由中国科学院东北地理与农业生态研究所遥感与地理信息研究中心地理景观遥感学科组提供, 空间分辨率均为1km。生物气候因子数据来源于WorldClim (https://www.worldclim.org/), 空间分辨率为2.5arcmin。其中, 当前(1997—2000年)的环境数据为wc2.1, 未来(2021—2040、2041—2060、2061—2080、2081—2100)气候情景下的环境数据为第六次国际耦合模式比较计划(CMIP6)中BCC-CSM2-MR的数据, 选取共享社会经济路径(SSPs)中的3种情景模式：可持续发展(SSP1_2.6)、中度发展(SSP2_4.5)和以化石燃料为主的常规发展(SSP5_8.5)^{[15, 22]}。利用ArcGIS软件提取中国东北地区环境变量数据, 并将其从*.tif格式转换为*.asc格式。

由于19个气候因子间具有相近的属性含义, 可能存在多重共线性, 直接将19个气候因子导入MaxEnt软件进行模型构建, 可能会造成数据过拟合, 影响模拟结果^[23]。本研究通过刀切法检验各气候因子对建模的重要性, 运用ENMTools工具对各气候因子进行相关关系分析, 属性相同的气候因子相关系数绝对值≥0.8, 结合物种生理生态特性选择贡献率最高的一个, 最终选取9个气候因子用于模型构建(表 2)。

1.3 研究方法 1.3.1 构建MaxEnt模型

将物种分布数据(CSV格式)和环境变量(ASC格式)导入MaxEnt软件, 随机选取75%的分布数据作为训练集用于建模, 25%的分布数据作为测试数据集用于模型验证^[24], 正则化乘数及其他参数均选择默认值^[25]。通过刀切法(Do jackknife)测定各环境变量对模型构建的贡献率, 创建环境变量响应曲线(Creat response curves)以期获取物种分布概率对各环境变量的响应情况。利用受试者工作特征曲线ROC(Receiver Operating Characteristic)对模型预测的准确性进行评价^[26]。将ROC曲线与横纵坐标围成面积最大值AUC(Area Under Curve)作为模型精准度评价标准, 取值范围为0—1, AUC值越接近于1, 表明模型预测精度越高。判断标准, 当AUC值为0.5—0.6时, 模拟精度过低, 不可用进一步分析；AUC值为0.6—0.7时, 结果较差, 勉强可用；AUC值为0.7—0.8时, 结果一般, 可用；AUC值为0.8—0.9时, 结果良好；AUC值为0.9—1时, 结果优秀。基于上述标准, 对所模拟树种进行分析, 确定拟合效果。

按照文献^[27]方法, 提取12个树种各环境因子贡献率, 我们将贡献率之和超过60%的环境因子作为影响树种分布格局的主导环境因子。

1.3.2 潜在适生区分类

将MaxEnt模型模拟产生的树种分布图(*.asc文件)通过ArcGIS软件(ArcToolbox: Conversion Tools-To Raster-ASCII to Raster)转换为栅格数据(*.tif文件)并进行可视化处理。采用重分类(ArcToolbox: Spatial Analyst Tools-Reclass-reclassification)定义的间隔, 间隔大小为0.25, 将树种的分布区划分为非适生区(≤0.25)、低适生区(0.25—0.5)、中适生区(0.5—0.75)、高适生区(≥0.75)4个等级^[26], 并绘制各树种在东北地区潜在适生区分布预测图。

1.3.3 物种灭绝、损失及周转分析

将MaxEnt模型模拟预测当前及未来气候变化下获得的值为0(最低概率)—1(最高概率)的树种概率图导入ArcGIS软件中, 以敏感度和特异度阈值为标准, 利用con函数, 将树种概率图转换为0/1的二值图, 0代表不适宜分布区, 1代表适宜分布区^[5]。对于每个树种, 均生成当前分布二值图(Cmap)和未来分布二值图(Fmap)。为了区分气候变化情景引起的树种分布范围的变化, 使用公式Fmap×2+Cmap重新计算3种气候变化情景(SSP126、SSP245和SSP585)下的二值图, 得到树种分布的范围变化图(RCmap)；上述地图代码的整合过程, 使得RCmap栅格值为0、1、2和3, 可视化了每个气候情景下各树种的栖息地变化动态；0表示不适宜分布区, 1表示适宜分布区的丧失, 2表示适宜分布区的扩大(普遍迁移), 3表示稳定存在的适宜分布区。

为评估各树种的脆弱性, 在ArcGIS中统计RCmap中各树种的损失网格单元数(Nloss=sum(RCmap==1))、普遍迁移时潜在获得的网格单元数(Ngain=sum(RCmap==2))和稳定存在的网格单元数(Nstable=sum(RCmap==3)。在无迁移(各个树种不随环境因素及人为因素的变化而发生自然迁移)和在完全迁移(各个树种均随环境因素及人为因素的变化而发生适度自然迁移)两种假设下, 分别使用Nloss/(Nloss+Nstable)、(Ngain-Nloss)/(Nloss+Nstable)计算物种范围变化百分比(Prc)。范围变化的负值表示在两种迁移假设下树种栖息地的丧失。根据国际自然与资源保护联盟的标准, 将树种濒危等级划分为灭绝(-100%)、极危(-100% < Prc < -80%)、濒危(-80% < Prc < -50%)、易危(-50% < Prc < -30%)和无危(Prc>-30%)。这种方法虽然只考虑了未来气候变化的影响, 但它提供了气候变化对树种特有威胁的全面概述^{[5, 28]}。

在R语言中应用Raster包, 使用L/SR计算树种的损失率, 其中SR是当前气候条件下的物种丰富度层, L代表未来气候条件下单个网格内的树种损失总和。利用G/SR计算普遍迁移情况下的树木增益百分比, 其中G为未来气候条件下每个网格单元内获得的树种总量。利用(L+G)/(SR+G)公式计算了普遍迁移条件下网格单元的树种周转率^{[5, 28]}。

2 结果与分析 2.1 模型预测精度评价与主导因子筛选

依据MaxEnt模型预测东北地区建群树种分布格局, 样本数量及AUC值见表 3, 结果显示训练样本及验证样本AUC值在0.860(辽东桤木)—0.967(水曲柳、胡桃楸)之间, 均大于0.8, 说明MaxEnt模型预测精度良好, 结果可靠, 能够较好地模拟气候变化下东北地区建群树种的生态适宜区。

树种不同, 影响的主导环境因子有所不同(表 4)。兴安落叶松、樟子松的适宜性分布主要受到降水因素的影响。降水量季节性变化是影响二者的关键气候因子, 贡献率分别为36.4%、48.3%。红松、红皮云杉、山杨、春榆、白桦、黄檗、水曲柳、胡桃楸、蒙古栎的适宜性分布则受温度和降水因素的共同影响。其中, 红松、水曲柳、胡桃楸受年均温制约较大；山杨、白桦受降水量季节性变化因素影响较多, 尤其是白桦, 降水量季节性变化的贡献率高达38.7%；而春榆、黄檗适宜性分布受最暖季节降水量制约大, 均达到35%左右。辽东桤木的适宜性分布则由温度、降水、土地利用共同决定, 尤其是降水量季节性变化对其适宜性分布起关键作用, 贡献率达36.6%。

2.2 东北地区建群树种适宜性分布

预测结果表明(表 5、图 1), 当前12个树种均在长白山及小兴安岭地区有所分布, 白桦、蒙古栎、山杨、兴安落叶松、樟子松、辽东桤木等树种在大兴安岭地区也有分布。东北地区总面积为152万km², 兴安落叶松、樟子松、山杨、白桦、蒙古栎、辽东桤木适生区面积均能占到整个地区的40%左右。樟子松适生区面积最大, 占整个地区总面积的近45%；水曲柳、胡桃楸适生区面积较小, 占整个地区的15%左右。

东北地区建群树种的潜在适宜性分布会受到气候变化的影响, 不同气候情景对不同植物的影响存在差异(表 5、图 1)。主要表现在：红松潜在适生区面积变化幅度随气候变化情景强度的增强发生明显增加；兴安落叶松、山杨、春榆、白桦、水曲柳、胡桃楸、蒙古栎、辽东桤木潜在适生区面积变化幅度随气候变化情景强度的增强发生明显减少；红皮云杉、樟子松、黄檗潜在适生区面积变化幅度较小, 增加或下降的幅度在5%以内。

在SSP126和SSP245气候情景下, 各树种适生区变化较小。在SSP585气候情景下, 树种适生区发生明显变化, 红松潜在适生区面积占整个研究区域面积的42%左右, 其中中高适生区面积占总面积的近26%；樟子松潜在适生区总面积虽有小幅度下降, 但高适生区面积却大幅增加, 由当前的0.9%增加到11.8%。山杨、白桦、蒙古栎潜在适生区面积明显减少, 减少幅度超过20%。到21世纪末期, 兴安落叶松、春榆、白桦、水曲柳、胡桃楸、蒙古栎的适生区几乎消失殆尽。

2.3 东北地区建群树种受威胁状况

未来气候变化对东北地区关键乔木树种造成了威胁(图 2), 不同气候情景强度下, 其受威胁程度不同。且在无迁移和完全迁移两种假设下, 关键乔木树种面临的受威胁程度也存在差异。

SSP126气候情景下, 树种的受威胁程度相对较小, 无迁移假设下, 大约有37%的树种处于易危状态；完全迁移假设下, 约有12%的树种处于易危状态。在SSP245和SSP585两种气候情景下, 树种的受威胁程度整体表现为随年份的增长而增强的趋势。SSP245气候情景中, 在无迁移假设下, 到21世纪中期有58%的树种处于易危或濒危境地, 21世纪末将有8%的物种面临极危风险；在完全迁移假设下, 物种受威胁程度较小, 到21世纪末期, 有67%的物种处于低风险, 没有树种处于极危状态。SSP585气候情景中, 树种受威胁程度相对较高, 在无迁移假设下, 从21世纪中期8%树种处于极危状态, 到21世纪末面临极危风险的树种比例将达到40%；在完全迁移假设下, 到21世纪中后期有16%左右的树种处于易危状态。

整体看来, 在无迁移假设和SSP585气候情景下, 东北地区乔木树种的受威胁程度相对较高。在完全迁移假设及SSP126气候情景下, 树种大多处于低风险状态。

2.4 东北地区建群树种的损失和周转

在不同气候情景下, 东北地区乔木树种的损失、增益和周转情况存在一定差异(图 3、图 4)。在各个情景下, 树种增益变化不大, 与年份及气候情景间均无明显关系, 未来树种增益率高的区域集中分布在大兴安岭地区；树种损失和周转的变化范围(0—100%)较广, 整体看来, 树木损失和周转率随气候强度增强呈上升趋势。SSP126气候情景下, 损失率和周转率与年份间没有明显关系, 但SSP245和SSP585气候情景下, 东北地区树木损失率和周转率随年份的增长而增大。SSP126气候情景下, 2021—2100年东北地区树种平均损失率(31.5%左右)及周转率(49%)最低, 环松嫩平原损失及周转较大, 其中松嫩平原北部地区树木损失比例可超过80%。SSP245气候情景下, 2021—2100年东北地区树木平均损失率约38.9%, 平均周转率约59.7%。SSP585气候情景下, 2021—2100年东北地区树木平均损失率(约54.2%)和周转率(68.8%)最高。到21世纪中后期, 3种气候变化情景下东北地区树种平均损失率达60%以上, 平均周转率达80%以上, 其中大兴安岭地区、松嫩平原北部地区的树木损失比例可超过80%。

3 讨论

本研究对未来气候变化下东北地区乔木树种适宜生态区进行系统全面的预测分析, 为未来东北地区植物物种保护提供了新科学依据。应用MaxEnt模型结合ArcGIS软件模拟预测了东北地区12种建群树种未来气候变化下的生态适宜区。12个树种的MaxEnt模型预测精度均达到了0.8以上, 表明MaxEnt模型能够较为准确的预测东北地区乔木树种生态适宜区。未来气候情景适生区面积减小的关键种有：兴安落叶松、山杨、春榆、白桦、水曲柳、胡桃楸、蒙古栎、辽东桤木, 减小幅度达到10%—30%；适生区面积变化不大的关键种有：红皮云杉、樟子松、黄檗；适生区增加的关键种有：红松。可以看出不同树种具有迥异的适应能力。前人的研究也发现了类似的趋势, 如：张喜娟、于景华等^[29]应用MaxEnt模型对东北兴安落叶松林空间分布研究发现东北地区兴安落叶松林适生区衰减成为必然趋势。贾祥等^[30]对气候变化下阔叶红松林的4类主要建群树种红松、紫椴、水曲柳和蒙古栎潜在地理分布区的影响进行了分析, 结果显示：未来气候情景下, 到2080年, 大兴安岭的部分地区也将会成为阔叶红松林的高度适宜区。但整体看来阔叶红松林适生区面积有缩减的趋势, 而低度适宜区的面积有增加的趋势。

树种迁移与否对树种未来受威胁具有重要影响, 可以从物种损失、周转与增益率来解释。树种在一定区域内可能会发生适宜性扩散及迁移, 为反映这一情况, 我们做了两种极端迁移假设^[5]。发现树种迁移对树种损失可以起到缓解作用, 在一定程度上可降低树种的受威胁程度。2081—2100年, 在无迁移的假设下, 气候变化将威胁到60%的研究树种(3种SSPs气候情景的平均值)；在完全迁移的假设下, 36%的研究树种(3种SSPs气候情景的平均值)受到气候变化的威胁。这就表明, 物种适宜性迁移扩散的发生可能是树种在一定区域内适应气候变化影响的结果。这与Dyderski等^[31]假设欧洲森林树种发生迁移的范围有限, 发现大多数研究物种的适宜栖息地面积将显著减少的结论一致。未来, 东北地区树种增益率较高的区域主要发生在大兴安岭地区。气候变化会减缓区域气候从而引起区域陆地生态的改变^[32—33], 海拔及地形异质性对气候变暖具有缓冲作用^[34], 东北地区物种增益率高的地区大多为山地, 气候变化速度较慢, 更适合树种的生存。周转率代表的是物种组成的变化, 本研究中东北地区较高的周转率主要由物种损失所致。未来东北地区树种发生较高的周转情况, 在一定程度上可能会改变东北地区树木分布格局及森林结构, 从而可能对森林生态系统的功能产生重大影响^{[5, 31]}。未来需要开展物种自然迁移能力、人为辅助迁移等系统研究, 提高受威胁树种的适应能力。

本研究中树种分布数据, 是基于实测样方调查(科技部基础调查专项2000多个实地调查数据), 并结合各类标本、文献、图片等记录数据获得的, 相较于仅从各标本数据数据库获取物种分布记录, 数据来源更加全面可靠。但这也存在缺陷：尽管我们已经通过ArcGIS软件对数据点进行了筛选, 对25km²内的数据点只保留一个, 但也可能出现某一区域数据点过于密集的问题, 造成模拟结果的过拟合。以往的研究, 气候数据大都采用CMIP5中的典型浓度路径(RCP)中的情景模式, CMIP5采用的RCP排放路径均假设人为气溶胶排放在未来会大幅减少, 导致气溶胶情景的差异很小, 不能体现气溶胶影响近期气候变化的多种可能性^[35]。而本研究采用的(CMIP6)中的3种共享社会经济路径(SSPs)充分考虑了这一不足, 能够为区域气候预估以及减缓适应研究提供更加合理的模拟结果。SSPs还能更好地反映土地利用的不同变化对区域气候的影响^[36]。此外, 植物物种的生长、繁殖很大程度上会受到环境条件(生物因子和非生物因子)的制约。较大的环境条件变化会导致其在当前生境下难以生存, 进而影响其地理分布区域发生变化, 导致区域内物种数量增加或减少^[37—38]。

尽管本研究已经考虑了影响MaxEnt模型模拟效果的数据来源及环境变量, 但是也存在一定的局限性。目前很多研究使用多个模型的对比验证是物种分布最大限制因子和预测未来分布。如：Hernandez等^[39]基于四种生态位模型(Bioclim, Domain, GARP和MaxEnt), 对不同物种进行了不同尺度的模拟和评价, 研究发现, 对于地理范围小且分布数据极少的稀有物种(数据点≥5)而言, MaxEnt模型的预测能力最优。其次, MaxEnt模型较其他模型而言运行速度快、运算结果稳定、预测能力强^{[22, 40]}。本文仅仅使用了MaxEnt模型对数据进行分析, 开展多模型对比分析研究是未来研究的重要方向。其次, 本研究虽然我们已经考虑了气候、土地利用、地形等因素, 但未考虑土壤、历史因素、物种间相互作用等其他因素对植物物种分布格局的影响^[41—42], 这对于物种未来潜在分布的预测有一定的局限性, 结果可能存在较大差异。

4 结论

本研究表明MaxEnt模型能够准确预测东北地区建群树种的生态适宜区。结合当前和未来共享的社会经济路径(SSP126、SSP245、SSP585)的气候情景, 能够精确识别驱动东北地区建群树种当前和未来潜在生态适宜区的主导环境变量。降水量是影响东北地区建群树种的主导环境因子, 尤其降水量季节性变化对大部分乔木树种分布格局起关键作用。未来气候情景适生区面积减小的树种有：兴安落叶松、山杨、春榆、白桦、水曲柳、胡桃楸、蒙古栎、辽东桤木, 减小幅度达到10%—30%；适生区面积变化不大的树种有：红皮云杉、樟子松、黄檗；适生区增加的树种有：红松。未来东北地区建群树种的高周转率会对东北地区森林的可持续发展造成威胁, 但我们发现树木的适度迁移会减缓气候变化对树种损失的威胁。相关研究结果有助于提出更加科学有效的策略适应气候变化, 实现东北地区森林的可持续发展。