气候变化背景下我国长柄扁桃潜在适生区预测*

张淑宁 陈俊兴 敖 敦 红 梅 张雅茜 包福海 王 淋乌云塔娜 白玉娥 包文泉

（1.内蒙古农业大学 呼和浩特 010018；2.中国林业科学研究院经济林研究所 郑州 450003）

植被分布格局是植物对包括气候在内的环境条件长期适应的结果（Fanget al.，2006）。据IPCC 第五次评估报告，过去100 年全球气温上升近1 ℃，预测到2100 年仍将升温1.4～5.8 ℃，这将对陆地生态系统产生巨大影响（Stockeret al.，2014）。研究表明，全球气候变化对植物的物候、生长和生殖产生显著影响，进而引起植物分布范围的迁移，使植物分布由低纬度和低海拔向高纬度和高海拔迁移，或造成其分布范围缩小（Thomaset al.，2004；Matiaset al.，2017）。尤其是濒危植物，因其遗传多样性较低，对气候变化更敏感，若不采取有力保护措施，将大规模灭绝，全球生态系统会遭到严重破坏（Zhuet al.，2012）。研究濒危植物与气候因子的关系和预测其适生区受气候变化的影响，可为其保护利用提供重要理论依据。

生态位模型可基于物种分布数据及其对应的环境参数来获取物种的生态需求，进而预测适生区域（Schorret al.，2012）。近年来，生态位模型在植物适生区预测、植物保护和入侵植物管理以及探究气候变化对物种分布格局的影响方面得到了广泛应用（Elithet al.，2009） 。 目前常用的生态位模型有MaxEnt、CLIMEX、GLM、HABITAT、GARP 和BIOCLIM 模型，其中，MaxEnt 模型预测准确、可靠，稳定性强（Elithet al.，2009），已在互花米草（Spartina alterniflora）（张静涵等，2020）和朝鲜崖柏（Thuja koraiensis）（兰雪涵等，2021）等众多物种中应用，是目前预测物种分布区及其变迁的首选模型。

长柄扁桃（Amygdalus pedunculata）属蔷薇科（Rosaceae）落叶灌木，自然分布于内蒙古、河北、陕西等干旱半干旱区的沙地或山地，是我国特有的木本油料树种（苏贵兴，1987）。其核仁营养丰富，富含油酸、亚油酸和蛋白质等成分，是食品、化妆品和药品等的重要原料，具有极高的开发利用价值（郭春会等，2005）。但气候变化、基础建设及乱砍滥伐等使长柄扁桃生境严重破坏，分布范围和数量逐年减少，现已被列入我国重点保护的濒危植物，长柄扁桃野生资源的保护与修复工作已刻不容缓。然而，目前对长柄扁桃的研究主要集中于人工栽植、成分提取、生理生化、遗传多样性评价和良种选育等方面（慕晓炜等，2014；孙树臣等，2020）；而关于长柄扁桃适生区预测及受气候变化影响的研究较少。本文将优化MaxEnt 模型的参数，确定长柄扁桃当前适生区域，推测末次间冰期、末次盛冰期和全新世中期3 个历史时期以及未来气候情景SSP126、SSP245 和SSP585 下的长柄扁桃适生区变化趋势，以期为长柄扁桃资源的保护与人工栽植提供科学依据。

1 研究方法

1.1 地理分布数据的收集

在2015—2019 年，对长柄扁桃自然分布区进行了系统调查，从内蒙古、陕西、河北、山西等地共获取61 个分布数据；结合相关研究文献、国家标本馆、全球生物多样性信息网络数据库，获取长柄扁桃187 个分布数据，对其删除错误及重复分布点，最终获得148 个分布点（图1）。

图1 我国长柄扁桃地理分布点信息Fig.1 Distribution information of A.pedunculata in China

1.2 环境因子的获取及筛选

通过分析长柄扁桃分布区的气候、地形和土壤因子，初步筛选出34 个环境因子作为建模变量，包括19个生物气候因子、14 个土壤因子和1 个地形因子（海拔）；其中，气候因子数据取自世界气候数据库（www.worldclim.org）的末次间冰期（last inter glacial）、末次盛冰期（last glacial maximum）、全新世中期（mid holocene）、当前（1960—1990 年）、2050s（2041—2060年）和2090s（2081—2100 年）6 个时期，包括最暖季平均气温、月平均昼夜气温差、最热月最高气温、最湿月降水量、年均气温、最湿季平均气温、最干季平均气温、气温变异系数、气温年较差、年降水量、最暖季降水量、等温性、最冷季平均气温、最干季降水量、最冷季降水量、最冷月最低气温、降水量变异系数、最湿季降水量和最干月降水量。土壤因子和地形因子数据来源于世界土壤数据库的中国土壤数据集（http://www.fao.org/faostat/en/#data.）；土壤因子包括0～30 cm 土层中的砾石体积百分比、阳离子交换量、土壤质地类型、砂砾百分比、黏粒百分比、粉砂百分比、黏性成分阳离子交换量、土壤密度、pH 值、交换性盐基总和、盐基饱和度、可交换钠盐百分比、电导率和有机碳含量。

将气候因子和海拔因子数据按1∶100 万中国行政地图矢量图进行剪裁，并利用ArcGis 将其转换为ASC 格式；运用ArcGis 提取MU_GLOBAL 图层的表层土壤因子栅格图层，将其转为ASC 格式；最后将各环境因子的空间分辨率统一设置为30″，使环境因子数据栅格图像范围和像元大小完全一致。

本研究涉及的地图数据下载于自然资源部官网（http://www.mnr.gov.cn/）。以社会经济共享路径（shared socio-economic pathways, SSP）的低浓度（SSP126）、中浓度（SSP245）和高浓度（SSP585）排放情景代表未来气候系统模式，假设未来50 年土壤和地形因子不发生明显变化（曹丽格等，2012；翁宇威等，2020），对长柄扁桃进行适生区预测。为避免各环境因子间的共线性而导致模型预测结果过度拟合，本研究利用方差膨胀因子VIF 和Pearson 相关性检验，筛选相关性小于0.8 且VIF 值小于10 的因子（张天蛟等，2017）；同时，对相关系数大于0.8 的因子保留其生态学意义更重要的因子（Yeet al.，2020）。最终，筛选出8 个气候因子、4 个土壤因子和1 个地形因子（表1）。

表1 预测长柄扁桃地理分布的环境因子Tab.1 The environmental variables for prediction of geographical distribution of A.pedunculata

1.3 物种分布模型的构建与评价

将长柄扁桃148 个分布数据和所筛选的13 个环境因子数据导入MaxEnt 模型，随机选取75%的样点作为训练数据集构建模型，将剩余25%的样点作为测试数据集验证模型，设置10 次重复；并用Jackknife 评估各环境因子的权重，结合其贡献率确定主导环境因子。

为准确评价模型的可靠性，采用刀切法分析环境因子的贡献率，利用受试者工作特征曲线ROC 下的面积AUC 值评估所建模型精度（Fandet al.，2014）。通常认为AUC 值越大预测越精确，其中，AUC 值在0.5～0.7 时表示预测较差，0.8～0.9 时表示较好，而0.9～1.0 时则表示非常好（Freitaset al.，2019）。

1.4 模型的优化

应用ENMeval 程序包优化MaxEnt 模型，将调控倍频设置为0.5～8.0，每次间隔0.5，共设置16 种调控倍频；并采用9 个特征组合，即：L、LQ、H、LQH、LQHP、LQHPT、QHP、QHPT 和HPT。ENMeval 程序包将上述144 种参数组合进行测试，最终依据Akaike信息量准则（AICc），以delta.AICc 作为测试模型的复杂度和拟合度，选取参数组合，其中delta.AICc 值越低，表示模型预测越准确（李安等，2020；赵光华等，2021）。

1.5 数据处理

根据MaxEnt 模型预测的长柄扁桃适宜性阈值，结合自然断点法划分长柄扁桃生境适宜性指数，并运用ArcGis 软件将其可视化。将长柄扁桃生境适宜性划分为4 个等级，即不适生区（0～0.1）、较不适生区（0.1～0.3）、一般适生区（0.3～0.5）和高度适生区（0.5～1.0）。比较不同时期的适生区变化，获得气候变化背景下长柄扁桃未来空间分布格局变化；调用SDMToolbox工具，计算长柄扁桃当前和未来气候情景下的适生区质心位置，使用dists 函数计算其质心迁移距离，以质心位置的变化表示其适生区空间分布的迁移方向。用ArcGis 将当前和未来适生区数据转换为栅格数据并对其进行模糊叠加，得到长柄扁桃潜在地理分布图层，再进行重新划分适生等级，获取长柄扁桃适生区划图。

2 结果与分析

2.1 MaxEnt 模型的优化与评价

将148 个分布点和13 个环境因子导入MaxEnt 模型，预测长柄扁桃的潜在分布区。在MaxEnt 默认参数下（RM=1，FC=LQHPT），delta.AICc 值为225.87；而当RM=1.5、FC=LQHPT 时，delta.AICc 值为0；依据AIC 信息准则，此时模型最优（图2），模型的10 次预测训练集和测试集AUC 值也均大于0.90，均值为0.967，其训练遗漏率较低（图3），说明模型预测较好，结果可靠。

图2 不同参数设置下MaxEnt 模型的评估结果Fig.2 Evaluation results of MaxEnt model under different settings

图3 最优参数设置下MaxEnt 模型ROC 曲线和遗漏率曲线Fig.3 ROC curve and Omission curve of MaxEnt model under optimal parameter setting

2.2 长柄扁桃在我国当前潜在地理分布

基于MaxEnt 模型，构建长柄扁桃当前适生分布区（图4），其当前总适生区（一般适生区和高度适生区）面积为447 274.3 hm2，主要分布于内蒙古高原和黄土高原地区，其中高度适生区为219 166.7 hm2，占我国国土总面积的2.3%（表2）。当前我国长柄扁桃高度适生区集中于内蒙古的“呼包鄂”地区，此外，赤峰、乌兰察布和陕西榆林也少有分布。一般适生区则以高度适生区为中心向四周分布，从西至东依次分布于甘肃、宁夏、内蒙古、陕西、山西和河北等省区。

表2 长柄扁桃不同时期的适生区面积①Tab.2 The suitable living area for A.pedunculata under different climates scenarios/years

图4 基于MaxEnt 模型预测的我国长柄扁桃当前适生区分布Fig.4 Distribution of current suitable areas of A.pedunculata in China based on MaxEnt

2.3 长柄扁桃在我国历史潜在地理分布

通过对我国长柄扁桃历史潜在地理分布的推测发现，在末次间冰期，长柄扁桃适生分布区较广，西可达西藏，东可达山东，该时期长柄扁桃的高度适生区面积为620 431.9 hm2，为当前高度适生区的2.83 倍（图5，表2）；而在末次盛冰期，长柄扁桃的适生区严重收缩，空间丧失率高达81.44%（图6）；并且这种收缩趋势延续至全新世中期，此时高度适生区面积仅为末次盛冰期的44.8%（图6，表3）；但全新世中期有内蒙古东部的通辽，山西太原以北，以及北京、河北等较高纬度地区已成为新增适生区。因此，我国长柄扁桃从末次间冰期至末次盛冰期，其空间分布格局变化显著，主要由高纬度和低纬度向中高纬度收缩；而从末次盛冰期至全新世中期，适生区则由东向西进一步收缩（图5）。

表3 长柄扁桃不同时期适生区的变化①Tab.3 Change of A.pedunculata suitable area in different periods

图5 历史、当前及未来不同气候情景下长柄扁桃在中国的适生区分布Fig.5 Distribution of suitable areas of A.pedunculata under historical, current and future climate scenarios in China SSP：社会经济共享路径Shared socio-economic pathways.

图6 历史、当前及未来不同气候情景下长柄扁桃在中国的空间分布格局变化Fig.6 The spatial pattern of A.pedunculata in China under historical, current and future climate scenarios SSP：社会经济共享路径Shared socio-economic pathways.

2.4 未来气候情景下的长柄扁桃适生区预测

在未来不同气候情景下长柄扁桃适生区除在2050s-SSP126 情景下较当前有增加外（增加42.18%），其余未来气候情景下均有不同程度的缩减（图6）。在2050s-SSP126 情景下，适生区丧失面积最少（92 100 hm2），新增率最大（63.32%），新增区可分为东、西部两区，东部区为内蒙古锡林郭勒、赤峰、通辽以及内蒙古与吉林省交界区域，而西部区为甘肃、宁夏与内蒙古交界带（表3）。该情景下，长柄扁桃适生区面积约为635 954.8 hm2，是当前适生区的1.42 倍，说明此时将有大面积非适生区转变为适生区。

2090s-SSP245 情景下，长柄扁桃适生区面积约为345 833.3 hm2，较当前适生区收缩了近1/3，适生区新增率最小（4.33%），新增区域分布于天津、河北、山西等地。在2090s-SSP126 情景下，长柄扁桃适生区丧失面积最大（241 200 hm2），丧失率高达53.83%，而新增率仅为9.15%；此时，吉林、辽宁、通辽及赤峰地区的适生区将大面积丧失，仅有内蒙古锡林郭勒以及陕西延安的部分地区新增为适生区。

通过预测未来适生区变化可知，长柄扁桃对不同气候情景下气候变化的响应存在一定差异（图5）。在SSP126 情景下，高度适生区面积具有先增后减的趋势；而在SSP245 和SSP585 情景下，高度适生区面积均随时间呈缩减趋势，说明长柄扁桃对气候变化的响应较敏感。

2.5 长柄扁桃适生区质心迁移路线

通过对空间分布格局变化的分析发现，在未来不同气候情景下，长柄扁桃适生区的迁移距离存在一定差异，但其迁移趋势基本一致（图7，表4）。当前长柄扁桃适生区质心位于内蒙古乌兰察布市凉城县（112.52°E，40.55°N）；而当气候情景为SSP126-2090s时，适生区质心向西北迁移78 638 m，此时适生区质心位于呼和浩特市土左旗（110.91°E，40.95°N）；当气候情景为SSP585-2090s 时，适生区质心向西迁移340 468 m，迁移至内蒙古自治区乌拉特前旗（108.50°E，40.49°N）。在未来气候变化背景下，适生区质心整体往西北迁移，并且随着温室气体排放浓度的加大，迁移距离具有向西扩大的趋势。

表4 历史、当前和未来不同气候情景下长柄扁桃适生区质心的迁移距离①Tab.4 Centroid migration distance of suitable areas of A.pedunculata under historical, current and future climate scenarios/m

2.6 环境因子分析

刀切法分析结果显示（表5），最暖季度降水量、降水量变异系数、气温变异系数、年均温、表层土壤盐基饱和度、年降水量、最湿季度平均气温等7 个环境因子对长柄扁桃地理分布的贡献率大于4%，其累计贡献率和重要置换值分别为87.9%和87.2，其中，最暖季度降水量的贡献率最大（22.7%），而降水量变异系数的贡献率次之（20.3%），表明上述7 个因子是影响长柄扁桃适生分布的主要环境因子，是构建模型的最重要环境变量。

表5 参与建模的环境因子贡献率及置换重要值Tab.5 Percentage contribution and permutation importance of environment variables for A.pedunculata in the MaxEnt model

依据环境因子响应曲线可知（图8），长柄扁桃对6 个主要环境因子的响应曲线均呈单峰（图8），表明所建模型的预测较为准确，其中，年均温、气温变异系数、降水量变异系数、年降水量、最暖季度降水量和表层土壤盐基饱和度对长柄扁桃分布的高度适宜范围分别为2.3～9.8 ℃、1 090～1 380、97%～119%/130%～150%、180～500 mm、148～296 mm 和92.9%～99.4%。

图8 长柄扁桃对6 个主要环境因子的响应曲线Fig.8 Response curves of A.pedunculata to five major environmental factors

3 讨论

3.1 优化模型的长柄扁桃适生区预测

本文基于气候、土壤和地形因子，参照最新的Checkerboard2 法，调用R 语言ENMeval 程序包，优化了MaxEnt 模型的参数。模型参数的优化可限制背景数据与校准区域，使预测分布区覆盖当前的分布点，改变正则化水平，提高模型预测与实际分布区的拟合度，并通过检查地理预测图对模型准确性进行衡量（Phillipset al.，2017）。因此，优化后的模型可有效降低其复杂性，能够显著提高模型预测的准确性（Phillipset al.，2017；赵光华等，2021）。屈振江等（2016）和叶学敏等（2019）曾采用MaxEnt 模型预测我国红富士苹果（Malus pumila）和南酸枣（Choerospondias axillaris）的适生区研究，但均未对MaxEnt 模型进行参数优化，导致预测误差较大，较难解释其研究结果。本文采用优化参数后的MaxEnt 模型，将默认参数下MaxEnt 模型的delta.AICc 值由225.87 降低至0，形成最优模型，降低了模型的复杂度，显著提高了模型预测准确性，其预测结果也与课题组前期调查结果相似（包文泉等，2018）。这与Rong 等（2019）预测青海云杉（Picea crassifolia）适生区时提出的MaxEnt 模型参数优化可使预测结果更具说服力的结论一致。可见，采用MaxEnt 模型预测物种适生区，其参数优化至关重要，是决定模型预测准确的关键因素之一。

长柄扁桃根系发达，小灌木，具有极强的抗旱、抗寒和耐贫瘠性，喜生于干旱及荒漠区的固定沙地、石砾质阳坡、山麓等地（Heet al.，2021）。通过本研究团队前期对长柄扁桃的野外调查可知，我国长柄扁桃当前实际分布区主要在内蒙古、陕西、山西、河北、宁夏和甘肃（包文泉等，2018），而本研究应用优化后的MaxEnt 模型所预测长柄扁桃适生区也恰好主要位于内蒙古、陕西、河北、甘肃、山西等地，预测与实际分布区基本一致，表明MaxEnt 模型可用于长柄扁桃分布区预测，并且优化后的MaxEnt 模型预测可靠、准确。但由于本研究仅采用了MaxEnt 模型，预测结果可能仍存在精度不够高或一些不确定性，因此后续还应采用Bioclim、Domain、Biomod2 和Garp 等模型同MaxEnt相结合，以组合模型预测加以验证本研究结果，所得结果可能会更精确。

3.2 环境因子对长柄扁桃适生分布的影响

分析表明，最暖季度降水量、降水量变异系数、气温变异系数、年均温、表层土壤盐基饱和度、年降水量和最湿季度平均气温是影响长柄扁桃分布的主要环境因子。环境因子的响应曲线显示，长柄扁桃适生区的年均温为2.3～9.8 ℃，降水量变异系数为97%～119%/130%～150%，最暖季度降水量为148～296 mm，这与内蒙古地区长柄扁桃适生区的年均温与降水特征相符（马梓策等，2019），也与长柄扁桃分布广泛的陕西榆林地区的年均温与降水量变异系数特点基本相符（刘晓琼等，2017）。前期调研发现，长柄扁桃年周期中不同物候期对土壤水分含量及环境湿度的要求存在较大差异，在果实快速生长期与花芽分化期要求降雨量较高，而花期及营养生长期对降水要求相对较低（包文泉等，2016），这与本研究中降水量变异系数和最暖季度降水量等是影响长柄扁桃分布的主要环境因子这一结果一致。此外，本研究筛选出4 个土壤因子进行建模，但仅有表层土壤盐基饱和度的贡献率大于4%（8.4%），这与长柄扁桃根系发达、抗旱、耐瘠薄，对土壤条件要求较低，能在荒山、荒地、沙地、岩石等生态脆弱区正常生长的生物学特性相一致（乌云塔娜等，2014）。

本研究表明，长柄扁桃主要适生区分布在我国北方地区，温度和降水是影响长柄扁桃适生的主要因子，其贡献值显著大于土壤和地形因子。但由于土壤和地形因子是较复杂的环境因子，并且二者还可对温度和降水进行再分配，因此，温度、降水、土壤、地形因素对长柄扁桃适生分布的影响均不能忽视，针对四者间存在何种关系以及具体影响机制是什么，后期还需进一步研究和分析。

3.3 气候变暖对长柄扁桃适生分布的影响

通过对长柄扁桃历史、当前和未来潜在分布区的模拟表明，气候变暖对长柄扁桃适生区分布的影响显著。从长柄扁桃历史分布区的推测可知，在末次间冰期，长柄扁桃分布于我国北方大部分地区，面积约为当前适生区的1.5 倍；而在末次盛冰期，长柄扁桃适生区从东部向西部严重收缩，集中于黄土高原和华北平原地区，表明末次盛冰期的气候变化严重缩小了长柄扁桃的分布范围，使其仅能在中高纬度地区的适生区得以留存；此后，全新世中期是距今最近的一次大暖期，气候变暖使长柄扁桃适生区较末次盛冰期缩减近一半，并向高纬度区迁移，新增区集中于内蒙古东北部。这一结果与众多植物历史分布区研究相符，物种的适生区范围与气温升高呈负相关（郑卓等，2013），在气温升高时向高纬度地区迁移（刘金陵等，2004；白伟宁等，2014）。本研究中长柄扁桃在其历史上呈动态分布，从末次间冰期至末次盛冰期，长柄扁桃的分布格局变化最显著，由西南往东北收缩；而从末次盛冰期至全新世中期，则由东向西收缩；可见，冰期与间冰期交替的历史气候事件促进了长柄扁桃当前分布格局的形成，但其影响程度和作用机制还尚未清楚，后期还需更深入探究。

与当前相比，未来高温高湿环境使长柄扁桃适生分布区不断向西北扩张，但扩张面积远低于丧失面积，表明长柄扁桃对气候变化响应较敏感，尤其在高浓度温室气体情景下更敏感。这一结果与绝大数北方落叶林木对气候变化的响应一致（Zhaoet al.，2021）。

长柄扁桃空间分布格局的变化分析也表明，随着气候变暖，适生空间有向西北迁移的趋势；并且在温室气体的低、中排放浓度下，适生区质心的迁移距离较小；而在温室气体高排放浓度下，适生区几何中心的迁移距离较大。这与温带树种适生分布区在气候变暖下往高纬度、高海拔区迁移的趋势相符（Yeet al.，2020；张殷波等，2019）。因此，未来由大气CO2浓度增加引起的全球气候变化同样是影响长柄扁桃分布的主要因素，尤其是低纬度、低海拔区的长柄扁桃分布面积较当前锐减，具有适生生境丧失的可能。

4 结论

本研究对MaxEnt 模型默认参数进行优化，得到最优模型的特征组合为线性、二次型、片段化、乘积型与阈值性特征，其调控倍率为1.5。优化后模型的训练遗漏率和复杂度较低，预测结果与实际分布一致，准确可靠。当前我国长柄扁桃适生分布区主要在内蒙古高原和黄土高原地区，属于中高纬度区。在历史、当前和未来气候情景下，长柄扁桃对气候变化的响应较为敏感；在未来不同气候变化情景下，长柄扁桃适生区有缩减的趋势，由低海拔和低纬度向中高纬度和高海拔地区迁移，尤其在温室气体高排放浓度下迁移距离更大。

在全球气候增温增湿背景下，长柄扁桃在吉林、辽宁、赤峰、通辽等高纬度地区的适生区将大面积消失，而在中高纬度和高海拔（内蒙古呼和浩特、鄂尔多斯、锡林郭勒，陕西延安）有新增适生区出现。随着未来温室气体排放浓度增加，长柄扁桃当前适生区的东南部区可能不再适宜。因此，长柄扁桃产业的今后发展，应选择当前和未来均适宜的区域，如重点选择内蒙古的‘呼包鄂’，以及内蒙古、陕西、山西、河北与宁夏交接区域，可避免未来气候变暖带来的损失。