西南地区松属乔木对气候变化的响应

陆双飞，陈禹衡，周斯怡，殷晓洁

(1.西南林业大学林学院，云南 昆明 650224；2.南京林业大学生物与环境学院，江苏 南京 210037)

热量和水分是决定区域尺度上植物地理分布的主要环境因素[1]。以温暖化为主要特征的气候变化，将导致热量和水分在数量和空间上的变化，进而引起植物分布的变化。这一变化过程可能发生种群的局部灭绝或高度分散[2]。政府间气候变化专门委员会(Intergovernmental Panel on Climate Change，IPCC)第五次评估报告指出，相对于1986—2005年，2081—2100年的地球平均温度将升高0.3～4.8 ℃[3]。已有研究表明，未来气候下物种分布区可能发生不同程度的收缩或扩张[4-5]。

松属(Pinus)是松科(Pinaceae)中最大的属，全球100余种。我国松属资源丰富，自然分布22种10变种，引种栽培16种2变种[6]。松属植物具有极高的生态价值[7]、经济价值和医药价值[8]。中国西南是全球25个生物多样性热点地区之一[9]，也是中国裸子植物分布的中心[10]。然而，由于特殊的地理位置和地形地貌，使西南地区成为对气候变化最为敏感的地区之一[11]。因此，对该地区松属乔木对气候变化响应的研究意义重大。

采用最大熵模型(maximum entropy model，MaxEnt model)对西南地区12种松属乔木当前和未来气候情景下的潜在地理分布进行模拟，以期回答：①影响12种松属乔木分布的主导环境因子是什么？②当前和未来气候情景下12种松属乔木的潜在分布格局怎样？③潜在分布区的环境特征如何？在区域尺度上，为气候变化下西南地区松属乔木资源的保护和可持续经营提供理论参考。

1 材料与方法

1.1 研究区概况

西南地区包括了以云贵高原为中心的云南省、贵州省、重庆市、四川省及西藏自治区的一部分，总面积113.74万km2，占我国国土面积的11.84%，境内海拔高差较大，地形复杂多样，因此有着极其丰富的针叶林资源[13]。该地区有着众多我国特有的松属乔木如云南松(PinusyunnanensisFranch.)、细叶云南松(Pinusyunnanensisvar.tenuifoliaCheng et Law)、思茅松[Pinuskesiya(A.Chev.) Gaussen.]、毛枝五针松(PinuswangiiHu et Cheng)等树种，同时也有巴山松[Pinustabuliformisvar.henryi(Masters) C. T. Kuan]、华南五针松(PinuskwangtungensisChun ex Tsiang)、乔松(PinuswallichianaA. B. Jackson)等南方特有树种及华山松(PinusarmandiiFranch.)、白皮松(PinusbungeanaZucc. et Endi)等全国广布树种。但由于对森林长期的过度采伐、不合理的耕作方式和矿产资源开发等导致该地区的森林资源破坏严重，大面积的以松属乔木为建群种的林分遭到了破坏[14]。因此确定西南地区松属物种的适宜区对当地的天然林保护工程建设有着重要意义。

1.2 松属物种分布数据

松属树种的地理分布信息主要通过4个途经获取：①中国数字植物标本馆(http:∥www.cvh.ac.cn/)，②国家标本资源共享平台(http:∥www.nsii.org.cn/)，③全球生物多样性信息网络平台(http:∥www.gbif.org)，④中国期刊全文数据库(https:∥www.cnki.net/)。对于缺少准确坐标信息的分布记录，需精确到县级以下，通过百度坐标拾取系统(http:∥www.gpsspg.com/maps.htm)确定经纬度信息，并通过使用Google Earth(http:∥earth.google.com)卫星地图依据植物分布的具体生境来检验信息的准确性。每个1 km×1 km栅格中只随机保留1个分布点[12]。最终整理得到共计1 329条有效数据(图1)，将其汇总保存为CSV格式用于MaxEnt进行物种分布模型的建立。

1.3 环境数据及基础地理数据

气候数据来源于世界气候数据库(http:∥www.worldclim.org)，空间分辨率为30″(约1 km×1 km)。采用当前气候(1950—2000年)及未来气候情景2070年(2060—2080年)数据。未来气候数据由中国国家气候中心开发的BCC_CSM1-1模型产生，其对温度和降水等各变量的预测有较高精度[15]。未来气候情景选择代表性浓度路径4.5(representative concentration pathway 4.5，RCP4.5)，该模式情景与目前减排趋势相对应[16]。土壤数据来源于中国科学院寒区旱区科学数据中心(http:∥westdc.westgis.ac.cn)，其基于第二次全国土地调查南京土壤所提供的1 ∶100万土壤数据制作，栅格大小约为1 km2。数字高程模型(digital elevation model，DEM)数据来源于中国科学院计算机网络信息中心国际科学数据镜像网站(http:∥www.gscloud.cn)，空间分辨率为90 m，利用ArcGIS的表面分析工具提取得到海拔、坡度、坡向数据。利用ArcGIS将环境因子数据的分辨率重采样为30"，数据格式转换为ASCII(表1)。中国矢量地图来源于国家基础地理信息中心(http:∥www.ngcc.cn/ngcc/)的公益性地图。

图1 松属乔木的样点分布图Figure 1 Sample points distribution of Pinus species

1.4 环境因子筛选

为避免模型的过度拟合，提高模拟精度和实际生产中的可用性，需对环境因子按类别(温度、温度变化、降水、降水变化、地形和土壤6类)进行筛选，将贡献率低和性质相近的环境因子剔除[12]。首先利用37个环境因子和8个松属乔木地理分布数据，通过MaxEnt模型进行初步模拟，筛选出贡献率大于1.0%的因子；利用ArcGIS的波段集统计(band collection statistics)工具，对环境变量进行Pearson相关检验分析(图2)[17]。若同类环境影响因子中有2个或以上的环境影响因子相关性系数|r|≥0.8，则保留贡献率最高的因子参与最终模拟，筛选结果如表1所示。

1.5 最大熵模型分析

最大熵模型以最大熵理论为基础[18]，该理论可解释为，一个物种在没有约束的情况下，会尽最大可能扩散，接近均匀分布[19]。物种分布点所蕴含的气候、海拔和土壤等约束信息降低了熵，最大熵模型选取熵最大的分布作为最优分布[20]。在众多物种分布模型算法中，MaxEnt因其相对于其他建模方法具有预测精度高，在物种分布点较少的情况下仍能得到较准确的结果等优点而被广泛应用于入侵生物管理、造林区规划、生物多样性保护等研究[21-22]。研究表明，对于只有物种“存在”数据的情况下，MaxEnt模型的预测能力优于GARP、ENFA、Bioclim等其他物种分布模型[23-24]。因此采用MaxEnt模型进行西南地区松属乔木地理分布模拟。

将12个松属物种的分布样点数据与筛选出的12个环境影响因子数据导入MaxEnt模型进行物种地理分布模拟。随机选择75%的分布数据作为训练集，剩余25%的分布数据作为测试集，重复运算10次取平均值。若为未来气候下的潜在分布预测，需在“Projection layers directory/file”中添加未来某时段相对应的环境因子数据，以获得构建的当前关系模型在未来气候情景下的投影结果，即未来地理分布。使用受试者工作特征曲线(receiver operating characteristics curve，ROC curve)下的面积(the area under the ROC curve，AUC)进行模型的模拟精度检验[25]。AUC值的范围为[0.5，1]，其值越大表示模型准确性越高，当AUC值大于0.8时，表明模型的预测精度“很准确”[26]。根据MaxEnt模型输出的环境因子百分比贡献率评价各环境因子对模型的相对贡献，筛选影响松属乔木地理分布的主导环境因子。将累积贡献率超过90%的参与累积的因子作为主导因子，通过主导因子的响应曲线分析西南地区松属乔木潜在分布区的环境特征。

最大熵模型以ASCII格式输出物种在每个栅格内的存在概率(P)，取值范围为[0，1]。采用最大的训练敏感性和特异性(maximum training sensitivity plus specificity，MaxSS)方法确定阈值(y)，当植物在某一地区的存在概率P

表1 环境因子Table 1 Environmental factors

图2 环境因子相关性系数矩阵Figure 2 Pearson correlation coefficients (r) among environmental factors

2 结果与分析

2.1 模型准确性评价

西南地区12种松属乔木的地理分布模拟结果AUC均值如表2所示。AUC均值最大的是思茅松(0.987，标准误差为0.003)，最小的是油松(0.896，标准误差为0.023)。模型的AUC均值都大于0.8，模拟精度均达到“很准确”的水平，且标准误差值较低，反映出模型在不同重复之间具有较好的稳定性。因此，MaxEnt模型可用于12种松属乔木物种的潜在分布预测研究。

表2 模型的AUC值和MaxSS逻辑阈值以及松属植物适宜区质心移动距离Table 2 Model AUC values, MaxSS logic thresholds, and migration distance with centroids of suitable regions of Pinus

2.2 不同气候情景下松属乔木的潜在地理分布

模型计算的西南地区12种松属乔木的MaxSS逻辑阈值如表2所示，利用ArcGIS重分类工具根据适宜区划分标准进行潜在分布的可视化，得到当前(2000)和未来(2070 RCP4.5)松属乔木在中国境内的适生性分区(图3～图5)。

注：A表示当前气候情景，B表示2070 RCP4.5气候情景，C表示适宜区质心变化，S1表示低适宜区面积(104 km2)，S2表示高适宜区面积(104 km2)。Note: A refers to current climate scenario, B to the 2070 RCP4.5 climate scenario, C to the changes of centroids of suitable regions, S1 to the area of low suitable regions (104 km2), and S2 to the area of highly suitable regions (104 km2). 图3 细叶云南松等4种松属乔木的潜在适宜区及其质心Figure 3 Potential suitable regions and centroids of 4 Pinus species including Pinus yunnanensis var. tenuifolia

当前气候下的潜在地理分布范围方面。细叶云南松的分布范围最窄，主要集中在云南、贵州和广西3省交界地区。其低适宜区面积为15.73×104km2，高适宜区面积为2.56×104km2(图3)。马尾松的分布范围最广，分布于秦岭-淮河以南的四川东部、贵州、重庆、广西、湖南、湖北、广东、江西、安徽、福建及浙江等地。其低适宜区面积为138.45×104km2，高适宜区位于重庆东部、湖北西部和贵州东北部地区，面积为9.07×104km2(图5)。思茅松适宜区主要集中于除滇西北、滇东北以外的云南大部和海南省，高适宜区仅云南普洱和西双版纳地区(图3)。乔松主要分布于西藏南部、云南西北部以及四川西部和南部地区(图3)。云南松、高山松和华山松适宜区范围较为相似，主要分布于四川盆地周边以及云南大部，西藏东南部也有分布(图3和图4)。巴山松主要分布于大巴山东段和巫山地区，高适宜区集中于重庆、湖北和贵州交界地区(图4)。毛枝五针松分布于云南南部和东南部、贵州南部以及南岭山脉以南和武夷山以东地区，高适宜区集中于云南东南部的文山州，仅占总适宜区面积的4%(图4)。华南五针松潜在地理分布范围包括重庆、贵州、江西、福建、江苏、华南3省以及湖南西部和南部地区，高适宜区位于广西和广东两省北部。白皮松和油松的潜在分布区纬度较高，呈西南-东北走向，沿四川东北部、重庆、贵州北部、甘肃南部、陕西、山西和河北北部一直延伸辽宁南部，其中油松分布较广，在西藏东南部等地也有分布(图5)。

注：A表示当前气候情景，B表示2070 RCP4.5气候情景，C表示适宜区质心变化，S1表示低适宜区面积(104 km2)，S2表示高适宜区面积(104 km2)。Note: A refers to current climate scenario, B to the 2070 RCP4.5 climate scenario, C to the changes of centroids of suitable regions, S1 to the area of low suitable regions (104 km2), and S2 to the area of highly suitable regions (104 km2). 图4 巴山松等4种松属乔木的潜在适宜区及其质心Figure 4 Potential suitable regions and centroids of 4 Pinus species including Pinus tabuliformis var. henryi

未来气候下潜在地理分布面积变化。2070 RCP4.5气候情景下，乔松、思茅松、细叶云南松3个树种的适宜区面积增加最大，分别增加了87.06%、54.23%、47.18%(图3)。其中，思茅松的高适宜区面积(S2)增长比例较大，增长了137.32%。特别地，细叶云南松虽然适宜区总体呈现扩张趋势，但高适宜区却呈现收缩趋势，由当前的2.56×104km2减少为2070 RCP4.5气候情景下的2.04×104km2，减少了20.31%。云南松、高山松和华山松3个树种的适宜区面积增加也较大，分别增加了29.50%、26.21%和14.25%(图3和图4)。其中，云南松的高适宜区面积增加最大，由当前的5.84×104km2增加至2070 RCP4.5气候情景下的16.46×104km2，增加了181.85%(图3)。油松的适宜区扩张范围最小，适宜区总面积增加了1.27%，但高适宜区面积减少了6.54%。2070 RCP4.5气候情景下，马尾松、华南五针松、白皮松、毛枝五针松和巴山松5个树种的适宜区出现了收缩(图4和图5)。其中，巴山松的适宜区范围收缩最大，低适宜区和高适宜区面积分别减少了59.36%和67.47%(图4)。马尾松的低适宜区面积增加了1.83%，但高适宜区面积却下降了86%，而白皮松低适宜区面积减少了26.14%，但高适宜区面积增加了48.95 %(图5)。

注：A表示当前气候情景，B表示2070 RCP4.5气候情景，C表示适宜区质心变化，S1表示低适宜区面积(104 km2)，S2表示高适宜区面积(104 km2)。Note: A refers to current climate scenario, B to the 2070 RCP4.5 climate scenario, C to the changes of centroids of suitable regions, S1 to the area of low suitable regions (104 km2), and S2 to the area of highly suitable regions (104 km2). 图5 华南五针松等4种松属乔木的潜在适宜区及其质心Figure 5 Potential suitable regions and centroids of 4 Pinus species including Pinus kwangtungensis

未来气候下潜在地理分布空间变化。2070 RCP4.5气候情景下，西南地区12种松属乔木不同等级适宜区的质心沿不同方向发生变化，且移动距离差异较大(图3～图5)。除思茅松和毛枝五针松外，整体上呈现出向西偏北或东偏北的高海拔高纬度地区移动。思茅松的高适宜区质心由于2070 RCP4.5气候情景下高适宜区面积的增加而向西移动128 km(表2，图3)，毛枝五针松的高适宜区质心则由于云南境内适宜区的锐减而向东南方向移动759 km(表2，图4)。

2.3 松属乔木潜在分布区的环境特征

MaxEnt模型输出的各环境因子对西南地区12种松属乔木潜在分布的百分比贡献率(percent contribution，PC)如图6所示。年降水量(BIO12，PC=27.61%)>最冷月最低气温(BIO6，PC=18.42%)>气温的季节性(BIO4，PC=14.81%)>海拔(ALT，PC=8.26%)>等温性(BIO3，PC=7.20%)>土壤质地(T_USDA_TEX，PC=5.64%)>气温日较差(BIO2，PC=4.38%)>降水的季节性(BIO15，PC=4.09%)8个环境因子的平均百分比贡献率之和为90.41%，为影响西南地区12种松属乔木潜在地理分布的主导环境因子。

注：各环境因子含义见表1，各松属乔木代码含义见表2。Note：the meaning of each environmental factor is shown in Table 1 and the meaning of each Pinus species code is shown in Table 2. 图6 西南地区松属乔木环境因子的百分比贡献率Figure 6 Percentage contributions of environmental factors for the Pinus species in southwestern China

不同物种对环境因子的要求存在差异，根据环境因子与物种存在概率的关系，得到主导环境因子在松属乔木适宜区的变化范围如表3所示。影响西南地区12个松属乔木分布的主导环境因子年降水量(BIO12)在适宜区的平均变化范围为647～2 600 mm，最冷月最低气温(BIO6)平均变化范围为-8～10 ℃，气温的季节性(BIO4)的范围为358～756，海拔(ALT)平均变化范围为281～3 054 m，等温性(BIO3)平均变化范围为25～45，气温日较差(BIO2)平均变化范围为6～11 ℃，降水的季节性(BIO15)平均变化范围为43～99。其中，分类型环境因子土壤质地(T_USDA_TEX)主要为3(轻黏土)、7(粉壤土)、10(砂黏壤土)。

表3 西南地区松属乔木主导环境因子变化范围Table 3 Variation range of dominant environmental factors for the Pinus species in southwestern China

表4 土壤质地对照表Table 4 Soil texture comparison table

将百分比贡献率最高的环境因子作为最主要环境因子，其与松属乔木存在概率之间的响应曲线如图7所示。乔松的存在概率与等温性(BIO3)呈负偏态关系，峰值约为41；思茅松的存在概率与气温的季节性(BIO4)呈负相关关系，分布上限为300。而云南松的存在概率与气温的季节性则近似呈正偏态关系，峰值为450，上限为680。白皮松、马尾松和巴山松的存在概率与最冷月最低气温(BIO6)近似呈正态分布，峰值分别为-6、0和-4 ℃；年降水量(BIO12)与华山松、华南五针松、油松、油松、毛枝五针松和细叶云南松的存在概率近似呈现正偏态关系，峰值分别为900、1 600、510、1 400和1 250 mm；高山松的存在概率与海拔(ALT)近似呈正态分布关系，最适宜的海拔段为2 400～3 400 m。

图7 西南地区松属乔木最主要环境因子的响应曲线Figure 7 Response curves of the most important environmental factors for the Pinus species in southwestern China

3 讨论与结论

采用最大熵模型对西南地区12种松属乔木当前和未来气候情景下的潜在地理分布进行模拟。研究结果表明，乔松、思茅松、细叶云南松、云南松、高山松、华山松和油松的适宜区总体呈现扩张趋势，而细叶云南松和油松高适宜区面积则有所减小；而巴山松、毛枝五针松、白皮松、华南五针松和马尾松的适宜区呈缩小趋势。除思茅松和毛枝五针松两个分布纬度较低的树种在未来气候下，适宜区质心分别向西南和东南方移动，其余树种的适宜区质心大致向东北或西北方向移动。西南松属植物适宜区年降水量的平均变化范围为647～2 600 mm，最冷月最低气温的范围为-8～10 ℃，气温的季节性为358～756，海拔均值为281～3 054 m，等温性范围为25～45，气温日较差范围为6～11 ℃，降水的季节性范围为43～99。土壤质地主要为轻黏土、粉壤土和砂黏壤土。该结果对松属植物资源的保护、造林区规划以及可持续经营等具有一定指导意义。

3.1 松属乔木潜在分布区变化的原因

植物地理分布是多种因素综合作用的结果，各类环境因子对不同物种的作用不同，每一物种均有决定其物种分布范围的主导因子，随着主导因子的变化，物种的分布面积将会随之变化。确定主导环境因子有助于定量描述环境与目标物种地理分布之间的关系，在简化关系的同时不至于丢失过多有价值的信息[29]。本研究根据12个松属树种的平均百分比贡献率，筛选出8个主导环境因子，其中地形和土壤则被认为是影响物种丰富度小尺度分布格局的重要因子[30-31]。

气候被认为是影响植物分布大尺度格局的主要因素[32]。气候类因子对西南地区松属乔木潜在分布的累计贡献率最高，达76.51%，这与有关西南常见针叶树分布研究的结果相符[16]。气候变化对物种分布甚至存亡具有重要的影响，未来气候下物种分布范围的收缩和扩张已被大量研究证实[33]，关于针叶树种的研究显示物种分布范围变化方向和程度不尽相同[5-16]。随着全球气候变暖与环境更加湿润，所有物种的未来适宜区质心均有向北移动的趋势。巴山松、毛枝五针松、华南五针松和马尾松的适宜区与高适宜区面积的缩小可能与其对气温日较差较为敏感有关，未来气候情景中气温日较差的变化可能会通过改变植物内源激素的分泌从而影响植物繁殖，进而改变植物的适宜区[34]。乔松、思茅松、细叶云南松、云南松、高山松、华山松和油松的适宜区总体呈现扩张趋势则可能与它们在降水量季节变异系数与海拔上有着较强的适应能力有关，未来降水量的季节性变化可能会深刻改变各个地区的地下水径流与土壤含水量，进而对更为耐旱的针叶树种提供了机遇[35]。同时更大的海拔适宜范围意味着这些树种得以在云贵高原有着更大的发展空间。

3.2 对后续研究的展望

松属乔木潜在分布模拟的准确性受物种样点数据的采样精度和完整程度、环境变量的选取以及模型参数设置等因素的影响。研究所用物种的地理分布数据主要来源于数字标本馆，由于受到地形和交通等条件的限制，采样点多集中于交通便利的道路附近、自然保护区等特殊区域[29]。这些采样点通常不代表环境梯度的随机样本，导致许多物种的实际分布信息不完整，可能影响模拟结果[36]。巴山松的样点数最少(19个)，而AUC值达0.970，印证了MaxEnt模型受样本容量影响小，在分布点数据较少的情况下也能得到满意的预测结果[37]。马尾松的AUC值较低，为0.914，而其分布范围最广，这与YANGetal[38]的研究相符，说明对于分布较广的物种建立高精度关系模型比较困难。有研究表明土地利用类型对植物分布的影响也十分显著[19]，后续的研究可综合考虑此类对目标物种分布有影响的生物因素数据集，以期获得更加准确的预测结果，同时对于分布范围较广的物种，后续研究可以采用多模型混合应用的方式以提高预测准确度。