末次间冰期以来及未来气候情景下红皮云杉适生分布区模拟*

王 亚 王军辉 王福德 刘逸夫 谭灿灿 袁艳超 聂 稳刘建锋 常二梅 贾子瑞

（1.林木遗传育种国家重点实验室 中国林业科学研究院林业研究所 北京 100091；2.中国林业科学研究院森林生态环境与自然保护研究所 北京 100091；3.黑龙江省林业科学研究所 哈尔滨 150081）

植物与气候的相互关系一直是植物学和生态学研究的热点（安芷生等，2005），气候是决定物种分布的主要因素之一，物种分布格局变化能直观反映气候在不同时期的变化。基于物种分布模型（species distribution models, SDMs）和物种分布点的经纬度数据与环境因子，通过数学方法可得到各环境因子的贡献率，构建物种潜在分布范围，并以此推测未来和过去时期的适生分布区。物种分布预测在生态学、生物地理学等领域十分重要，已在土地利用变化、未来气候变化情景（Faleiroet al.，2013）中得到广泛应用，用来研究生物多样性损失（Bertrandet al.，2012）、濒危物种管理和保护（Guoet al.，2016）、生物热点区域确定（Walterset al., 2019）等。物种分布预测的可靠性受环境变量选取和模型参数的影响（Fourcadeet al.，2018;Sillero, 2011），环境变量之间的相关性会引入冗余信息并导致模拟结果偏差（Verbruggenet al.，2013）。因此，去除多余变量可消除变量间重线性关系，利于精确模拟物种的可能适生分布，并为生态保护等理论问题提供重要参考数据（Ashcroftet al.，2011）。

全球气候变化背景下的物种适生分布区不断变化。第四纪末次间冰期以来，全球气候频繁振荡，极大影响物种种群动态变化（张爱平等，2018）。物种经历盛冰期和间冰期的不断回旋，逐步形成当前的地理分布格局（邹旭等，2018）。物种分布模型研究表明，一些物种的分化和形成与历史事件相关。对濒危物种奇瓦瓦云杉（Picea chihuahuana）的最大熵模型分析发现，自末次盛冰期以来其适生区一直减少（Pinedo-Alvarezet al.，2019）；由于森林碎片化和气候变化，其孤立种群在未来面临消失的危险（Aguilar-Sotoet al.，2015）。对紫果云杉（Picea purpurea）、丽江云杉（Picea likiangensis）、青杄（Picea wilsonii）的潜在分布区预测发现，等温性、最暖季降水量和年均温是这3 个物种生态位分化的主要因子（张爱平等，2018），而且这3个物种在末次盛冰期发生扩张。对云杉属34 个物种的生态位分析表明，气候对云杉形态和物种分化的影响是由地理隔离、生态位保守性和生态适应等许多生物和非生物因素介导的，温度和降水分别是云杉系统发育和形态初级分化的主要驱动因素（Wanget al.，2017）。对东亚孑遗植物中全新世、末次盛冰期和未来的分布区预测结果发现，孑遗植物物种丰度和气候避难所相对应，中国西南部和越南北部是许多第三纪孑遗植物的长期避难所（Tanget al.，2018）。

红皮云杉（Picea koraiensis）又名臭松，为松科（Pinaceae）云杉属（Picea）植物，是中国东北重要的针叶树种和古老的孑遗树种，主要分布于大兴安岭北部、小兴安岭、长白山、完达山、张广才岭、老爷岭等地区，地理坐标范围在41°40'—53°15'N、122°—131°E 之间，垂直分布范围为250～1 100 m，属当地主要造林与用材树种。红皮云杉分布于针阔混交林中，并处于针叶林到阔叶林的过渡带，是研究植物分布响应气候变化的理想植物之一。红皮云杉具有涵养水源、防风固沙、水土保持等生态效益以及较高的经济价值，在木材生产和生态服务方面均有重要作用。目前，关于红皮云杉的研究主要涉及生态学系统特征和生物量（张含国等，2003）、遗传多样性（王秋玉等， 2004）、矿质元素含量（姜兴林等，2002）等，不同时期气候变化下的红皮云杉地理分布格局变化尚未见报道。鉴于此，本研究基于最大熵模型（Merowet al.，2013），结合GIS 技术，利用末次间冰期（last interglacial, LIG, 120～140 ka BP）、末次盛冰期（last glacial maximum, LGM,21 ka BP）、全新世中期（mid-holocene, MH, ca.6 000 years ago）、现代（current, 1960—1990 年）和未来2 个时期（2050s 和2070s）的气候数据和物种分布点数据，推测红皮云杉过去、现代和未来的适生分布区，识别主要气候影响因子，揭示气候对红皮云杉地理分布格局形成的影响，推测种群迁移的历史动态和生物避难所，并预测未来分布变化趋势，以期为该物种的资源管理和多样性保护提供科学依据。

1 材料与方法

1.1 物种分布数据收集

通过3 种方法收集红皮云杉地理分布数据：第一，实地调查收集红皮云杉的分布地点记录，每个调查点之间的距离不小于20 km，记录分布点的位置名称、经纬度和海拔；第二，检索查询相关文献并记录红皮云杉分布的地理信息；第三，搜索中国虚拟植物标本馆（http://www.cvh.ac.cn/）和全球生物多样性信息网络（https://www.gbif.org/）等数据库，收集红皮云杉的分布点数据。汇总分布点信息，删除地理信息不完整和重复的分布记录，最终确定59 个红皮云杉分布点。根据MaxEnt 模型的输入文件格式，将物种分布点数据按照物种名、分布点经度和纬度的顺序进行整理。

基于ENMTools v1.4 软件（Warrenet al.，2010），使用Trim duplicate occurrences 函数对整理好的物种分布数据进行冗余筛选，以避免地理聚集导致的空间自相关。最后共获得55 个红皮云杉有效分布点数据，用来构建模型。从自然资源部网站（http://www.mnr.gov.cn/）下载中国地图，使用ArcGIS v10.4.1 软件将红皮云杉分布点数据在地图上展示出来（图1）。

图1 当前红皮云杉在中国的分布记录Fig.1 Current distribution of P.koraiensis based on occurrence records

1.2 环境变量

云杉属植物生长发育和繁殖受多种气候因子影响，如温度、降水等。本研究选取与温度和降水相关的19 个环境变量进行最大熵模型分析,其中，现代19个与温度和降水有关的生物气候变量（Bio1—Bio19）来自WoldClim version 1.4（Ficket al.，2017），精度为2.5 arc-minutes（弧分）。为避免因高度相关的气候变量之间的相互影响而导致模拟结果过度拟合，利用MaxEnt v3.3.4 的刀切法（jackknife）将19 个气候变量按贡献率大小排序，筛选出贡献率大于1 的变量。利用ENMTools 1.4 在19 个气候变量之间进行相关系数分析，一对相关系数大于|0.9|的气候变量被视为非独立变量。对每对显著相关的变量，只选1 个有较大贡献的变量用于分析。

1.3 过去和未来气候数据

过去和未来时期的气候变量数据从Worldclim version 1.4（Ficket al.，2017）下载，其中LGM 和MH 时期的数据精度为2.5 arc-minutes（弧分），LIG 时期的数据精度为30 arc-seconds（弧秒）。大气环流模型选择CCSM4（community climate system model four），CCSM4 是能够模拟过去和未来气候变化对植物分布影响的最有效的气候模型之一（Abdelaalet al.，2019），已成功用于类似研究（Araret al.，2020; Liet al.，2020）。利用ArcGIS v10.4.1 将所有栅格文件转换成ASCII 格式，将LIG 数据重抽样为2.5 arc-minutes。下载未来2050s 和2070s 两个时间段的气候数据，选择RCP4.5（representative concentration pathway）和RCP8.5 两种排放气候场景，用于模拟红皮云杉未来地理分布格局。

1.4 最优模型选择

由于要素类型（feature class，FC）和正则化乘数（regularization multiplier，RM）2 个变量对模型结果影响较大，为获得更精准的适生区预测结果，本研究利用MaxEnt v3.3.4 软件、R 3.6.2 程序和kuenm 程序包（Coboset al.，2019）对这2 个变量进行优化。首先利用R 软件的kuenm 程序包将每个物种csv 格式的数据点分布文件拆分，70%的地理分布点用于训练集构建、30%用于测试集验证。基于拆分后的地理分布点数据和筛选后的环境变量，利用R 语言的kuenm 软件包从MaxEnt 模型的40 个正则化常数（0.1～4.0）和10种要素类型（L、Q、P、H 的2 种组合和3 种组合）中选择最优模型组合，要素类型含义如下：L（linear features）为线性，Q（quadratic features）为二次型，P（product features）为乘积型，H（hinge features）为片段化。利用kuenm 程序包逐一检验上述400 组参数，根据最佳模型选择标准（Coboset al.，2019），选择统计学上显著、遗漏率（omission rate）低于阈值（0.05）且ΔAICc（delta Akaike information criterion）值不高于2 的模型，具有ΔAICc 最小值的模型被作为最优模型。

1.5 红皮云杉潜在分布预测

应用MaxEnt 最大熵模型预测红皮云杉在过去、当前和未来气候情景下的潜在分布。将选定的物种分布点数据和环境变量导入MaxEnt v3.3.4 软件。采用10 倍交叉验证法对模型进行训练和验证，输出格式为Logistic，输出文件类型为asc，最大迭代次数设置为5 000，线程设置为10，重复次数为10 次。MaxEnt建模完成后，采用受试者工作特征曲线（receiver operating characteristic，ROC）评估模型，并使用曲线下面积（area under the ROC curve，AUC）测量模型可信度（Alloucheet al.，2006），AUC＞0.8 时表明物种适生区预测精度较好，AUC＞0.9 时表明预测精度非常好。本研究中，使用ArcGIS v10.4.1 软件创建当前、过去（LIG、LGM 和MH）和未来（2050s 和2070s）的适生区预测图。MaxEnt 生成结果导入ArcGIS v10.4.1 软件，并以中国地图为底层进行分析。使用自然断点法划分红皮云杉的适生区等级：非适生区（0～0.1）、低适生区（0.1～0.3）、中适生区（0.3～0.5）和高适生区（0.5～1），计算各时期的适生面积。

1.6 适生分布区质心转移

为明确红皮云杉适生分布区的变化趋势，在ArcGIS v10.4.1 软件中加载SDMToolbox （Brown，2014）工具箱，利用Centroid changes command 命令，计算不同时期的分布质心。将这些质心在ArcGIS v10.4.1 软件中利用矢量箭头投影，表示适生区变化的迁移距离和方向。

2 结果与分析

2.1 模型选择和精度评价

为精准预测适生区，本研究筛选出6 个环境变量，并选择最优的要素类型和正则化乘数，用于最大熵模型分析。6 个环境变量为年均气温（Bio1）、昼夜温差与年温差比值（Bio3）、气温季节性变动系数（Bio4）、降水量季节性变化（Bio15）、最干季度降水量（Bio17）、最暖季度降水量（Bio18），最优模型要素类型为PH（product features, hinge features），最优正则化乘数为2.3。10 次重复训练数据集和测试数据集的平均AUC 分别为为0.995 8 和0.995 7，均接近于1，说明模型有较高稳定性和精度，模拟结果可信度高（图2）。

图2 基于MaxEnt 模型预测红皮云杉分布的受试者工作特征曲线Fig.2 The ROC curve predicting the distribution of P.koraiensis based on MaxEnt

2.2 环境变量筛选

根据MaxEnt 模型的输出环境因子贡献率和刀切法分析结果，比较各环境变量对红皮云杉潜在分布区变化的贡献率（表1、图3），最暖季度降水量（Bio18）、气温季节性变动系数（Bio4）和降水量季节性变化（Bio15）为3 个主要变量，贡献率分别达40.6%、28.9%和24.6%，年均气温（Bio1）、最干季度降水量（Bio17）和等温性（昼夜温差与年温差比值）（Bio3）的贡献率分别为4%、1.5%和0.4%。

表1 影响红皮云杉适生分布的环境变量贡献百分率Tab.1 Contribution percentage of environmental variables affecting the suitable distribution of P.koraiensis

图3 环境变量刀切法检验Fig.3 Jackknife tests for environmental variables

重要变量响应曲线（图4）显示出对红皮云杉分布影响较大的变量（贡献率＞2）适生范围：最暖季度降水量（Bio18）为332.9～527.9 mm，气温季节性变动系数（Bio4）为129.4～173.5，降水量季节性变化（Bio15）为87.6～111.2，年均气温（Bio1）为-1.4～5.6 ℃。

图4 红皮云杉对重要环境变量的响应曲线Fig.4 Response curves of P.koraiensis to important environmental variables

2.3 目前、过去和未来适生分布区

红皮云杉当前适生分布区包括中国的大兴安岭、小兴安岭、长白山脉及其支脉等地区，涉及黑龙江、吉林、辽宁等省，其中高适生区的分布范围与实际分布相符（图5）。红皮云杉当前的总适生区面积达133.58×104km2，其中，低、中、高适生分布区面积分别为40.93×104、42.36×104和50.29×104km2（表2）。

表2 红皮云杉在多种气候条件下的适生面积①Tab.2 Prediction of suitable area of P.koraiensis under different climate scenarios

图5 当代、过去和未来气候情景下红皮云杉在中国适生分布区Fig.5 Suitable distribution areas of P.koraiensis in China under contemporary, past and future climate conditions

历史气候条件下，红皮云杉在末次间冰期（LIG）、末次盛冰期（LGM）和中全新世（MH）的适生分布区面积分别为目前的102.59%，105.20% 和100.53%；高适生区和总适生区的变化趋势一致，末次盛冰期种群扩张，全新世中期及以后种群持续收缩（图5、表2）。

未来气候情景下，红皮云杉的总适生区和高适生区面积均呈减少趋势，高适生区边缘呈片段化分布，而小兴安岭腹地、长白山腹地、老爷岭以及完达山等地的种群较为稳定（图5、表2）。红皮云杉在未来气候情景中2 种温室气体排放模式下（RCP4.5-2050s、RCP4.5-2070s、RCP8.5-2050s 和RCP8.5-2070s）的总适生区面积分别为当前的99.06%、96.88%、98.05%和95.08%。在高排放（RCP8.5）模式下，红皮云杉适生区面积缩减更大。RCP4.5 模型预测2050s 时期高适生区缩减的区域在大兴安岭北部边缘（塔河县）到中段（嫩江市）、长白山西部（辉南县和抚松县）；2070s 时期这2 个区域持续缩减，大兴安岭北部（蒙克山地区）和中部（高峰地区）、张广才岭北部（勃利县和桦南县）和东部（海林市）的高适生区消失。RCP8.5 模型预测2050s 时期高适生区缩减的区域在张广才岭北部（勃利县和依兰县）和东部（海林市）；2070s 时期，红皮云杉高适生区边缘持续碎片化，长白山西部（辉南县和抚松县）的高适生区消失。

2.4 适生分布区质心转移

通过重建现代、过去和未来不同时期气候条件下红皮云杉的适生区质心坐标，分析质心变化。红皮云杉LIG 时期的分布质心坐标在46°14'27.38''N,127°1'46.92''E；LGM 时期分布质心（46°10'54.01''N,127°14'49.34''E）向东南转移21.62 km；MH 时期分布质心（46°15'0.44''N, 127°18'54.83''E）向东北转移8.96 km；当前时期分布质心（46°6'14.27''N, 126°14'28.81''E）向西南转移104.43 km。未来气候条件下，当气候情景为中排放浓度（RCP4.5）时，2050s年代分布质心（46°25'2.47''N,127°17'54.03''E）向东北转移105.72 km；2070s 年代分布质心（46°15'33.81''N, 127°3'15.60''E）向西南转移27.4 km；当气候情景为高排放浓度（RCP8.5）时，2050s年代分布质心（45°49'27.43''N, 126°57'0.64''E）向东北转移105.72 km；2070s 年代分布质心（46°11'31.04''N,127°1'25.55''E）向西南转移27.4 km（表3、图6）。

表3 不同时期红皮云杉适生区质心坐标和迁徙距离①Tab.3 Centroid coordinates and migration distance of suitable growth areas for P.koraiensis in different periods

图6 不同时期红皮云杉适生区质心的迁移路线Fig.6 Migration route of centroid in the suitable growth area of P.koraiensis in different periods

3 讨论

3.1 影响红皮云杉分布的环境因子

根据自然选择学说，物种形成是由于环境不同而引起的适应性改变。环境是物种形成中最重要的因素，而全球气候变化在环境变化中占有重要地位。本研究采用19 个环境变量，计算Pearson 相关系数，排除高相关性变量中贡献率较低的变量，并通过AUC评价建模精度。本研究训练集和测试集的AUC 均高于0.99，说明模拟结果具有较高可靠性，当前适生分布区的预测更接近实际情况。

气候是决定物种分布变化的关键因素，其中温度因子通过影响植物光合作用、呼吸作用来控制物种分布；而植物生长发育也受降水条件制约（Abeliet al.，2020）。本研究采用MaxEnt 自带刀切法，不仅可提高模型预测可信度（De Cauweret al.，2014），而且能客观、精确地评价19 个气候变量。最暖季度降水量（Bio18）、气温季节性变动系数（Bio4）、降水量季节性变化（Bio15）对红皮云杉的分布影响最大，其次为年均气温（Bio1）、最干季度降水量（Bio17）、昼夜温差与年温差比值（Bio3）等环境因子。其中，与水分相关的因子（Bio18、Bio15、Bio17）贡献率为66.7%，说明与温度因子相比，红皮云杉适生区分布受降水影响更大，与王改妮（2018）、段国珍（2019）对红皮云杉及其近缘种分布格局的模拟结果相似。这是因为红皮云杉在云杉属物种分布区的最北端，北方河流冬季有封冻期，降水量较少，对于耐阴、喜湿润的云杉属植物来说，生长季降水量不足会限制生长，长期缺水则会导致植物的适应性分化。与其近缘种白杄（Picea meyeri）和沙地云杉（Picea mongolica）相比，红皮云杉适应湿冷环境的习性十分明显，无论是大气湿度还是土壤水分供给状况均高于前二者。与云杉组（Sect.Picea）内其他类群相比，红皮云杉是对水分需求较高的种群，其分布区内的年均降水量仅次于云杉（Picea asperata）和青扦（Picea wilsonii）；与同域种鱼鳞云杉（Picea jezoensis）（云杉属鱼鳞云杉组sect.Omorica）相比，红皮云杉更适应低湿的土壤环境，在排水良好的山坡，红皮云杉的竞争力不及鱼鳞云杉，分布数量较少。此外，影响擦木（Sassafras tzumu）、柳杉（Cryptomeria japonica）等植物适生区分布的主要限制因子也是降水（刘想等，2018；戴旻峻等，2022），这些植物稳定生长在山地沟谷和水分肥沃的地区，同样，红皮云杉主要生长在河谷溪旁、河岸阶地和山麓。综上表明，红皮云杉对水分的需求较高，降水是影响红皮云杉分布和扩散潜力的主要限制因子，因此红皮云杉在推广种植或异地保护时，可选择在降水充沛地区。

温度在影响红皮云杉分布上也起到较重要作用。根据刀切法结果，年均气温（Bio1）对红皮云杉分布的影响不可替代，Bio1 的适宜范围在-1.4～5.6 ℃，最适值为4.3 ℃，在-1.4～4.3 ℃时，分布概率随年均温升高而增大，在4.3～5.6 ℃时，分布概率随年均温升高而降低，该因子可能限制红皮云杉继续向南扩张。气温季节性变动系数（Bio4）的贡献率达28.9%，Bio4 的适宜范围在129.4～173.5，最适值为155.6。Bio4 代表温度在一年中的变化幅度，数值越大温度变幅越大，因此均一稳定的生长环境可能对促进红皮云杉种群扩张有积极作用。

3.2 红皮云杉在不同气候情景下的分布格局变化和质心转移

植物在与气候的长期交互作用下形成一定的地理分布格局。通常，拥有更广生态位的物种具有更强的环境适应能力，可充分利用其周边的自然资源（Yinget al.，2016），因此在预测红皮云杉林分更新状况时必须考虑其对气候变化的适应情况。本研究采用数据库数据与实地调查数据相结合的方法，利用古气候数据、现代气候数据和未来气候变化情景，构建红皮云杉物种分布模型，并重建LIG、LGM 和MH 3个历史时期的地理分布格局。

研究发现，红皮云杉种群适宜生境曾经历多次适应性转换，包括末次盛冰期的种群扩张以及全新世中期及之后的种群持续收缩。末次间冰期的陆地平均气温比现代高1～2 ℃（何韬，2019），与现代更相似；末次盛冰期的陆地平均气温比现代低5～10 ℃，而中高纬度地区低20 ℃以上（何韬，2019）。本研究结果和Stewart 等（2010）提出的理论是吻合的，即适应寒冷的植物在间冰期的分布范围比盛冰期小，如白云杉（Picea glauca）（Stewartet al.，2010 ；Bannisteret al.，2001），以往对云杉属其他物种的生态位模拟分析中也有一致的结果（Shenet al.，2019；王改妮，2018；段国珍，2019）。相对于末次间冰期，末次盛冰期气候较冷，红皮云杉的质心向东南方向转移；全新世中期气候再次回暖，红皮云杉的质心向东北方向转移。红皮云杉的质心转移与气候变暖条件下物种向温度更低的区域迁移的规律是相符的。

温度适当升高可以加速植物生长发育，并影响种群数量变化。随着全球温度持续上升，在未来气候变化2 种温室气体排放情景下，红皮云杉的适生区面积均呈现持续减少趋势，高度适生区呈片段化分布，全球气候变暖不利于红皮云杉适生区扩张。对附子（Aconitum carmichaeli）（韦洁敏，2020） 、 金钱松（Pseudolarix amabilis）（张文秀，2020）、桫椤（Alsophila spinulosa）（杨启杰等，2021）等植物的潜在分布预测分析也发现，随着未来全球气候变暖，这些物种适生区面积会持续减少，与本研究结论一致。未来气候变化情景中，与RCP4.5 排放路径相比，RCP8.5 排放路径下的适生区面积缩减速度较大，红皮云杉在高排放情景下适生区面积变化更显著，对气候变化响应明显。随着未来气候变暖，红皮云杉总适生区质心先向高纬度迁移，后逐渐向低纬度回落，出现这种情况的主要原因是红皮云杉潜在适生区在向西南方向海拔更高的长白山地区扩张。全球变暖影响植物的多样性、分布和迁移模式，大多数物种迁移到高海拔和高纬度地区（Parmesan，2006）。木梨（Pyrus xerophila）、四合木（Tetraena mongolica）未来的潜在分布范围有向西北方向的高纬度地区迁移的趋势（刘超等，2020；段义忠等，2019）。未来气候情景下，红皮云杉的高适生区缩减，边缘呈片段化分布，缩减区域包括大兴安岭的塔河、嫩江、蒙克山和高峰，长白山西部的辉南和抚松，张广才岭的勃利、桦南、依兰和海林等地，建议对这些区域的种群加强就地保护、建立自然保护区、减少人为破坏与干扰。小兴安岭腹地和长白山腹地的种群生境较为稳定，是建立红皮云杉异地保护的最佳候选区域。

3.3 红皮云杉避难所推测

当全球经历第四纪冰期的漫长考验时，众多生物依靠生物避难所的保护而得以幸存，并在冰期后重新扩散。识别生物避难所有3 个主要途径：一是基于生物地理模式，通过收集目标生物体的古生物学、生态学和遗传学数据，量化生物避难所的时空纬度；二是调查物种避难所栖息地的产生过程；三是通过模拟物种分布格局对环境变化的响应来识别避难所（胡忠俊等，2013）。冰期避难所理论表明，避难所在环境波动期间的生态条件相对稳定，可促进物种多样性的积累。因此物种具有较高水平遗传多样性的地理分布区域，可以作为寻找避难所的依据之一（石甜，2021）。

红皮云杉形成于第三纪中晚期（4.12～7.64 Ma, 百万年）(Shaoet al.，2019)，说明最先引起红皮云杉群体遗传分化的因素是第三纪的地质和气候事件。在第三纪中晚期，青藏高原出现快速抬升（8～10 Ma），加速了亚洲大陆内部的干燥度和亚洲季风的产生（Anet al.，2001）。第三纪末发生地壳运动，使得大兴安岭沿东侧断层翘起；第四纪初期的构造运动又引起火山喷发，在小兴安岭地区形成大片玄武岩熔岩台地，同时有轻微翘起，沿北西方向断裂隆起（樊祺诚等，2015；刘建强等，2017）；第四纪中晚期至全新世，长白山发生多次火山喷发，长白山主体于全新世形成后，地壳相对稳定（王团华，2006）。红皮云杉的分化时间段发生于第三纪末至第四纪期间，东亚季风气候增强和地壳变动的时间基本一致。与其分化时间接近的还有第三纪植物枫香（Liquidambar formosana）（10.30 Ma）（孙荣喜，2017）、水青树（Tetracentron sinense）（9.6 Ma）（Sunet al.，2014）、青冈栋属（Cyclobalanopsis）（9.07 Ma）（Denget al.，2018）等。

植物在冰期会退缩到避难所，耐寒植物在冰期往往会扩张且在间冰期发生种群收缩，故生物避难所还应包括间冰期避难所（Stewartet al.，2010）。红皮云杉属于耐寒植物，通过分子证据、适生分布区模拟结果以及孢粉化石证据可推测红皮云杉的间冰期避难所。分子证据显示，在红皮云杉不同地理种群中，长白山腹地的临江种群（π= 0.002 35）和天桥岭种群（π =0.002 32）均具有较高的核苷酸多样性（nucleotide diversity, π）(Wanget al.，2023)；此外，由于多个历史时期的物种分布重叠区往往是植物潜在避难所（Carnavalet al.，2008），结合对末次间冰期和末次盛冰期的适生分布区模拟结果推测，长白山腹地可能是红皮云杉第四纪间冰期避难所，且在该种群进化过程中发挥了重要作用。从生境特征上看，长白山地区属于红皮云杉分布区内海拔最高的区域（海拔100 0～274 9 m），较大的海拔落差为红皮云杉等山地植物适应气候变化提供了充足的缓冲环境，且该区域末次间冰期时的气候湿冷（孙广友等，2012），为红皮云杉提供了生长所需的充沛水分。此外，小兴安岭腹地也可能是红皮云杉的避难所之一。分子证据显示小兴安岭腹地的种群（沾河）核苷酸多样性较高（π= 0.002 31）（Wanget al.，2023）；此地又是多个时期适生分布区的重叠区域；从孢粉数据来看，松辽盆地发现大量中更新世晚期的云杉属孢粉化石（石宁，1996），与末次间冰期的时间吻合；小兴安岭处于华力西褶皱带，第三纪时山体北部与东北方的结雅河低地、西南方的松嫩平原联成一片，在上新世末更新世初沿北东向和北西向断裂抬升为山岭，与结雅河低地及松嫩平原分离（王欢等，2018），其复杂多变的山地峡谷地貌为红皮云杉提供了丰富的避难场所。因此，结合分子证据、适生区模拟、孢粉化石和地形地貌，推测红皮云杉在间冰期的生物避难所在小兴安岭腹地和长白山腹地。

由于气候数据是根据数学方法估计所得，与古代实际气候存在差异，且未考虑非生物因素的影响，因而MaxEnt 模型在适生分布区预测方面存在一定局限性。然而，通过与其他物种分布模式的比较，MaxEnt具有较高的可靠性和稳定性（Phillipset al.，2006；Elithet al.，2010；朱耿平等，2016）。在模拟物种当代适生分布区时，其结果和物种实际分布范围具有高度一致性（高文强等，2016；赵光华等，2021；王东升等，2022）。该方法可以呈现物种在不同时期的地理分布格局，判断环境因素对物种分化和形成的影响，推测物种的生物避难所，预测物种在未来气候下的分布范围变化趋势并为物种保护提供依据（李国庆等，2013；朱耿平等，2013）。

4 结论

MaxEnt 模型预测红皮云杉适生分布区结果可靠，训练集和测试集的AUC 均大于0.99。与温度因子相比，降水对红皮云杉的地理分布影响更重要，在其推广种植或者异地保护时，可选择降水充沛地区。温度因子可能限制红皮云杉向南扩张，均一稳定的生长环境对促进其种群扩张有积极作用。由于第三纪末至第四纪期间的地质和气候事件，红皮云杉的适生区经历多次转换，包括末次盛冰期种群扩张、全新世中期及之后的种群持续收缩。红皮云杉的适生分布区质心转移与气候变暖条件下物种向温度更低的区域迁移的规律相符。全球气候变暖不利于红皮云杉的种群扩张，未来情景下红皮云杉的适生区呈缩减趋势。应考虑优先保护大兴安岭北部到中部、长白山西部和张广才岭北部等高适生区的边缘种群。根据分子证据、适生区模拟，结合孢粉化石和地形地貌，推测小兴安岭腹地和长白山腹地是红皮云杉的间冰期避难所。