紫果云杉(Picea purpurea)是我国青藏高原东北部特有的亚高山暗针叶林建群树种，不仅具有生长快、材质好、耐寒耐贫瘠等优点，也是长江及黄河上游地区重要的水源涵养林和水土保持林优势树种，在西部生态安全屏障和林业生产实践中具有重大意义(中国森林编辑委员会，1998)。基于细胞质DNA片段、核基因序列以及花粉性状的研究表明紫果云杉起源于丽江云杉(P. likiangensis)和青杄(P. wilsonii)的自然杂交作用，是裸子植物中稀有的同倍体杂交成种案例(孙永帅，2012)。生态位分化在同倍体杂交物种形成过程中起到重要作用，大多数杂交物种与亲本种间均具有明显的生境分化，且前者往往占据更极端的生境(Gompert et al., 2006; Schumer et al., 2014; Yakimowski et al., 2014)。例如:向日葵属(Helianthus)中3个同倍体杂交种H. anomalus、H. deserticola和H. paradoxus分别占据沙地、沙漠和盐沼生境(Rieseberg et al., 2003); 高山松(Pinus densata)占据2个亲本种都不能正常生长的高海拔地带(Wang et al., 2011)。这种对新生境的占领使杂交物种避免了与亲本种可能的渐渗杂交和对可利用资源的竞争，同时对极端环境的适应则进一步塑造了杂交物种的生物学特征，促进了其与亲本生态位的分化(Gompert et al., 2006; Rieseberg，1997)。与其他同倍体杂交物种类似，紫果云杉与亲本种在地理分布上也存在明显分化，即:紫果云杉主要分布于青藏高原东北缘的四川西北部、甘肃南部和青海东南部的中纬度高海拔(2 600~3 800 m)地带，而丽江云杉主要分布于青藏高原东南缘的云南西北部和四川西南部的低纬度高海拔(2 300~3 800 m)地带，青杄则广布于我国中西部地区的中纬度低海拔(1 400~2 800 m)地带(郑万钧，1938)。尽管紫果云杉沿海拔和纬度与亲本种间存在着明显的生境分化，但决定其生态位分化的环境因子目前并不清楚，鉴别和分析杂交物种独特的生态位特征有助于阐明其占领极端新生境以及物种形成的生理生态基础(Mallet，2005)。

气候因子是同倍体杂交物种生态位分化的重要组成部分(Rundle et al., 2005)。Mao等(2011)结合物种分布数据、GIS技术和生态模型分析对同倍体杂交物种高山松的研究表明，其生境具有较低的生长积温和较高的冷冻频率等，揭示了生态选择在高山松物种形成中的重要作用。而近百年来，随着温室气体排放量的迅速增加，以气候变暖为主要特征的全球气候发生着剧烈变化，而这种变化对物种的分布将产生深远的影响(Forster et al., 2007; Pearson et al., 2003)。气候变暖对全球中低纬度高山地区和高纬度地区的暗针叶林将造成巨大威胁，如: Mantgem等(2009)发现区域性增温及其造成的水分匮缺可能是美国西部松属(Pinus)、冷杉属(Abies)、铁杉属(Tsuga)等森林树种死亡率在近几十年呈显著增加的主要原因; Garzón等(2008)研究表明，到21世纪80年代欧洲赤松(Pinus sylvestris)和欧洲冷杉(Abies alba)等树种的潜在分布将出现大幅下降。紫果云杉分布区不仅是我国生态脆弱区亦是气候变化敏感区，预测其潜在分布区的变化对未来如何应对和减缓气候变暖对其生态屏障功能产生的影响具有重大现实意义。鉴于此，紫果云杉是研究未来气候变化对我国西部高山地区(青藏高原东北部)森林分布影响的理想材料，本文研究目的为：1)定量分析紫果云杉与亲本种间生境差异，以期揭示其生态位分化的主导环境因子; 2)预测云杉属3个树种在未来气候变化背景下潜在分布区的变化，深化未来气候变暖对其潜在分布区变化影响的认识。

以《中国植物志》中所述云杉属3个树种的分布范围为据，结合前人研究以及中国数字植物标本馆(CVH; http://www.cvh.org.cn/englishindex.asp)的分布点，筛去重复和不确定点，最终确定可代表 3个树种核心分布的284个样点，包括58个紫果云杉样点、127个青杄样点和99个丽江云杉样点(图 1)。

图 1 树种分布样点 Figure 1 Species distribution of the occurrence sites

为全面准确地量化云杉属3个树种的环境生态位特征，初始收集了26个环境变量，其中19个气候因子变量来自Worldclim网站(http://www.worldclim.org/)，3个土壤特性相关变量(pH、湿度和有机碳含量)来自可持续发展和全球环境研究中心(Center for Sustainability and the Global Environment)，其他4个环境变量(大气压强、生长积温、湿天频率和地面结霜率)来自政府间气候变化专门委员会数据库(Intergovermental Panel on Climate Change database)。为减少各变量间的自相关性干扰，最后选取了9个具有低相关性(种内种间Pearson相关系数r＜0.7)的环境变量用于后续分析(表 1)。

表 1 本研究所用环境变量名称及来源 Tab.1 The environmental variables and their sources used in this study

表 1 本研究所用环境变量名称及来源 Tab.1 The environmental variables and their sources used in this study 环境变量 Environmental variables 简写 Abbreviations 空间分辨率 Spatial resolutions 来源 Sources 最暖月最高温Max temperature of warmest month / ℃ bio5 2.5 min WORLDCLIM(http://www.worldclim.org/) 最冷月最低温Min temperature of coldest month / ℃ bio6 2.5 min WORLDCLIM(http://www.worldclim.org/) 降雨季节性变化Precipitation seasonality(%) bio15 2.5 min WORLDCLIM(http://www.worldclim.org/) 最暖季降雨量Precipitation of warmest quarter / mm bio18 2.5 min WORLDCLIM(http://www.worldclim.org/) 最冷季降雨量Precipitation of coldest quarter / mm bio19 2.5 min WORLDCLIM(http://www.worldclim.org/) 土壤有机碳密度Soil organic carbon density /(kg·m-2) SOC 5 min The Nelson Institute Center for Sustainability and the Global Environment, University of Wisconsin-Madison(http://www.sage.wisc.edu/atlas/index.php) 土壤pH Soil pH SpH 5 min The Nelson Institute Center for Sustainability and the Global Environment, University of Wisconsin-Madison(http://www.sage.wisc.edu/atlas/index.php) 土壤水分Soil moisture / mm SM 5 min The Nelson Institute Center for Sustainability and the Global Environment, University of Wisconsin-Madison(http://www.sage.wisc.edu/atlas/index.php) 地面结霜率 Ground-frost frequency/ (d·month-1) FRS 5 min Intergovermental Panel on Climate Change database(http://www.ipcc-data.org/)

首先通过非参数Kruskal-Wallis多重秩和检验比较了各环境变量在物种水平的分化变异，并绘制了各物种对应各环境变量的核密度图；其次，通过判别式分析(DFA)评估了各环境变量对两两物种对分化的贡献度，用Wilks’s λ检验物种对无差异的零假设；最后，通过主成分分析(PCA)综合评估了3个物种的环境变量差异，每个环境变量对每个物种空间分布的相对贡献通过绘制距离载荷图(biplot)表示。核密度图、非参数检验、PCA分析和距离载荷图(biplot)利用R3.0.3(R Development Core Team，2014)完成，DFA分析利用SPSS 16.0完成。

参照Ren等(2016)的方法，从CMIP5计划中选取了对中国温度和降雨模拟效果较好的3种大气环流模型(general circulation models，GCMs: BCC-CSM1-1、CCCma_CanESM2和CSIRO-Mk3.6.0)和3种最新的代表性浓度路径情景(representative concentration pathways，RCPs: RCP2.6、RCP4.5和RCP8.5)。其中RCP2.6代表缓和场景, 指采用缓和策略如提高能源使用效率、使用可再生型新能源等方法使温室气体排放量大幅降低; RCP4.5代表维持温室气体排放中等平稳场景; RCP8.5代表温室气体高度排放场景。3种情景下, 2100年中国平均地表温度增幅分别为1.87、2.88、5.51 ℃(Meinshausen et al., 2011; 梁玉莲等，2016)。从CCAFS网站(http://www.ccafs-climate.org)获取了以上模型及情景下未来2个时间段，即21世纪50年代(2050s)和80年代(2080s)对应的共18套模拟数据。从Worldclim网站下载了用于构建物种当前潜在分布区的环境气候数据。以上数据均采用2.5 min的空间分辨率。此外，为减少变量间的共线性对模型精确度的影响，选取了种内种间Pearson相关系数r＜0.7的环境变量。

云杉属3个树种潜在分布区模型的构建通过软件MaxEnt v3.3k来完成，该软件基于最大熵模型(MaxEnt)对物种分布有良好的预测，是目前认可度高且应用较为广泛的物种分布模拟软件之一(殷晓洁等，2013；Phillips et al., 2008)。将上述筛选后变量及各物种分布点导入MaxEnt，使用默认参数设置(Phillips et al., 2008)。为评估模型的有效性，设置了75%的物种分布点作为训练数据集，剩余25%作为交互检验数据集，模型精度通过受试者工作特征曲线(receiver operating characteristic curve, ROC)进行评估，以ROC曲线下面积(area under the curve, AUC)表示，其值常介于0 ~ 1.0，理论上，AUC值为0.5 ~ 0.7时模型预测效果较差，0.7~0.9时为一般，高于0.9时预测效果较好。

为研究每个树种潜在分布区的变化，本研究使用第10分位训练存在阈值(10^th percentile training presence threshold)将具有适生概率等级划分的分布图转换为二元分布图像(适生/非适生)。对应3种GCMs，每个物种每个时间段在每个情景下均得到3张二元分布图，通过Arcgis 10.3的图像叠加，将其合为1张，用于后续分析和潜在分布区面积计算。以上图像再经图像挖填处理，以对比各物种在各情景下未来与当前分布区的变化。

Kruskal-Wallis多重秩和检验结果表明，紫果云杉有2个环境变量比亲本更极端，1个介于亲本之间，剩余6个变量表现为与其中1个亲本相近(图 2)：在热能状态方面，尽管最暖月最高温(bio5)与丽江云杉相近，而紫果云杉的最冷月最低温(bio6)显著低于其亲本种，其生境总体仍表现出比亲本种更低的热能特征; 在水分状态方面，紫果云杉的生境总体降雨适中水分充足，主要表现为其最冷季降雨量(bio19)与青杄相近高于丽江云杉，而最暖季降雨量(bio18)高于青杄而低于丽江云杉，且具有高于亲本的土壤湿度(SM)。

图 2 紫果云杉与亲本种9个环境变量的核密度分布 Figure 2 Kernel density plots for 9 environmental parameters of Picea purpurea and its parent species “=”表示2个物种间无差异(P=0.05)，“＞”或“＜”表示物种间存在显著性差异以及值的高低(P＜0.05)。图中bio5，bio6和FRS为实际数据×10。紫果云杉、丽江云杉和青杄在图中分别用p，l，w表示。 The equal sign means there is no significant difference(at the level of P=0.05)between two species, while the significant difference are coded as higher or lower(P < 0.05). Besides, the data of bio5, bio6 and FRS showed in this figure are the actual data ×10. Picea purpurea, P. wilsonii, and P. likangensis are donated p, w, and l, respectively.

判别式分析(DFA)结果(表 2)表明，所有树种对间均有显著的环境分化(Wilks’s λ，P＜0.001)。bio5、地面结霜率(FRS)、土壤有机碳含量(SOC)和SM对紫果云杉与青杄间环境变量分化贡献率较高，而bio6和降雨量季节性变化(bio15)对紫果云杉与丽江云杉环境变量分化贡献率较高。另外，bio5、SM和bio18对丽江云杉与青杄间环境变量分化贡献率较高。综合两亲本的环境数据后发现，bio5、bio6、和SM对紫果云杉与亲本种环境变量分化的贡献率较高。

表 2 9个环境变量分别对云杉属两两树种判别(DFA)的相对贡献^① Tab.2 The contributions of environmental variables in discriminant function analysis(DFA) between pair-wise comparisons of Picea species

表 2 9个环境变量分别对云杉属两两树种判别(DFA)的相对贡献① Tab.2 The contributions of environmental variables in discriminant function analysis(DFA) between pair-wise comparisons of Picea species 环境变量 Environmental variables 紫果云杉-青杄 P. purpurea-P. wilsonii 紫果云杉-丽江云杉 P. purpurea-P. likiangensis 青杄-丽江云杉 P. wilsonii-P. likiangensis 紫果云杉-青杄和丽江云杉 P. purpurea-P. likiangensis & P. wilsonii bio5 -0.577 0.037 -0.291 0.344 bio6 -0.195 0.273 0.055 0.343 bio15 -0.062 0.456 0.134 0.232 bio18 0.100 0.209 0.222 0.096 bio19 -0.100 -0.089 -0.092 0.051 FRS 0.339 -0.023 0.152 -0.236 SOC 0.339 -0.058 0.139 -0.288 SM 0.430 -0.073 0.227 -0.339 SpH -0.143 -0.052 -0.133 0.062 DFA(Wilks’s λ) 0.211 0.101 0.069 0.401 DFA(P) ＜0.000 ＜0.000 ＜0.000 ＜0.000 ①每组比较中贡献较高的环境变量以加粗示出。Values in boldface are for the most significant variables for each comparison.

基于9个环境变量的PCA分析共提取出3个主成分(特征值>1)，可解释77.5%的分布点变异，且主成分1和2(PC1和PC2)可解释65%的变异。PC1和PC2可将284个分布点清晰地分为3部分(图 3)，分别对应紫果云杉、丽江云杉和青杄3个树种，且每个树种均表现出不同的生态位环境特征。紫果云杉占据一个以高FRS、SOC和SM及低bio6为主要特征的潮湿低温生境，而青杄以高bio5、bio19和低bio18、FRS和SM为主的夏季温暖干旱、冬季潮湿的生境，丽江云杉的生境具有温和潮湿(高bio6和bio18)和强季节性变化(高bio15)的特征。

图 3 基于284个物种分布点9个环境变量的主成分分析(PCA) Figure 3 Principal component analysis of distance biplot for 284 species occurrence data on 9 environmental variables

对应3种大气环流模型(GCMs)和3种代表性浓度路径情景(RCPs)以及2个时间段，每个树种通过模拟均得到18个未来潜在分布区模型。每个树种的平均模型测试AUC值均高于0.99(紫果云杉0.997、丽江云杉0.997和青木千0.992)，表明所有模型都具有极好的预测精度。对比潜在分布区面积变化(表 3，表 4)发现，紫果云杉仅在RCP2.6下2080s的潜在分布区面积与当前相比出现了缩减(总体减少约5%)，而在此情景下2050s的潜在分布区面积和其余2个RCPs下未来2个时间段的潜在分布区面积都高于当前，总面积平均增加了17%以上; 而青杄在所有情景下未来潜在分布区面积均低于当前，总面积平均减少了21%以上；丽江云杉除在RCP4.5和RCP8.5下2080s的潜在分布区面积稍高于当前(增加不到2%)外，其余情景下的适生面积也都低于当前，总面积平均减少5%左右。

表 3 云杉属3个树种当前及未来潜在分布区面积预测 Tab.3 Predicted potential distribution area of three Picea species

表 3 云杉属3个树种当前及未来潜在分布区面积预测 Tab.3 Predicted potential distribution area of three Picea species 104 km2 时间段 Periods RCPs 紫果云杉 P. purpurea 丽江云杉 P. likiangensis 青杄 P. wilsonii 当前适生区 Current 16.22 23.17 142.04 2050s RCP2.6 17.32 21.53 103.32 RCP4.5 18.16 20.42 119.49 RCP8.5 18.39 22.94 118.21 2080s RCP2.6 15.46 22.53 106.59 RCP4.5 24.91 23.56 114.13 RCP8.5 16.72 23.25 105.60

表 4 云杉属3个树种未来潜在分布区与当前相比面积的变化 Tab.4 The change in suitable environmental area for three Picea species in the future compared with the current

表 4 云杉属3个树种未来潜在分布区与当前相比面积的变化 Tab.4 The change in suitable environmental area for three Picea species in the future compared with the current % 时间段 Periods RCPs 紫果云杉 P. purpurea 丽江云杉 P. likiangensis 青杄 P. wilsonii 增加区 Gained 减少区 Lost 增加区 Gained 减少区 Lost 增加区 Gained 减少区 Lost 2050s RCP2.6 37.21 31.70 16.83 23.93 13.30 38.01 RCP4.5 31.15 19.09 14.69 26.45 15.65 31.26 RCP8.5 32.21 19.07 20.87 21.94 14.31 30.97 2080s RCP2.6 16.48 21.19 18.52 21.27 6.66 33.68 RCP4.5 69.52 15.97 24.21 22.29 9.47 28.86 RCP8.5 26.20 23.61 21.95 21.92 15.62 41.33

与当前分布区相比，RCP2.6下紫果云杉未来在青藏高原东北部的潜在分布区减少，而增加了青藏高原东南部边缘的潜在分布; 而丽江云杉减少了其在青藏高原东南部的分布，青杄则在其原有分布区基础上减少了其在四川南部、内蒙古东部等边缘的分布(图 4)。在RCP4.5和RCP8.5下，紫果云杉潜在分布区由青藏高原东南部边缘向内部扩张的趋势更为明显; 而丽江云杉在青藏高原东南部的潜在分布区有轻微减少，但增加了青藏高原东部和南部的分布; 青杄潜在分布区的变化与RCP2.6下一致(图 5，图 6)。

图 4 在温室气体最低排放情景(RCP2.6)下，云杉属3个树种当前及未来的潜在分布区 Figure 4 Current and future potential distributions of three Picea species under the lowest greenhouse gas emission scenario (RCP2.6) A1-A3:紫果云杉; B1-B3:丽江云杉; C1-C3:青杄。所有未来潜在分布区分布的预测均基于3个大气环流模型(GCMs: BCC-CSM1-1、CCCma_CanESM2和CSIRO-Mk3.6.0)。下同。 A1-A3: Picea purpurea; B1-B3: P. likiangensis; C1-C3: P. wilsonii. All the future distribution projections are based on an ensemble of predictions from three general circulation models(GCMs: BCC-CSM1-1, CCCma_CanESM2 and CSIRO-Mk3.6.0). The same below.

图 5 在温室气体中度稳定排放情景(RCP4.5)下，云杉属3个树种当前及未来的潜在分布区 Figure 5 Current and future potential distributions of three Picea species under the moderate greenhouse gas emission scenario (RCP4.5)

图 6 在温室气体高度排放情景(RCP8.5)下，云杉属3个树种当前及未来的潜在分布区 Figure 6 Current and future potential distributions of three Picea species under the highest greenhouse gas emission scenario (RCP8.5)

同倍体杂交物种生态位分化的产生是初始杂交后代中出现超亲分离和随后生态选择共同作用的结果(Yakimowski et al., 2014)。环境对物种的生态选择作用可表现在多个方面，如生境、气候、资源等(Rundle et al., 2005)，因此为评估同倍体杂交物种形成过程中生态选择的作用，需对杂交物种及其亲本种主要生态因子的差异做详尽的量化分析(Mao et al., 2011)。基于前人对紫果云杉生态位分化的研究(孙永帅，2012)，本研究通过对紫果云杉生境的量化分析，明确了冬季低温和高土壤湿度是主导其与亲本种生态位分化的关键环境因子。

同倍体杂交物种与亲本种产生生态位分化并占领极端新生境，与初始杂交居群超亲分离产生超亲性状密切关联(Rieseberg et al., 1999)。超亲分离现象在分离的杂交居群中普遍存在，如向日葵属同倍体杂交物种Helianthus anomalus叶片肉质化程度、瘦果长度、表皮导度等一系列与耐旱相关的形态、生理性状表现出超亲分离特征(Schwarzbach et al., 2001); 同倍体杂交物种H. paradoxus具有与耐盐相关的超亲性状如叶片肉质化、元素摄入等(Lexer et al., 2003)。此外，高山松作为木本植物同倍体杂交物种的典型案例也表现出如高的总干物质积累量、长期水分利用效率以及低温胁迫下维持高光合的能力等性状的超亲分离，这些超亲性状的出现促进了高山松对高海拔极端生境的占领(Ma et al., 2010; 2013)。基于核基因的种群遗传数据表明紫果云杉的物种形成始于约130万年前最大冰期时，随后在75万年前更新世冰期出现了大规模的居群扩张(Sun et al., 2014)。因此，推测在紫果云杉物种形成初期，一些杂交居群也可能产生了耐冬季低温相关的超亲性状，这些居群逐渐向青藏高原东北部高海拔低温区域扩张，随后在生态选择作用下与亲本种间逐步产生生态位分化。

在未来气候显著增温情景(RCP4.5和RCP8.5)下，紫果云杉与亲本种的潜在分布区变化并不一致：紫果云杉潜在分布区面积显著增加，而丽江云杉总体波动较小，但青杄分布区域则出现了大幅减少，此结果可能与不同区域对气候变化的响应程度不同有关。梁玉莲等(2016)对我国21世纪末期气候变化预估发现，在RCP4.5和RCP8.5下年降水和夏季降水增大值均出现在青藏高原东南部地区，而我国中西部地区夏季降水量减少趋势明显。进一步表明青杄分布区未来夏季干旱将加剧，造成了该物种未来潜在分布区的减少，而青藏高原东南部降雨量的增大则促进了紫果云杉的扩张及丽江云杉的迁移。此外，紫果云杉和丽江云杉在同一情景下出现了2050s和2080s潜在分布区面积变化趋势不一致，这可能与大气环流模型存在地区气候变幅差异有关(Xu et al., 2012; 2001)。

与其他气候变化对物种分布预测研究类似，本研究基于MaxEnt模型对云杉属3个树种未来潜在分布区的预测也有局限性，主要由于模型中未考虑种间相互作用，如竞争、病虫害侵染等对树种分布的影响(Guisan et al., 2005)。首先森林中不同树种间的竞争会影响林木生物量的分配，进而影响树木的生长、形态特征和存活状态，最终可能影响树种的分布(Szewczyk，2009); 再者病虫害也是影响森林存活的重要因素，尤其随着经济全球化，外来病虫入侵几率增大，对本土林木生长也会造成威胁(Gregory et al., 2010)。而且，一些不确定因素如火灾、人为干扰等对物种分布的影响无法用模型进行预测(张雷等，2011)。此外，模型所预测的当前物种分布仅指理论上物种可能分布的区域，并非物种的实际分布，如本文中青杄的潜在分布甚至包含了一些国外地区，事实上，青杄为我国特有树种，在国外是没有分布的。尽管如此，本研究仍为同倍体杂交物种紫果云杉与其亲本种在未来气候变化下潜在分布的变化提供了参考：假设未来气候变化与模型预测一致，气候变暖将有助于同倍体杂交物种紫果云杉分布范围的扩张，推测其在未来西部高山地区将发挥更高的经济价值和更重要的生态安全屏障作用。

紫果云杉与亲本种产生了明显的生态位分化，冬季低温和高土壤湿度是主导其生态位分化的关键环境因子; 在未来气候显著增温背景下，基于物种分布模型预测紫果云杉潜在分布区面积将显著增加，而其亲本种潜在分布区面积将显著减少或波动较小，推测其在未来将发挥更重要的生态安全屏障作用。