气候变化背景下药用植物牡丹的潜在适生区分析

黄雪梅 ,黄怡 ,贾泽旭 ,陈德朝*

1.四川省林业科学研究院森林和湿地生态恢复与保育四川省重点实验室，四川 成都 610081；

2.绵阳师范大学生命科学与生物技术学院,四川 绵阳 621000；

3.西藏大学,西藏 拉萨 850000

全球变暖趋势日益显著，中国气候变化趋势近似全球气候变化[1]。中国国家气候中心统计表明，极端气候出现的频次和强度均显著升高[2]，在此背景下，农业、林业、畜牧业等行业均受到严重影响[3]，药用植物也不例外。几十年来，植物与环境对药用植物生长备受研究人员关注[4]。研究表明，环境因素和人类活动不仅影响药用植物分布，在其活性成分形成上，也有着极其重大作用[5]。气候严重影响了植物的生长和繁殖，成为决定植物分布的主要变量[6]。受气候全球变暖威胁，许多药用植物地理分布面急速下降缩小，部分药用植物将面临灭绝的风险[7-9]。到目前为止，已有研究采用科学建模的方法对药用植物的生境分布进行预测[10-11]。因此，有必要利用生态模型预测适生区分布、适宜栽培区域以及决定药用植物生境分布和活性成分增强的环境因子。

在2004 年，Phillips 等[12]通过以最大熵理论作为基础，开发物种分布模型软件MaxEnt。MaxEnt 模型自问世以来，被国内外学者广泛应用，特别是在植物适生区预测效果方面在业内获得了极高的评价，例如桫椤（Alsophila spinulosa）[13]、梓叶槭（Acer amplumsubsp.catalpifolium）[14]和藤枣（Eleutharrhena macrocarpa）[15]等等。目前，常用物种分布模型有10 余种，由于不同模型的理论基础不同，因此它们的模拟效果和预测表现差异较大，Petltpierre等[16]和江一帆[17]研究表明，在众多物种分布模型中，MaxEnt 模型预测精度最高，模拟效果更好。因此，选取MaxEnt 作为模拟软件预测目标物种的潜在分布并分析气候变量对物种分布的影响。

牡丹（Paeonia suffruticosa）是一种落叶亚灌木，在生态、经济和观赏上也有着独特价值，是重要的植物资源；除此之外，牡丹也作为中药材被广泛使用。野生牡丹大都分布在偏远地带，当地居民肆意挖掘现象泛滥，导致野生牡丹数量不断减少[18]。所以通过探究牡丹的适宜分布将有效指导引种栽培，对生产实践具有重要意义。牡丹作为药用植物也是一种特殊的植物资源，它的种植、生长均受到气候变化的影响。如何有效应对未来气候变化和人为干扰，是药用植物引种栽培规划中的关键问题。从生态位理论的角度对牡丹地理分布进行模拟，分析目标物种地理分布与气候和人类活动的关系，研究环境变量和人类活动对牡丹分布的影响，预估气候变化对牡丹潜在分布区的影响，研究结果可为牡丹应对气候变化导致适生区的改变提供数据支撑，并为其引种栽培规划提供理论依据。

1 材料及方法

1.1 物种数据来源

牡丹的分布信息：实地调查和查询中国数字植物标本馆（CVH）收集牡丹的分布点。参考Elith 等人[19]的研究方法，通过ArcGIS 的缓冲区分析功能，计算网格中心与分布点之间的距离，在同一网格内，只保留1 条离中心最近的记录，共得到牡丹分布点92 个（见图1）。

图1 牡丹样本分布点位图Fig.1 Species occurrence records of P. suffruticosa

1.2 变量的来源与处理

气候数据来源于Worldclim 数据库，该数据是以全球的气象记录信息为基础数据，整合插值生成的全球气候栅格数据，包括19 个生物气候因子。人类活动强度数据（HAI）网格，从2009 年《人类足迹》获得[20]。气候系统模型经过相互比较选择耦合模式为第6 次（CMIP6）北京气候中心气候系统模型作为未来气候模型，因为其更接近真实情况[21]。上述数据的空间分辨率均为2.5 arc-minutes。

1.3 建模变量筛选

采取以下步骤筛选建模变量：首先，将所有初始变量数据和物种的分布数据全部导入到MaxEnt，统计出各变量之间的平均初始贡献，去除贡献率相对极高低的变量。其次，使用ArcGIS 软件提取剩余变量属性值。然后，通过Spearman 对其进行相关分析。最后，比较变量之间相关系数，对于相关系数高于0.8 的变量，先根据变量的生态学意义筛选，仍无法确定时则比较在初始模型中的贡献率。最终针对牡丹筛选出8 个环境因子（见表1）。

表1 参与牡丹建模运算的环境变量Tab.1 Environmental variables involved in P. suffruticosa model operation

1.4 模型精确性验证

选 用 ROC（receiver operating characteristic）曲线法评价MaxEnt 模拟结果的准确性。ROC曲线其曲线下面积（area under curve，AUC）的大小是评估模型预测是否准确的指标[22]。AUC 值从0—1 与预测效果好坏呈正相关，其越接近1 时模型的预测精度越高。将评估模型的精度划分为4 个等级，其中较差等级（0.60≤AUC＜0.70）、一般等级（0.70≤AUC＜0.80）、好等级（0.80≤AUC＜0.90）、很好等级（0.90≤AUC≤1.00）。

1.5 适生等级划分

应用ArcGIS 软件将MaxEnt 软件输出的ASCⅡ格式文件转为Raster 格式文件。再用ArcGIS 的统计分析功能计算在各省（区、市）的分布面积。MaxEnt模型默认适生等级为10 级[22]，根据牡丹在中国的实际发生情况并结合相关文献，利用ArcGIS 的Reclassify功能对连续分布概率（P）进行重分类，其中浅黄色为不适生区，0.0≤P＜0.1；浅绿色为低度适生区，0.1≤P＜0.3；淡蓝色为中度适生区，0.3≤P＜0.6；深蓝色为高度适生区，0.6≤P。

1.6 未来高适生区质心中心计算

统计高适生区在未来不同时期面积变化情况以及其中心点的位移，计算公式如下：

式中：t 表示不同的时段；I 为适生区的单位栅格的数目；si(t)为t 时段单位栅格的面积；s(t)为t 时段适生区的总面积；(Xi(t)，Yi(t))为t 时段适生区单位栅格的质心坐标；(X(t)，Y(t))为t 时段适生区的质心坐标。

式中：D 为t 时段到t+1 时段适生区质心的位移距离；θ为t 时段到t+1 时段适生区的质心位移方向，0°＜θ＜90°表示位移方向为东北，90°＜θ＜180°表示位移方向为西北，180°＜θ＜270°表示位移方向为西南，270°＜θ＜360°表示位移方向为东南。

2 研究结果

2.1 模型精准度

基于92 条牡丹分布记录利用MaxEnt 模型对牡丹在中国的潜在适生区进行模拟预测，依据牡丹测试集数据的AUC 值为：0.969 和训练集数据AUC 值为：0.979，可以得出牡丹模拟的效果很好。其ROC曲线如图2 所示。

图2 牡丹Maxent 模型的ROC 曲线Fig.2 Receiver operating characteristic curve of P. suffruticosa Maxent model

2.2 影响牡丹潜在适生区的主导影响因子

目前研究者对于主导因子数量的确定观点是不一致的，根据贡献率确定主导因子是目前多数人的看法，将达到一定程度的环境因子作为主导因子，但是在程度的选择上是主观的，所以出现了不同的标准。选择贡献率排前三的环境因子作为主导环境因子。经过运算，从牡丹的运算结果来看8 个环境因子中年降水量（Bio12）、人类活动强度（HAI）和最冷月最低气温（Bio6）贡献率分别为32.80%、30.70%和20.8%，三者累计值达84.30%，排在前三位，作为主导环境因子（见表2）。

表2 影响因子对于牡丹的贡献率Tab.2 Contribution rate of influencing factors to P. suffruticosa

通过刀切法（Jackknife）进行检验可知，在使用单一的环境因子变量时，年降水量（Bio12）、最冷月最低气温（Bio6）和人类活动强度（HAI）三个环境因子是对牡丹正规化训练增益影响最大，说明了这些环境因子变量在制约牡丹生长上确实占据了主导地位（图3）。

图3 环境变量应用刀切法对牡丹的检验结果Fig.3 The jackknife test result of environmental variables for P. suffruticosa

2.3 牡丹潜在分布及分布适宜性评价

牡丹潜在适生区预测模型的拟合效果结果很好。从预测结果可知牡丹潜在适生区结果如图4 所示。

图4 牡丹当代潜在分布及分布Fig.4 Potential distribution of P. suffruticosa under contemporary climate conditions

牡丹适生区主要位于华南地区北部、西南地区的东部、华北地区的南部和华中华东地区的绝大多数区域的各省，总适生区面积为324.35×104km2。高适生区面积占总适生区的7.37%，面积为23.90×104km2，主要呈饼状分布在川东地区，另一部分呈细条状或者斑块分布在渝西南、黔西北、湘西南等地区。中适生区面积占总适生区的31.89%，面积为103.44×104km2，主要包围在高适生区的外围呈饼状分布在川东、宁南、陕东南、黔和渝等地区，另一部分呈细条状或者斑块分布在京南、湘南等地区。低适生区面积占总适生区的60.74%，面积为197.01×104km2，主要分布在川西、滇北、桂北、藏东南、豫、鄂等地区，超过总适生区面积的一半。综合来看，牡丹的高适生和中适生区呈饼状分布于川、渝、黔等地区，与其实际分布十分吻合，进一步表明模拟效果准确度高。

2.4 未来气候变化对牡丹地理分布的影响

利用Maxent 模型预测了牡丹在现代气候条件下2050 年和2070 年SSP1-2.6、SSP2-4.5、SSP5-8.5 不同排放情景下牡丹潜在适生区，从而得到了气候变化情景下的牡丹潜在适生区（图5）和气候变化情景下牡丹不同等级潜在适生区的面积（见表3）。在三种SSPs（SSP1-2.6、SSP2-4.5 和SSP5-8.5）情 景下，牡丹低适生区面积在2050 年代、2070 年代均有不同程度的增加，牡丹中适生区面积除在2070 年代SSP2-4.5 情景下有微小减少外，在其他情景下也均有不同程度的增加，但是牡丹高适生区面积在2050 年代、2070 年代均有不同程度的减少，减少最大的是在2070 年代，SSP2-4.5 情境下，高适生区面积减少了47.91%，减少最小的是在2070 年代，SSP2-2.6 情境下，高适生区面积也减少了22.22%，说明气候变化对于牡丹高适生区影响较大。牡丹总适生区面积在2050 年代、2070 年代均有不同程度的增加，说明除了一部分高适生区转变为中、低适生区外，新出现了很多适生区，尤其是低适生区（见图5 和表3）。

表3 不同时期牡丹潜在适生区面积/104 km2Tab.3 Suitable areas of P. suffruticosa under different climate change scenarios/104 km2

图5 气候变化情景下的牡丹潜在适生区域分布Fig.5 Potential distribution of P. suffruticosa under future climate change scenarios

与目前的分布相比，在2050 年代三种SSPs（SSP1-2.6、SSP2-4.5 和SSP5-8.5）情景下总适生区面积将分别增加13.54%、12.25%和11.64%，高适生区面积将分别减少30.13%、23.35%和37.91%，中适生区面积将分别增加4.52%、16.67%和5.24%，低适生区面积将分别增加23.57%、14.24%和21.01%（见表3）。与目前的分布相比，在2070 年代三种SSPs（SSP1-2.6、SSP2-4.5 和SSP5-8.5）情景下总适生区面积将分别增加12.33%、17.90%和26.15%，高适生区面积将分别减少22.22%、47.91%和42.05%，低适生区面积将分别增加17.59%、35.73%和38.15%（见表3），中适生区面积在2070 年代SSP1-2.6 和SSP5-8.5 情景下，分别增加10.3%、19.04%，在SSP2-4.5 情景下，减少0.85%（见表3）。

综上所述，在未来气候变化情景下牡丹的高适生区面积总体呈现减少趋势，中、低适生区面积总体呈现增加趋势，低适生区面积增加较大，且高适生区面积在2050 年代、2070 年代，这3 种SSPs（SSP1-2.6、SSP2-4.5 和SSP5-8.5）情景下均显著减少，高适生区面积减少最多得情境下超过了47%，接近一半。

2.5 未来气候变化下牡丹高适生区质心移动轨迹

对牡丹在不同时段及其气候变化情景下高适生区质心移动轨迹图（见图6）可知，从基准气候到2 050时段再到2 070 时段，牡丹高适生区质心，在3 种浓度排放情景下，均呈向高纬度和东北方向迁移的趋势（见图6），在从基准气候到2 050 时段质心移动幅度均较大，所以高适生区面积减少也较为显著。从2 050时段到2 070 时段，在高浓度排放情景下，牡丹高适生区质心迁移幅度较大，所以高适生区面积减少也较为显著，在中排放情景下，牡丹高适生区质心迁移幅度变化不大，所以高适生区面积较2 050 时段减少不明显。

图6 气候变化情景下的牡丹高适生区质心移动轨迹Fig.6 Movement trajectory of the mass center of highly suitable areas for P. suffruticosa under climate change scenarios

3 讨论

3.1 牡丹当前的潜在分布区

在过去几十年里，全球气候和环境变化研究到科学领域受的重视程度日益加强，生态学家已意识到生物多样性丧失所带来的危害，并采取生态保护作为环境研究领域的热点[23-24]。研究表明（见图5），牡丹的适生区总面积为324.35×104km2，主要位于华南地区北部、西南地区的东部、华北地区的南部和华中华东地区的绝大多数区域的各省，其中高适生区23.90×104km2，绝大部分呈饼状分布在川东地区，另一部分呈细条状或者斑块分布在渝西南地区、黔西北地区、湘西南地区、晋南地区等。根据植物标本数据库表明，牡丹在渝西南地区、川东地区、晋南地区、苏中地区、鲁北地区、湘西南地区等均有分布，上述区域均位于预测的适宜区域内，说明牡丹的预测结果可靠。

3.2 药用牡丹与环境变量之间的关系研究

在限制物种潜在分布上，气候因子是在宏观尺度上的关键因素，植物与气候之间相互作用研究是生态学的一个重要方向。通过MaxEnt 模型所预测结果可知道，制约牡丹潜在地理分布的主导环境因子是最冷月最低气温（Bio6）、年降水量（Bio12）和人类活动强度（HAI）。本研究表明，随最冷月最低气温的升高牡丹存在概率先呈现上升趋势，当最冷月最低气温持续升高时牡丹存在概率呈现下降趋势，说明最冷月最低气温过高或者过低到一定程度时会超出牡丹生存阈值，使其在部分区域生存概率下降，年降水量多少也制约牡丹存在概率，随着年降水量增加牡丹存在概率也随之正向增加，当随着年降水量的持续增加牡丹的存在概率又快速下降，这说明在降水量上，牡丹只能承受一定范围内的降水量，年降水量过低或过高都会降低牡丹存在概率，这也与牡丹的喜欢寒冷和气候湿润，耐寒，但不耐旱的生活习性一致，这与目前相关文献的描述一致，同时也说明其对温度和降水的要求较高，进一步验证了本研究结果的准确性。王茹琳等人[25]探究了中华猕猴桃（Actinidia chinensis）在中国的潜在分布以及应对气候变化的响应，赖文峰等人[26]探究了藏药桃儿七（Sinopodophyllum hexandrum）的潜在适生区，表明植物适宜分布区域受气温、降水量等气候影响剧烈，都会随着这些变量的变化而变化，一部分适生区将会消失，一部分新的适生区将会出现，表明温度和降水量对牡丹的适生区有着较为重要的影响。Xu 等人[27]的研究发现植物分布对于人类活动有着不同的反应。本研究表明，人类活动强度的增加降低了物种潜在分布概率，使其分布区域变得支离破碎，这意味着该物种对人类活动干扰比较敏感。虽然人工栽培牡丹可满足目前对于牡丹丹皮和根的部分需求，但野生采伐仍然是药物市场上最重要的来源，因此，积极开展牡丹人工种植工作。Thomas等人[28]研究表明，在气候变暖背景下，一部分物种将会灭绝，但是也有很大一部分物种的生长和分布受到不同程度的发展，这表明，气候变暖对于物种的生长和分布是一把双刃剑；显然牡丹在气候变暖背景下，其生长和分布得到了一定程度的发展。

3.3 未来气候变化情景下牡丹潜在地理分布的变化

通过MaxEnt 模型预测未来3 种不同浓度（SSP5-8.5、SSP2-4.5、SSP1-2.6）排放情景下牡丹潜在地理分布，对两者进行空间叠加，得出牡丹潜在分布在未来3 种不同浓度排放情景下的变化（见图7）。

图7 气候变化情景下的牡丹适生区变化Fig.7 Changes of the potential distribution of P. suffruticosa under climate change scenarios

预测结果可知，未来两个时段，牡丹潜在地理分布在任何排放浓度下，高适生区面积呈减少趋势，低适生区面积呈增加趋势，中适生区面积除在2070 年代中浓度排放情境下有微小减少外，在其他情景下也均有不同程度的增加（见图5）。牡丹高适生区面积在2050 年代、2070 年代均有着显著的减少，特别是在高浓度排放情境下高适生区面积减少更为显著，牡丹中、低适生区面积在2050 年代、2070 年代总体上呈现增加趋势，可能是由于气候变化背景下，特别是在高浓度排放情境下，气温值、降水量和人类活动强度在分布地区超过了牡丹适宜生存的阈值，使得牡丹在部分区域存在概率下降，高适生区面积减少显著，同时又有一部分气温值、降水量和人类活动强度达到了适宜牡丹生存条件，成为牡丹的潜在分布区，且呈增大的趋势（见图5和图7）。牡丹高适生区质心在3 种浓度排放情景下，均呈向高纬度和东北方向迁移趋势（见图7），总体来说迁移幅度都很大，只有在2 050 时段到2 070时段，在中排放情景下，牡丹高适生区质心迁移幅度变化不大，这就表明，未来将会出现很多适宜牡丹生长的新区域，同时也将有更多区域将会变得不再适宜牡丹生长；总的来说，新出现适宜牡丹生长区域比消失老区域多，但是高适生区新出现的区域比消失的老区域少很多。囊括全球表面达20%样本区域的生物灭绝风险研究表明，将会有一部分物种将在中等排放浓度下的2050 年灭绝，但是也有一部分物种的生长和分布均将会在这期间受到不同程度的发展。这表明了，并不是所有物种生长和分布都会随着气候变暖开始消亡或者开始增加，对于物种的生长和分布气候变暖是呈两面性的[29]。

预测结果可以得出，在2050 年代，牡丹在高浓度排放情景下丧失适生区面积最大，为12.00×104km2；在低浓度排放情景下新增的适生区面积最大，为49.95×104km2；但是在3 种浓度排放情景下，新增的适生区面积相差不大。在2070 年代，牡丹在高浓度排放情景下，丧失和新增的适生区面积最大，分别为41.52×104km2和125.92104km2（见表4）。尽管在未来气候情景下，牡丹的总适生区呈增长趋势，但是其高适生区面积却大幅度减少，所以对于牡丹的保护已刻不容缓。在低浓度排放情境下，2050 年代牡丹丧失的适生区主要位于桂北地区、粤北地区和晋中地区，新增适生区主要位于青东南地区、藏东南地区、蒙东地区、黑南地区和吉西北地区；在2070 年代，与2050 年代相比，在福建省西北部新增了牡丹适生区的丧失区域，与新增适生区大致相同，呈增大趋势。在中浓度排放情景下，2050 年代牡丹丧失的适生区主要位于滇东北地区、桂北地区、粤北地区和晋中地区，新增适生区主要位于青东南地区、藏东南地区、蒙东地区、黑南地区和吉西北地区；在2070 年代，与2050 年代相比，牡丹适生区丧失区域呈一幅度扩张，晋大部分地区、闽大部分地区、皖中地区都新增了牡丹适生区的丧失区域，与新增适生区大致相同，呈增大趋势。在高浓度排放情景下，2050 年代牡丹丧失适生区主要位于滇东北地区、桂北地区、粤北地区和晋中地区，新增适生区主要位于晋北地区、冀北地区、青东南地区、藏东南地区、蒙东地区、黑南地区和吉西北地区；在2070 年代，与2050 年代相比，牡丹适生区的丧失区域呈大幅度扩张，向北边快速扩张，桂、粤、晋等省基本丧失了牡丹的适生区域，新增适生区大致相同，范围呈大幅度向北增大趋势。牡丹的适生区有部分将会在未来消失，也会在未来出现更多新区域适宜牡丹生长，虽然新增的区域远多于消失区域，但是高适生区却大幅度减少。就纬度变化而言，牡丹在3 种浓度排放情境下均向高纬度东北方向迁移，这可能与牡丹喜欢寒冷和气候湿润，耐寒，但不耐旱的生活习性相关，这也与植物志上的描述以及赵宣的研究是基本一致的[30]。张华等人[31]下研究表明，随着气候变暖，物种会向着不同纬度地区进行迁移，本研究中的牡丹在不同浓度排放情景下，也在向着高纬度地区进行迁移，和上述的研究结果规律一致。牡丹同其他物种一样，在气候变暖背景下，其潜在地理分布区部分将会新增和消失，同时潜在地理分布也将会发生转移。

表4 气候变化情景下牡丹在适生区的变化/104 km2Tab.4 Future changes of the suitable habitat areas for P. suffruticosa under climate change scenarios/104 km2