1. 引言
物种的地理分布和环境紧密相关，环境因子常常制约着物种的生存及分布[1]。气候变化是影响物种适生区分布的重要因素之一[2]。特别是在全球变暖的背景下，气候变化引起的极端天气，叠加人类活动等多种因素，导致物种分布范围呈现向两极或更高海拔迁移的趋势，物种的空间分布格局发生剧烈变化，同时生境破碎化加剧，严重影响了生物多样性[3]-[6]。
水晶兰(Monotropa uniflora L.)，杜鹃花科水晶兰属，是一种多年生草本腐生植物，又称梦兰花、水兰草、银锁匙，广泛分布于中国(青海、浙江、安徽等省区)、俄罗斯、日本、印度、东南亚以及北美地区[7]。水晶兰具有重要的药用价值，性平和，味微咸，全草及根均可入药。临床主要用于补虚止咳，治疗肺虚咳嗽[8]。然而，由于过度采挖及生境破坏，水晶兰的野生种群受到严重威胁，2013年水晶兰在中国被列入国家近危(NT)珍稀植物名录[9]。作为一种珍稀的腐生草本植物，水晶兰对生境要求严格且特殊，大多生长在海拔800米至3200米的黑暗潮湿的混交林或落叶阔叶林里，依赖特定的真菌共生体系生存[10]。这种对环境的高度敏感性，使其极易受到外界环境波动的影响。然而，目前水晶兰已有的研究主要集中在其生态特征、药用价值等方面，而针对其地理分布格局及其对环境变化响应的研究则较为有限。
物种分布模型(SDMs)是通过分析物种分布数据和相关环境因素，预测物种适宜栖息地分布的数学模型。这些模型通过分析从物种分布点收集的环境信息，评估物种的生存需求，并将这些数据投影到选定的研究区域[11]-[13]。目前，物种分布模型主要包括生物气候分析系统模型(BIOCLIM) [14]、生态位因子分析模型(ENFA) [15]、广义线性模型(GLM)和最大熵模型(MaxEnt) [16]等。与其他模型相比，MaxEnt模型不仅具有良好的预测性和稳定性，而且具有简单、易于操作和样本量小的优点，与ArcGIS空间分析相结合后，广泛应用于分析物种的潜在适生区[17]-[19]。目前，MaxEnt模型已经成为生态学领域中应用最广泛、效果最好的物种分布模型之一。例如益争祝玛等人利用MaxEnt模型预测乌桕在中国的潜在分布[20]，唐燕等人利用MaxEnt模型对中华枸杞在全国的适宜性分布概率进行预测[21]，陈旭等人利用MaxEnt模型预测气候变化对小花老鼠簕潜在适生区的影响[22]。然而，这些研究在环境因素选择上多局限于气候因素，本研究在此基础上，综合考虑气候、土壤、地形、人类活动等多重因素，通过MaxEnt模型对水晶兰的潜在适生区进行更为全面、准确的预测。但是，水晶兰的分布必然受到其共生真菌(如红菇科真菌)及这些真菌的宿主树种(如壳斗科、松科等植物)分布的限制，故后续还需结合共生真菌和宿主树种等因素进行更为精确的预测。
2. 材料和方法
2.1. 数据获取与筛选
2.1.1. 物种分布数据
水晶兰的物种分布数据来自于全球生物多样性信息网络平台(GBIF, https://www.gbif.org)、中国植物标本馆数字化平台(CVH, http://www.cvh.ac.cn)，整合后共获取水晶兰分布数据坐标300个。在消除错误和重复的记录后，我们在SDMtoolbox中使用“Spatially-Rarefy-Occurrence-Data”工具，以确保每个1 × 1公里的网格仅包含一条出现记录，这种方法最大限度地减少现有空间自相关的出现记录数量并防止过拟合。最后，我们获取了131条水晶兰记录，并将它们保存在“.csv”文件中进行建模，如图1。
注：该图基于自然资源部标准底图服务网站下载的审图号为GS(2024)0650号的标准地图制作，底图无修改。
Figure 1. Distribution points of Monotropa uniflora L.
图1. 中国水晶兰分布点
2.1.2. 环境因子数据
本研究从气候、地形、土壤及人类活动四个维度出发，初步选取了包括19个生物气候变量、3个地形变量、10个土壤变量及1个人类活动数据在内的33个环境变量，作为模型构建的初始变量集。其中，气候数据来自世界气候数据库(http://www.worldclim.org)，包括当前和未来(2021~2040年、2041~2060年和2061~2080年) 19个生物气候因子，未来气候数据采用CMIP6的BCC-CSM2-MR气候模式，选取三种共享社会经济路径(SSP)情景：SSP126 (低排放情景)、SSP245 (中排放情景)、SSP585 (高排放情景)，空间分辨率均为30″；地形数据获取自世界气候数据库2.1版本(https://worldclim.org)；土壤数据由10种土壤特性组成，来自全球SoilGrids 2.0版数据库(https://soilgrids.org)；人类影响指数(Human Influence Index, HII)来自NASA社会经济数据与应用中心(https://www.earthdata.nasa.gov)。考虑到地形和土壤变量在短时间内相对稳定，本研究假设其在未来气候下保持不变。
所有环境数据均利用ArcGIS 10.8.1进行空间一致性处理，统一坐标系、空间分辨率与空间范围。为消除多重共线性对模型精度的影响，对环境变量进行筛选。首先，将所有环境变量带入MaxEnt模型预运行，去除贡献度为0的变量；其次，将环境变量导入SPSS中进行Pearson相关性分析，当环境变量的相关系数大于0.8时，选择MaxEnt模型试运行结果中贡献率较高的一个[23]。最后得到5个气候变量，2个地形变量，4个土壤变量，1个人类活动数据共12个环境变量用于后续建模，如表1。
Table 1. Environment variables used by the model
表1. 模型使用的环境变量
变量 |
环境因子 |
贡献率/% |
置换重要性/% |
bio14 |
最干月降水量 |
28.5 |
3.4 |
bio6 |
最冷月最低温 |
18.6 |
60.5 |
bio18 |
最暖季降水量 |
13.8 |
6.4 |
bio4 |
温度季节性 |
13 |
2.5 |
Elev |
海拔 |
8.2 |
10.9 |
pH.H2O |
土壤pH |
7.2 |
0.8 |
Slope |
坡度 |
3.4 |
2.4 |
SOC |
土壤有机碳含量 |
2.3 |
5.5 |
CLAY |
粘粒 |
1.8 |
5.5 |
HII |
人类影响指数 |
1.5 |
2.8 |
SILT |
粉粒 |
1.1 |
2 |
bio15 |
降水季节性 |
0.5 |
0.9 |
2.2. MaxEnt模型的构建
将水晶兰的分布点数据与筛选出的环境因子导入MaxEnt模型，以模拟预测中国水晶兰的潜在适生区。在模型训练中，采用25%的分布点作为测试集，75%的分布点作为训练集，并进行10次重复运行。在每次重复运行中，使用不同的随机数据划分方式，以减少单次划分带来的误差，最后对10次重复运行测试集的AUC值进行均值化处理。模型输出格式设置为Logistic。使用刀切法构建响应曲线，并通过受试者工作特征(Receiver Operating Characteristic, ROC)曲线验证模型精度。AUC值(即ROC曲线下面积)作为衡量模型预测准确度的重要指标，不受阈值影响。一般来说，AUC值在0.5~0.6时，预测结果较差；0.6~0.7时，预测结果较差；0.7~0.8时，预测结果一般；0.8~0.9时，预测结果较好；0.9~1.0时，预测结果非常好。
2.3. 适生区划分与面积计算
将MaxEnt模型输出的结果导入ArcGIS 10.8.1平台，采用重分类(Reclassify)工具对其进行等级划分。将生境划分为4个等级：非适生区(0 < X < 0.2)、低适生区(0.2 < X < 0.4)、中适生区(0.4 < X < 0.6)、高适生区(0.6 < X < 1) 4类生境。随后，通过GIS工具统计各适生等级的栅格数量并换算为面积(km2)，以量化不同等级适生区的空间占比。
2.4. 适生区动态变化与质心迁移计算
本研究基于ArcGIS的SDMToolbox工具，对水晶兰潜在适生区的变化进行了系统评估。首先，依据最大检验灵敏度与特异性逻辑阈值(设定临界阈值为0.2)，将MaxEnt模型生成的生境适宜度指数进行二值化处理，得到各时期的潜在适生区分布图：适宜度指数高于阈值的区域划为适宜区，低于阈值的区域则为非适宜区。随后，对不同时期的二值分布图进行叠加分析，依据像元变化将适生区动态划分为三类：稳定区、收缩区与扩张区，从而系统量化了研究时段内适生区面积的时空变化特征[24] [25]。此外，我们还计算了各时期水晶兰适生区的质心，并根据质心的经纬度坐标计算其迁移距离。最后，利用合并工具和数据管理工具将质心点数据转换为线要素，形成质心迁移轨迹，直观展示质心的时空迁移路径。该分析结果将与适生区变化分析相结合进行解释，以全面了解适生区动态的幅度和空间格局[26]。
3. 结果与分析
3.1. 模型预测结果
MaxEnt模型在当前气候条件下对水晶兰分布的预测结果显示，AUC值为0.905，如图2。可见MaxEnt模型在预测水晶兰潜在适生区方面有较高的准确性和可靠性。
Figure 2. ROC curve of MaxEnt results
图2. MaxEnt结果ROC曲线
3.2. 影响水晶兰分布的主导因子分析
基于MaxEnt模型运行自动生成的贡献率和刀切法检验结果，对影响水晶兰分布的主导环境因子进行综合评价[27]。由表1可知，贡献率最大的前6个环境变量为最干月降水量(bio14)、最冷月最低温(bio6)、最暖季降水量(bio18)、温度季节性(bio4)、海拔(Elev)、土壤pH (pH.H2O)，贡献率分别为28.5%、18.6%、13.8%、13%、8.2%、7.2%，累计贡献率为89.3%。刀切法检验结果如图3显示，对水晶兰分布影响较大的前6个环境因子依次为冷月最低温(bio6)、温度季节性(bio4)、土壤pH (pH.H2O)、最暖季降水量(bio18)、最干月降水量(bio14)、降水季节性(bio15)，训练增益分别为1.0、0.86、0.79、0.78、0.72、0.51。两者相差不大。综上，可得气候是主导变量，总贡献率为74.4%，其中降水总量占42.8%，对水晶兰分布的影响大于温度(31.6%)。
Figure 3. Jackknife method to test the impact of environmental variables
图3. 刀切法检验环境变量的影响
通过响应曲线分析影响水晶兰分布的环境变量的阈值，设置存在概率阈值为0.5，当其分布概率大于0.5时，则对应的气候因子值为适宜水晶兰生长[28]。如图4可得适合水晶兰生长的环境变量范围为：最干月降水量(bio14)为10.0~55.4 mm、最冷月最低温(bio6)为−5.8℃~5.5℃、最暖季降水量(bio18)为456~1530 mm、温度季节性(bio4) ≤ 703、海拔(Elev)为1167~3663 m、土壤pH (pH.H2O)为5.0~6.1。得出结论：水晶兰偏好分布在中高海拔、全年温度稳定且土壤呈酸性的森林生态系统中，这为其潜在适生区划分和野外保护提供了关键的环境参数依据。
3.3. 当前气候条件下水晶兰潜在分布
当前气候条件下，水晶兰的总适生区面积为201.26 × 104 km2，主要以低适生区为主，面积达103.67 × 104 km2，占总适生区面积的51.51%。中适生区与高适生区面积分别为69.67 × 104 km2与27.92 × 104 km2，占比分别为34.62%与13.87%。从地理空间上看(图5)，水晶兰的适生区分布在90˚E~125˚E、20˚N~35˚N。其中，高适生区集中于云南北部、四川西南部及西藏东南部的山区原始森林中；中适生区主要环绕高适生区，并向东延伸至贵州西部、湖南西部、广西西北部以及福建的中高海拔山区；低适生区分布最为广泛，分布在大部分南方省份。
3.4. 未来不同气候条件下水晶兰潜在适生区预测
在未来不同气候情景(SSP126、SSP245、SSP585)下，水晶兰的潜在适生区面积与空间格局均发生显著变化。在面积变化方面，与当前气候条件相比，水晶兰的适生区面积在未来各时期与情景下整体均呈
Figure 4. Response curves of dominant environmental variables
图4. 主导环境因子响应曲线
注：该图基于自然资源部标准底图服务网站下载的审图号为GS(2024)0650号的标准地图制作，底图无修改。
Figure 5. Current distribution of suitable areas for Monotropa uniflora L.
图5. 水晶兰当前适生区分布情况
扩张趋势(表2)。其中，2061~2080年SSP585情景下的扩张最为显著，总面积增长率达35.68%；2061~2080年SSP126与2041~2060年SSP245情景次之，增长率分别为34.12%与31.39%。在各级适生区中，高适生区面积的扩张最为突出，尤其是在2021~2040年SSP126路径下，高适生区面积激增至110.24 × 104 km2，增长率高达294.84%，显著高于当前水平。在空间分布方面，水晶兰呈现出以秦岭—淮河以南为核心、向北向西扩张的地理格局，整体上，未来气候条件将促使水晶兰的分布范围从传统的南方湿润山区向更高纬度和海拔的边缘地区扩散，形成更广泛的连续分布带。如图6所示。
3.5. 未来不同气候条件下水晶兰潜在适生区分布格局变化
未来不同时期不同气候情景下，水晶兰的适生区整体呈稳定状态。从空间分布来看(图7)，原有核心适生区得到高度保留，扩张区主要出现在当前适生区的北部(如秦岭一带)和西部(向青藏高原边缘或更高海拔地区延伸)，收缩区主要集中在当前适生区的东南沿海和南部边缘(如广东、广西南部)，生境变迁
Table 2. Potential suitable area distribution of Monotropa uniflora L. under different future climate scenarios
表2. 未来不同气候情景下水晶兰潜在适生区分布面积
气候情景 |
面积(×104 km2) |
低适生区 |
中适生区 |
高适生区 |
总适生区 |
增长率 |
当前 |
103.67 |
69.67 |
27.92 |
201.26 |
/ |
2021~2040 |
SSP126 |
62.84 |
59.37 |
110.24 |
232.45 |
15.49% |
SSP245 |
87.93 |
94.35 |
69.85 |
252.12 |
25.27% |
SSP585 |
98.19 |
75.66 |
53.23 |
227.08 |
12.83% |
2041~2060 |
SSP126 |
99.38 |
82.75 |
64.71 |
246.84 |
22.65% |
SSP245 |
100.19 |
95.17 |
69.08 |
264.44 |
31.39% |
SSP585 |
92.47 |
71.22 |
57.51 |
221.19 |
9.90% |
2061~2080 |
SSP126 |
82.28 |
85.55 |
102.10 |
269.94 |
34.12% |
SSP245 |
98.34 |
77.44 |
61.58 |
237.37 |
17.94% |
SSP585 |
103.51 |
79.78 |
89.78 |
273.08 |
35.68% |
注：该图基于自然资源部标准底图服务网站下载的审图号为GS(2024)0650号的标准地图制作，底图无修改。
Figure 6. Potential suitable habitat distribution of Monotropa uniflora L. under different future climate scenarios
图6. 未来不同气候情景下水晶兰潜在适生区分布
注：该图基于自然资源部标准底图服务网站下载的审图号为GS(2024)0650号的标准地图制作，底图无修改。
Figure 7. Changes in the potential suitable habitat distribution of Monotropa uniflora L. under different future climate scenarios
图7. 未来不同气候情景下水晶兰潜在适生区分布格局变化
表现为“南收北扩”。在面积变化方面(表3)，高排放情景在远期表现出最强的扩张潜力，2061~2080年SSP585情景下的扩张最为显著，扩张区面积为96.27 × 104 km2；2061~2080年SSP126与2041~2060年SSP245情景次之，扩张区面积为71.25 × 104 km2和65.66 × 104 km2。
Table 3. Changes in the potential suitable habitat area of Monotropa uniflora L. under different future climate scenarios
表3. 未来不同气候情景下水晶兰潜在适生区分布面积变化
气候情景 |
面积(×104 km2) |
扩张区 |
稳定区 |
收缩区 |
2021~2040 |
SSP126 |
33.22 |
199.02 |
2.28 |
SSP245 |
52.50 |
199.53 |
1.77 |
SSP585 |
33.14 |
193.79 |
7.52 |
2041~2060 |
SSP126 |
48.24 |
198.89 |
2.41 |
SSP245 |
65.66 |
198.16 |
3.14 |
SSP585 |
38.14 |
182.57 |
18.73 |
2061~2080 |
SSP126 |
71.25 |
198.96 |
2.35 |
SSP245 |
45.38 |
191.41 |
9.89 |
SSP585 |
96.27 |
176.81 |
24.49 |
3.6. 质心迁移
注：该图基于自然资源部标准底图服务网站下载的审图号为GS(2024)0650号的标准地图制作，底图无修改。
Figure 8. Migration of the centroid of Monotropa uniflora L. under different climate scenarios
图8. 不同气候情景下水晶兰质心迁移情况
如图8所示，未来情景下，SSP126路径质心从当前(Current)向东北方向迁移75.69 km至2030 s，随后在2050 s向西南方向折返46.92 km，至2070 s又向西北方向迁移65.58 km。SSP245路径下2030 s质心向东微移29.34 km，但到2050 s，质心大幅向东迁移了76.32 km，达到该情景的最东端；随后在2070 s，质心又向西北方向回迁了71.55 km。SSP585路径下持续向西北方向移动。从当前到2030 s迁移了73.58 km，2030 s至2050 s继续向北推进60.79 km，2050 s至2070 s迁移距离激增至115.90 km。
4. 讨论
4.1. 影响水晶兰适生区分布的环境因素
本研究表明，水热条件是限制水晶兰在中国潜在地理分布的关键因素。对于大多数生物而言，与温度和降水有关的环境因子往往是决定其生长发育以及物种地理分布、多样性和进化最重要的参数[29]。模型预测显示水晶兰最适生长范围：最干月降水量(bio14)为10.0~55.4 mm、最冷月最低温(bio6)为−5.8℃~5.5℃、最暖季降水量(bio18)为456~1530 mm、温度季节性(bio4) ≤ 703、海拔(Elev)为1167~3663 m、土壤pH (pH.H2O)为5~6.1。这些环境参数与水晶兰喜阴湿、偏好酸性土壤的生境特征高度吻合。这种环境选择性可能与其特殊的营养获取方式密切相关。由于缺乏叶绿素，水晶兰必须依赖与红菇科(Russulaceae)等真菌形成的共生体系获取碳源，而真菌群落的组成显著受降雨量、气温和湿度等因素影响[30]。另外土壤pH值(5.0~6.1)也进一步证实了这点，这一范围正是大多数森林土壤真菌的最适生长区间[31]。
4.2. 未来气候变化水晶兰适生区的影响
在未来气候变化情景下，水晶兰的潜在适生区面积整体呈扩张趋势，质心迁移轨迹总体向高纬度和高海拔方向移动。造成这种情况的原因可能是由于气候变暖导致植物原适生区的气温升高，原适生区温度逐渐不适宜植物生长[32]。而随着温度的升高和降水格局的改变，原本因寒冷而不适宜的区域也会逐渐转变为潜在的适生生境。因此，植物开始向高海拔地区迁移，以此寻找适宜自身生长的条件。这符合全球气候变化背景下物种分布迁移的普遍规律[33]-[36]。此外，我们发现高适生区大幅增长，适生区延伸至秦岭以北甚至华北部分地区。水晶兰作为菌异养植物，依赖深厚的腐殖质层和特定的红菇科共生真菌。北方地区降水季节性强、土壤有机质积累与南方湿润山地存在明显差异，共生真菌的迁移速度往往滞后于气候变化。因此，模型预测的“大幅增长”更多反映了气候限制的解除，而非物种真实生境适宜。
4.3. 影响水晶兰分布的生态因素及后期验证
尽管我们通过MaxEnt模型预测出了水晶兰的潜在适生区，但实际上水晶兰的分布不仅依赖气候因素，还受到共生真菌和宿主树种分布的限制。与水晶兰共生的真菌主要是红菇属(Russula)和乳菇属(Lactarius)的成员，这些真菌又与壳斗科(Fagaceae)和松科(Pinaceae)等树木形成菌根共生关系[37]。而本研究预测的主要适生区如云贵高原、四川西部等地区，正是壳斗科和松科植物的分布中心[38]。这一空间上的高度重叠，从侧面印证了我们模型预测的合理性。
未来的研究需要从以下方面进行验证和深入探讨：开展针对性的野外调查，在研究预测的适生区内设置样方，记录水晶兰的出现与否及其种群状况；在采集的土壤样本中检测红菇科等关键共生真菌类群的存在与丰度，整合真菌群落数据，与水晶兰的实际分布数据相结合；将气候、宿主植物分布和真菌群落数据结合，构建更精准的分布模型。
4.4. 水晶兰的保护对策
基于当前与未来适生区的时空演变特征，可以对水晶兰进行分级保护。核心区就地保护：在这些区域依托现有体系，维持林下枯落物层厚度，避免清理林底腐殖质，加强对生态系统的整体保护；收缩区抢救性保护，由于水晶兰难以人工栽培，建议开展基于共生真菌分离培养的迁地保护技术研究，或建立模拟原生生境的保护区；扩张区长期监测：对于未来潜在扩张区，应建立长期监测样地，关注水晶兰自然种群的扩散动态。
5. 结论
本文综合考虑了气候、土壤、地形、人类活动等环境因子，评估影响水晶兰地理分布的主导环境因子，预测当前以及未来不同气候情景下水晶兰的潜在适生区。MaxEnt模型AUC值达0.905，表明模型的预测结果具有较高可靠性。通过贡献率和刀切法检验结果得出，气候是影响水晶兰分布的主导环境因子，总贡献率为74.4%。当前气候条件下水晶兰总适生区总面积为201.26 × 104 km2，主要分布在秦岭–淮河以南。未来气候变化情景下水晶兰潜在适生区总面积呈现扩张趋势，高适生区增长显著，分布格局呈现出“南收北扩”，质心迁移轨迹总体往高纬度和高海拔方向移动。虽然预测模型显示适生区显著扩张，但受限于菌异养机制中真菌与宿主的迁移滞后，实际扩张可能受阻。未来，我们将结合更多生态因子和物种间相互作用(共生真菌分布、宿主树种限制)进行分析，进一步提高模型的预测精度，为水晶兰的保护和栖息地恢复提供理论支持。