Prediction of the potential distribution pattern of Pinus sylvestris var. mongolica in China under climate change

2018—2021年内蒙古自治区林业和草原监测规划院通过实地调查，在内蒙古自治区范围内樟子松分布区域共获取130个分布数据；结合中国数字植物标本馆(CVH)、中国国家标本植物平台(NSII)、全球生物多样性信息网(GBIF)及其他相关文献共获取樟子松分布点数据115个，经检索汇总后共保留245条记录。将文献记载及实地调研数据整合，配合ArcGIS软件对样点经纬度进行筛选剔除。为避免樟子松群体分布得不均匀，以经纬度2.5′×2.5′为1个网格单元，每个单元筛选1条信息记录，截至2021年最终保留200个样点进行后续拟合分析。

本研究选取空间分辨率均为30″的19个生物气候因子和1个地形因子作为环境因子(表1)。气候数据来源于世界气候数据库(https://www.worldclim.org)，地图数据来源于国家基础地理信息中心(http://www.ngcc.cn/ngcc/)的中国标准地图。

未来气候数据选取第6次国际耦合模式计划(CMIP6)的国家气候中心中等分辨率气候系统模式(BCC-CSM2-MR)。此模式为中国未来的气候模式。为更均衡合理模拟未来时期(2050s和2100s)樟子松潜在分布，同时考虑社会经济共享路径(SSP)和代表性浓度路径(RCP)，选择SSP-RCP2.6 (SSP126)、SSP-RCP4.5 (SSP245)和SSP-RCP8.5 (SSP585)等3个情景对樟子松未来分布区进行预测。SSP126情景使用可持续发展路径，代表温室气体排放量处于较低水平；SSP245情景使用中度发展路径，代表温室气体排放量处于中等水平；SSP585情景使用以化石燃料为主的发展路径，代表温室气体排放处于高水平^[18]。

为避免20个生物气候变量间的相似性过高导致模型过度拟合，通过方差膨胀因子(VIF)和Pearson相关性检验，结合MaxEnt模型预实验中各个因子贡献率，选择相关系数小于0.8、VIF小于10且对樟子松分布具有较高贡献率和生态学意义的因子。最终保留8个气候因子和1个地形因子进行建模，包括年平均气温、最热月最高气温、最冷月最低气温、年降水量、最冷季度平均气温、气温季节变动系数、最干季度平均气温、降水量季节性变化和海拔。

本研究利用MaxEnt模型预测当前和未来时期樟子松在中国的潜在分布。将樟子松200个分布点和9个环境因子导入MaxEnt 3.4.1软件，java环境下运行，随机选取25%的样本作为测试数据集，75%的样本作为训练数据集^[19]，以Bootstrap模式，重复迭代10次，建模结果以Logistic形式输出。

MaxEnt模型准确度可通过受试者工作特征曲线(ROC)进行检验。用该曲线下的面积(AUC)表示所建模型的精度，AUC通常为0.5~1.0，其值越大，表明模型预测越精准。其中AUC在0.5~0.6时表明模型预测失败，0.6~0.7时表明模型预测效果较差，0.7~0.8时表明模型预测效果一般，0.8~0.9时表明模型预测效果良好，0.9~1.0时表明预测效果极好^[20]，可以较准确地反映物种潜在分布。

本模型采用优化方法为Checkerboard 2法，调整调控倍率(RM)和特征组合(FC) 2个参数改变模型正规化水平。通过改变特征数量、组合方式以及RM数值，将MaxEnt模型中的5种特征：线性特征(L)、二次型特征(Q)、片段化特征(H)、乘积型特征(P)和阈值型特征(T)，与16个RM值(0.5~8.0，以0.5为间隔梯度)排列组合出148种组合方式。利用R语言中的ENMeval程序包对其进行优化测试，结果中各组合的ΔAICc值和10％测试遗漏率越低，表示模型预测精准度越高^[21−22]。

运用ArcGIS将模型运行后的数据进行适宜区划分和可视化处理，根据软件中提供的多种分类方法预测樟子松适宜性阈值。根据冉巧等^[7]对孑遗植物银杉进行适宜区划分的方法，本研究使用自然断点法划分樟子松适宜指数。参照张伟萍等^[8]对祁连圆柏适宜区划分的等级，将樟子松生境划分为不适宜区、较不适宜区、一般适宜区和高度适宜区(适宜指数分别为0~0.1、＞0.1~0.2、＞0.2~0.6、＞0.6~1.0)。对比樟子松不同时期的适宜区分布差异，获得樟子松不同适宜区空间分布图；利用ArcGIS 10.4.1软件中SDMtoolbox 2.0工具包，计算樟子松在不同时期的适生区质心位置得到樟子松分布区的迁移方向，并计算其质心迁移距离。

根据200个樟子松分布点和9个环境因子利用MaxEnt模型对樟子松的潜在分布区域进行模拟预测。在MaxEnt为默认参数设置时，RM为1.0，FC为LQHPT，∆AICc为29.77。当RM为1.5， FC为LQHPT时，∆AICc为0，此时模型最优，10%的训练遗漏率降低4.63% (表2)。因此选取RM为1.5， FC为LQHPT为模型最终参数。该参数下10次模拟训练AUC均值为0.964 (图1)，说明预测结果精确。

Figure 1.  ROC response curve under MaxEnt model

樟子松当前气候条件下的潜在适生区集中在53.55°N以南，135.09°E以西，总适生区面积(一般适生区面积与高度适生区面积之和)为645 104.2 km²，主要位于大兴安岭地区、内蒙古中部、黑龙江西北部，在山西、河北、吉林、辽宁也有少量分布^[23]，其中，高度适生区面积为147 638.9 km²，一般适生区面积为497 465.3 km²。樟子松高度适应区绝大部分在内蒙古东北部，黑龙江西北部呈零星分布。

未来时期不同气候情景下，樟子松适生区面积较当前时期有不同程度减少。2050s，樟子松总适生区在SSP126、SSP245和SSP585气候情景下降35.7%、52.5%和82.6%；2100s与2050s相比，总适生区面积在SSP126、SSP245和SSP585情景下呈下降趋势，降幅分别为62.1%、37.4%和100.0%。SSP126-2050s情景下，樟子松高度适生区面积相较当前大幅下降，面积约为当前的6.6%，仅存在大兴安岭较高海拔地区；SSP245-2100s气候情景下，樟子松高度适生区将完全消失。

在3种气候条件下樟子松适生区面积均下降，SSP245和SSP585气候情景下樟子松适生区对气候的响应最为明显 (表3和表4)。SSP126气候情景下，2050s樟子松适生区面积变化幅度最小，适生区面积下降230 225.7 km²，降幅为35.7%，主要变化区域位于呼伦贝尔市西南部；2100s适生区面积下降487 968.8 km²，降幅为75.6%，主要在兴安盟和锡林郭勒盟的东北部出现较大面积消失。SSP245气候情景下，2050s樟子松适生区面积下降338 449.4 km²，降幅为52.5%，其变化区域主要位于内蒙古自治区呼伦贝尔市西东两侧地区和赤峰市；2100s其适生区面积下降453 211.8 km²，降幅为70.3%，内蒙古自治区东部和黑龙江省的大部分地区樟子松适生区消失。SSP585气候情景下，樟子松适生区面积变化幅度最大，2050s樟子松适生区面积下降532 882.0 km²，比当前减少82.6%，仅在呼伦贝尔市中部有樟子松分布；2100s樟子松适生区将在中国完全消失。

在空间分布格局方面，未来时期不同气候情景下樟子松适生区的迁移位置存在差异，但总体上呈现出向西北和西南方向迁移(表5)。当前时期樟子松适生区的质心在内蒙古兴安盟阿尔山市(47.09°N，120.25°E)。气候情景为SSP126-2100s时，樟子松适生区向西北方向迁移，此时樟子松适生区质心位于内蒙古自治区呼伦贝尔市鄂温克族自治旗(48.41°N，120.17°E)，迁移距离为161.47 km；气候情景为SSP585-2050s时，樟子松向西北方向迁移，此时其适生区质心位于内蒙古自治区呼伦贝尔市新巴尔虎左旗(47.88°N，119.06°E)，迁移距离为126.36 km；气候情景为SSP245-2100s时，樟子松向西南方向迁移，此时其适生区质心位于内蒙古自治区锡林郭勒盟东乌珠穆沁旗(45.66°N，117.39°E)，迁移距离为270.80 km。未来时期不同气候变化情景下，全球增温增湿使得樟子松在中国适生区的质心整体向西北和西南迁移，质心随着适生区分布范围消减可能消失。

根据Jackknife法分析结果(表6)，樟子松潜在地理分布贡献率超过5%的环境因子分别为最冷季度平均气温(22.1%)、降水量季节性变化(19.8%)、最冷月最低气温(18.1%)、气温季节变动系数(17.8%)、最干季度平均气温(9.0%)、最热月最高气温(6.1%)。这6个环境因子总贡献率达92.9%。

樟子松分布主要环境因子与潜在分布区的关系如表7所示。SSP126气候情景下，2050s和2100s樟子松生境适宜度比当前降低了0.38和0.54；SSP245气候情景下，2050s和2100s樟子松适宜度比当前降低了0.46和0.55；SSP585气候情景下，樟子松生境适宜度急剧下降，2050s下降了0.55，到了2100s下降0.67，降幅为91%。

200个樟子松分布点的年平均气温与其生境适宜度变化正好相反，2050s，SSP126、SSP245和SSP585情景下年平均气温分别增加19.5%、56.6%和101.3%；2100s与2050s相比，SSP126、SSP245和SSP585情景下年平均气温分别增加15.9%、48.1%和90.7%。

200个樟子松分布点的年降水量在SSP126、SSP245和SSP585情景下均为递增趋势，且年降水量与年平均气温变化趋势相同，都与樟子松物种适宜度变化相反。2050s，SSP126、SSP245和SSP585情景下年降水量分别增加5.67、20.32和44.58 mm；2100s与2050s相比，SSP126、SSP245和SSP585情景下年降水量分别增加29.58、18.27和45.62 mm。

本研究基于气候和地形因子的研究，应用ENMeval数据包中的Checkerboard 2法来优化模型^[23]。作为北方主要造林树种，樟子松集中分布在中国北方地区，具有喜光、耐寒、耐旱、耐贫瘠和抗病能力强的生长特性。在中国的大兴安岭北部，西起莫尔道戈、金河、根河，东到新林、呼玛以北有连续的分布，在南伊图里河的银河、免渡河、阿尔山等地呈带状或块状分布着樟子松纯林以及其与兴安落叶松Larix gmelinii构成的混交林^[24]。自1955年樟子松被成功引种到辽宁省阜新市彰武县章古台以来，已推广到科尔沁沙地、浑善达克沙地、毛乌素沙地等地。本研究基于优化的最大熵模型预测结果显示：当前樟子松高度适生区也位于上述地区，预测结果与实际高度一致，表明经过优化后的MaxEnt模型用于樟子松适生区的预测结果精准可靠。

Jackknife法分析表明：年平均气温、最热月最高气温、最冷月最低气温、年降水量、最冷季度平均气温、气温季节变动系数、最干季度平均气温、降水量季节性变化和海拔是影响樟子松分布的主要环境因子。当樟子松存在概率大于0.6时，对应的环境因子最适合植物生长。主导环境因子分析表明：气温和降水相关因子对樟子松分布贡献率达到98.3%，对其分布起决定性作用。这与赵晓彬^[25]研究中樟子松适生于年平均气温−3.7~−1.5 ℃、年降水量375 mm的结论相符。本研究预测樟子松的最适宜海拔为800 m左右，与戴继先等^[26]研究发现樟子松天然分布在300~900 m的阳坡或向阳坡的结论一致。水热条件一定程度上决定了樟子松潜在地理分布，地形在一定程度上也会对水热资源进行二次分配，因此，气温、降水和地形对樟子松分布格局的相互作用都不能忽视。

2050s和2100s不同二氧化碳排放情景下，高温高湿使樟子松适生区不断减少。2100s樟子松适生区显著减少，其降幅明显大于2050s。未来不同的气候情景下，樟子松适生区的空间分布变化随着气候变暖加剧，整体向西北高纬度和西南季节性降水量增加的地区收缩。SSP126情景下，2050s适生区质心迁移距离最短，说明气候小幅变暖对适生区分布范围影响不大。SSP245情景下樟子松适生区质心点受季节性降水量影响整体向西南方向降水量充沛地区迁移，2100s与2050s相比其质心点迁移距离较大。SSP585情景下温室气体排放增加，2050s樟子松适生区质心点明显向西北高纬度地区迁移，2100s适生区在中国消失。张日升等^[27]研究表明：随着温室气体排放等环境问题加剧，年平均气温达7 ℃左右时樟子松将出现生长不良，导致大面积退化。吴祥云等^[28]也曾提出：樟子松人工林引种到年平均气温为7.1 ℃的章古台会发生“早衰”现象，导致面积减退。未来气候情景下，温室气体排放，气温升高导致樟子松适生区面积逐步北移，在中国分布减少甚至至消失。综上所述，樟子松迁移趋势与CHEN^[29]研究中国寒温带森林树种地理分布在未来全球变暖下往高纬度和降水量充沛地区迁移的趋势相符。

本研究通过对MaxEnt模型的优化，发现樟子松潜在分布最优模型的特征组合为片段化特征(LQHPT)、乘积型特征(P)和阈值性特征(H)，调控倍率为1.5。此模型复杂度低，预测结果较准确，预测樟子松适生区范围与实际分布情况基本一致。中国樟子松当前高度适生区集中在47.2°~53.3°N， 118.8°~121.9°E的大兴安岭北部高纬度地区。不同气候情景下，樟子松对水热变化最为敏感，特别是在年平均气温升高5.0 ℃左右时将导致其在中国的分布区发生显著变化。

通过樟子松不同适生区的划分，根据各自区划气候特点提出以下保护对策：①在不同地区选取生长良好的樟子松单株进行种质资源监测、保护、采集和培育，培育最适合当地气候条件的优良品种。②未来在不适合樟子松生长地区引种时，应避免引种樟子松。③提前进行系统的种源区划研究，培育出不同特性的樟子松并在本地区推广。应该在年平均气温为−3.7~1.6 ℃，降水量为375~425 mm的东北等高纬度地区引种樟子松，以发挥其最大生态效益。