生物结皮的功能特性对中国北方干旱地区未来气候-土地利用情景的不同反应 魏英欣（Yingxin Wei）、 Ju Mengchen（Mengchen Ju）、 Zou Yanuo（Yanuo Zou）、 Fan Jin（Jin Fan）、 Li Xinhao（Xinhao Li）、 Pang Jingwen（Jingwen Pang）、 Zhang Wenxin（Wenxin Zhang）以及 Bu Chongfeng（Chongfeng Bu）

《Land》：Biocrust Functional Traits Exhibit Divergent Responses to Future Climate–Land Use Scenarios in an Arid Region of Northern China Yingxin Wei, Mengchen Ju, Yanuo Zou, Jin Fan, Xinhao Li, Jingwen Pang, Wenxin Zhang and Chongfeng Bu

　　摘要 生物结皮是干旱生态系统中的关键组成部分，但同时也受到威胁。预测其在未来情景下的功能特征动态对于保护规划至关重要。在中国毛乌素沙地的三个降水带中，我们利用129个出现地点和84个特征采样点，结合MaxEnt栖息地建模和随机森林回归方法，预测了当前条件下以及12种未来气候-土地利用情景（四种社会经济路径（SSP）×三个时间段（2050年代至2090年代）下，苔藓和蓝细菌结皮的功能特征，包括覆盖度、厚度和总体积。研究表明，土壤氮含量、年降水量和土壤钾含量是影响生物结皮栖息地分布的关键环境因素。目前，苔藓结皮覆盖了研究区域的7.63%（厚度为10.56毫米），蓝细菌结皮覆盖了5.88%（厚度为4.88毫米），总体积为4629×10^4立方米。在不同排放和政策组合下，功能特征表现出不同的响应：覆盖度随情景变化而下降，而厚度保持相对稳定。在SSP126情景下，到2070年代苔藓覆盖度下降了3.32%，蓝细菌覆盖度下降了2.80%，总体积减少了2064.76×10^4立方米；到2090年代，苔藓覆盖度有所恢复（+0.26%）。相比之下，SSP370和SSP585情景下结皮持续减少且没有恢复。值得注意的是，蓝细菌的厚度即使在覆盖度下降的情况下也持续增加（+0.02至+0.23毫米），而苔藓的厚度在±0.13毫米范围内波动。高降水量过渡区（362–434毫米）的脆弱性最高，到2070年代苔藓覆盖度在SSP126情景下的下降幅度是SSP585情景下的3倍，并且体积损失一直持续到2090年代。这些发现为面对加速环境压力的过渡区提供了生物结皮特征的空间性预测和定量基准，有助于优先开展保护工作。

1. 引言
干旱地区存在生物土壤结皮（biocrusts）——这些由蓝细菌、藻类、地衣、苔藓及相关微生物组成的紧密表面群落覆盖了地球陆地表面的大约12% [1,2,3]。这些隐花植物群落作为生态系统工程师，估计贡献了全球陆地净初级生产力的7%和高达46%的生物氮固定 [4]，同时稳定了土壤以防止侵蚀 [5]，并调节地表水文 [6]。然而，生物结皮面临着气候变化和土地利用强度增加的快速威胁，预计在未来几十年内其全球覆盖面积将减少25–40% [2,7]。因此，预测生物结皮功能特征（特别是覆盖度和厚度）对环境变化的时空响应对于主动进行干旱地区保护至关重要。然而，实现这一目标仍受到三个相互关联的知识缺口的限制。首先，功能特征研究主要集中在特定梯度上的样地尺度调查，揭示了特征与环境之间的关系 [8,9,10,11,12]，但由于采样限制和尺度转换不确定性，将这些见解外推到连续的区域尺度预测仍然有限。其次，现有的建模方法主要关注栖息地适宜性（存在-缺失）和覆盖度预测 [2,13,14,15]，而作为控制土壤稳定性和碳储存能力的关键特征——生物结皮厚度——很少被纳入区域尺度的功能特征量化中 [16]。第三，气候预测通常针对单一的未来时间点，未能捕捉到不同时间范围内可能出现的非线性趋势或逆转模式 [17,18]。这些缺口共同阻碍了高保护价值区域的识别以及评估何时何地进行干预的紧迫性——这是设计长期环境变化下的具体恢复目标所必需的知识。
毛乌素沙地同时面临着解决这些挑战的迫切需求和独特机会。该地区位于干旱-半干旱过渡带，是黄河流域的关键生态屏障 [19,20]，其生物结皮动态直接影响区域内的固沙和碳循环。自20世纪90年代以来，大规模的复绿工作促进了灌木下的生物结皮广泛发展，但这些群落将如何应对预期的气候变化仍不清楚。尽管之前的研究已建立了该地区的基础样地尺度特征-环境关系 [10,11]，但将这些样地尺度的特征-环境见解与MaxEnt栖息地建模和随机森林回归结合，进行多情景下的区域尺度多特征预测还没有被探索。
为了弥合这一差距，我们开发了一个基于栖息地限制的特征预测框架，改编自 [2] 全球生物结皮分布建模方法。该框架整合了毛乌素沙地的129个出现地点用于MaxEnt栖息地建模，并使用84个功能特征采样点进行随机森林特征量化。具体来说，我们：（1）通过Jackknife方法和主成分分析确定了控制生物结皮出现的环境驱动因素，并绘制了毛乌素沙地的适宜栖息地分布图；（2）使用随机森林模型量化了适宜栖息地内的多种生物结皮功能特征（包括苔藓和蓝细菌结皮的覆盖度和厚度）的当前空间模式；（3）预测了四种气候情景（SSP126–585）下三个时间段（2050年代、2070年代和2090年代）的功能特征动态，以揭示长期趋势，包括潜在的非单调响应和覆盖度与厚度变化之间的对比模式。我们的发现揭示了生物结皮功能特征的当前分布模式及其在气候变化下的预期变化，为优先保护区域和设计适应快速环境变化的生物结皮群落及相关干旱生态系统的管理策略提供了定量基础。

2. 材料与方法
2.1. 研究区域
毛乌素沙地位于中国西北部的鄂尔多斯高原中部（北纬37°27.5′–39°22.5′，东经107°20′–111°30′），位于黄河大拐弯处（图1a）。作为历史上重要的风成沙源区，自20世纪90年代末以来，通过中国的“退耕还林”计划和系统的复绿措施，该地区经历了大规模的生态恢复，显著减少了风蚀和向黄河输送的沉积物 [19,20]。该区域面积为42,363平方公里，海拔在940–1611米之间，位于东亚季风带西北边缘的干旱-半干旱过渡带 [10]。明显的降水梯度（每年237–434毫米，向西北逐渐减少）为研究气候-生物结皮关系创造了理想条件。
图1. 中国毛乌素沙地生物结皮特征建模的研究区域和采样设计。(a) 毛乌素沙地在中国的沙漠带中的地理位置。(b) 沿降水梯度（每年200–450毫米）的采样点分布（补充图S1）。(c) 代表性的苔藓主导生物结皮样本。(d) 代表性的蓝细菌主导生物结皮样本。采样数据结合了29个环境预测因子（补充表S1），用于建模当前分布并预测12种未来气候情景下的变化。该地区的气候为温带大陆性气候，年平均气温为6.0–8.5°C，60–80%的降水量集中在7月至9月。年蒸发量（1800–2500毫米）远超过降水量，形成了水分受限的环境，其中生物结皮提供了关键的生态系统功能 [21]。植被从西部的沙漠草原向东部以 Artemisia ordosica、Salix psammophila 和 Caragana korshinskii 为主的森林草原过渡。生物结皮——包括苔藓主导的（Syntrichia、Pterygoneurum、Tortula）和蓝细菌主导的（Microcoleus、Chroococcidiopsis）群落——在灌木冠层下广泛发育（图1c,d）[11]。

2.2. 数据收集与预处理
生物结皮出现数据（n = 129个地理参照地点）是通过整合三次系统性的野外调查编制的：2016年6月至7月的初步调查，以及2016年的40个地点和2020年的44个地点的详细测量（图1b）。在去除重复坐标后，保留了129个独特的出现地点用于栖息地适宜性建模。这个样本大小超过了区域尺度上稳健MaxEnt建模推荐的最低阈值（n > 50）[22]。功能特征数据在84个采样点收集（大约每504平方公里一个站点），分层覆盖所有三个降水带，以确保对环境梯度的代表性 [22]。
为了进行功能特征描述，我们汇编了两次系统性调查的数据，共84个采样点，分布在三个降水带上。2016年的调查设置了40个30米×30米的样地，在每个样地内的十二个1平方米的正方形区域内目测估计了生物结皮的总覆盖量（苔藓和蓝细菌的总和），并在实验室中从完整的样本中测量了结皮厚度 [10]。2020年的调查在44个地点进行了测量，分别在两个分层微生境中测量了生物结皮的覆盖度和厚度：(i) 树冠下的微生境，使用点框法在两个随机选定的灌木上记录了四个方向（距灌木中心0.3米）的苔藓和蓝细菌覆盖度；(ii) 树冠间的开阔空间，使用相同的点框协议。厚度是在现场使用电子卡尺（MNT-150 T，上海Minitool集团有限公司，上海）测量的 [11]。样地级别的特征值是通过计算树冠下和空间的测量值的面积加权平均值得到的，权重由每个地点的灌木冠层覆盖率决定 [8,23]。尽管两次调查在样地设计和测量协议上有所不同，但都采用了标准的生物结皮调查方法 [8,23]：总覆盖度以百分比表示，便于直接比较2016年的目测值和2020年的点框计数。2016年的完整样本测量和2020年的现场电子卡尺测量代表了生物结皮研究中广泛使用的互补方法 [8]；虽然我们不能完全排除测量协议之间的系统差异，但任何此类偏差将表现为每次调查年份内的空间均匀偏移，而不会与驱动我们随机森林模型的环境梯度相关联，因此不会影响预测的特征变化的空间模式。两次调查的覆盖数据都统一转换为百分比值，厚度测量遵循既定协议进行了标准化 [23]。2016年和2020年84个采样点的年平均降水量分别为404.1毫米和422.9毫米（年际平均值差异小于5%），均在长期降水梯度237–434毫米范围内；由于建模框架使用了30年的气候平均值（BIO1–BIO17）而不是单一年份的值，因此特征与环境关系的学习是独立于年际气候波动的。

2.3. 建模生物结皮栖息地分布和识别环境驱动因素
我们使用了MaxEnt v3.4.4（美国自然历史博物馆，纽约）进行栖息地适宜性建模。我们在研究区域内随机选择了10,000个背景点来定义可用的环境空间 [30]。为了解决空间采样偏差[31]问题，我们通过使用MASS R包（版本7.3-60）中的kde2d函数进行核密度估计来生成了一个偏差校正层[32]，以确保背景数据和出现数据之间的空间偏差模式相似[33]。模型调优评估了208个候选模型，这些模型结合了29种特征类别组合（线性L、二次Q、Hinge H、乘积P、阈值T），并且正则化乘数从0.5增加到4.0（以0.5为增量），使用ENMeval 2.0进行评估[34]。模型性能通过AICc来评估，并根据最低的delta值确定了最佳特征类别（LH）和正则化乘数（3.5）（见补充表S3）。使用这些优化后的参数，我们在MaxEnt中实现了最终模型，并进行了10折交叉验证，以生成稳健的栖息地适宜性预测和性能估计。适宜的栖息地是通过最大测试敏感度加特异性阈值来划定的，该阈值在独立测试数据上最大化了真实技能统计量（TSS）[35]，从而平衡了遗漏误差和错误识别误差[36]。两个辅助阈值（第10百分位的训练存在；敏感性-特异性相等）使得敏感性分析成为可能（见补充表S4）。在当前条件下，模型在最佳阈值下的测试AUC为0.568 ± 0.075，TSS为0.236 ± 0.085（敏感性 = 0.618，特异性 = 0.618）（见补充表S4）。Boyce预测与期望值比率分析进一步确认了模型的良好校准（CBI = 0.855；见补充图S4b）[37]。这些中等程度的区分度指标与建模一个广泛分布的功能群（普遍性约0.55）是一致的，而不是一个分布狭窄的物种[38,39]。对于未来的情景预测，将训练好的MaxEnt模型应用到每个特定情景的环境表面上，并通过将保留的出现点映射到改变后的环境空间中，使用10折交叉验证来评估模型性能[40,41]。所有12个情景下的测试AUC范围从0.743到0.795，TSS范围从0.443到0.533，CBI范围从0.572到0.858（见补充表S5）。

为了量化各个环境变量的相对贡献，我们使用了AUC值进行了Jackknife分析，当每个变量单独使用时（AUC_with_only）[42]。该指标通过使用单个预测因子训练MaxEnt模型来隔离每个变量的独立预测能力，从而揭示出哪些环境因素包含了对生物结皮栖息地适宜性最具区分性的信息，而不考虑变量之间的相互作用[43]。生成了排名靠前的预测因子的响应曲线以确定阈值效应。为了确定最佳的环境范围，我们应用了95%最大适宜性阈值方法：将适宜性≥95%的最大预测适宜性的环境值范围确定为高适宜性区域（见补充表S6）。Jackknife结果和响应曲线（图2a,c）是从MaxEnt的内置输出中提取的；主成分分析（图2b）是使用stats包（prcomp函数）并通过R v4.4.2中的ggplot2和ggfortify进行的[44]。图2. Mu Us沙地中生物结皮栖息地适宜性的环境驱动因素。(a)雷达图显示了29个环境变量单独使用时（AUC_with_only）对MaxEnt栖息地适宜性的相对贡献。变量按类别分组：气候（红色）、植被（绿色）、地形（紫色）、土壤（蓝色）和土地利用（橙色）。(b)基于降雨梯度的生物结皮适宜栖息地像素（n = 244,544）的环境条件主成分分析（见补充图S1）。PC1和PC2分别解释了24.79%和20.35%的方差。(c)前五个环境预测因子对栖息地适宜性的响应曲线。额外的响应曲线见补充图S5。

在MaxEnt划定的适宜栖息地区域内（第2.3节），我们采用了随机森林（RF）模型[46]来预测四种关键生物结皮功能特征——苔藓覆盖率、蓝细菌覆盖率、苔藓厚度和蓝细菌厚度——在当前和12个未来气候-土地利用情景下的空间分布。对于每个情景下的每种特征（当前和12个未来预测），我们使用随机森林变量选择（VSURF）从第2.2节保留的29个环境变量中识别出简洁的预测因子子集。VSURF执行三个连续步骤：根据重要性得分消除无关变量，移除冗余预测因子，并生成一个优化的预测因子子集。每个情景的独立VSURF选择允许模型捕捉到环境条件变化时变量重要性的潜在变化；虽然考虑了从当前条件派生的固定预测因子集，但没有采用，因为它会阻碍在新环境组合下检测到新兴的预测因子[40]。核心预测因子——包括土壤有机碳（Soc；在52个特征-情景组合中的30个中被选中）、NDVI（20/52）和粘土含量（Btcly；17/52）——在所有情景中都得到了一致保留，证实了选择过程的稳健性（见补充图S6）。在变量选择之后，模型超参数——树的数量（ntree）、每次分割时随机采样的变量数（mtry）、每棵树的最大节点数（maxnodes）和每个终端节点的最小观测值数（nodesize）——通过R包‘rBayesianOptimization’[47]使用贝叶斯优化进行了优化。这种优化采用了10折分层交叉验证来最大化预测性能，同时防止过拟合。每种特征-情景组合的优化超参数和选定预测因子详见补充表S7。数据集被分为训练（70%）和测试（30%）两个子集，使用分层随机抽样来确保环境梯度的代表性。84个特征采样点分布在三个环境不同的降水区域（见补充图S1），平均每个504平方公里大约有一个采样点；这种采样配置对训练和测试部分之间的空间自相关性的影响在第4.3节中讨论。最终RF模型使用优化后的超参数进行训练，并在独立的测试数据集上进行了验证。模型性能通过决定系数（R2）、均方根误差（RMSE）、平均绝对误差（MAE）和皮尔逊相关系数（见补充表S7）来评估。在当前条件下，R2值范围从0.45到0.54，皮尔逊相关系数范围从0.52到0.71，表明跨特征的预测能力为中等到好。

为了评估功能特征对未来环境变化的响应，我们将训练好的RF模型应用于12个气候-土地利用情景（四个SSP × 三个时间段：2050年代、2070年代、2090年代）。对于每个情景，特征预测被限制在由MaxEnt预测确定的相应特定情景适宜栖息地区域内，确保覆盖率和厚度估计反映了气候适宜性变化以及环境对生物结皮发育的直接影响。每个特征的当前空间分布如图3所示，而在未来情景下的预测变化如图4和补充图S11-S20所示，相应的统计摘要见补充表S9和S10。图3. Mu Us沙地中生物结皮功能特征的当前空间分布。面板显示了250米分辨率的预测：(a) 蕨类主导的生物结皮覆盖率（平均=7.63%，面积=3231平方公里），(b) 蓝细菌主导的生物结皮覆盖率（平均=5.88%，面积=2493平方公里），(c) 总生物结皮覆盖率（平均=13.51%，面积=5724平方公里），(d) 苔藓厚度（平均=10.56毫米），以及(e) 蓝细菌壳厚度（平均=4.88毫米），(f) 总生物结皮体积（平均=4629 × 10^4立方米）。所有四个功能特征的空间预测不确定性见补充图S7。图4. 在Mu Us沙地中，对比2070年气候情景下生物结皮特征的预测变化。面板(a–f)显示了SSP126（低排放），面板(g–l)显示了SSP585（非常高排放）：(a) 苔藓覆盖率变化（?1406平方公里，?3.32%），(b) 蓝细菌覆盖率变化（?1187平方公里，?2.80%），(c) 总生物结皮覆盖率变化（?2593平方公里，?6.12%），(d) 苔藓厚度变化（?0.04毫米），(e) 蓝细菌厚度变化（+0.05毫米），(f) 总生物结皮体积变化（?2065 × 10^4立方米）；(g) 苔藓覆盖率变化（?225平方公里，?0.53%），(h) 蓝细菌覆盖率变化（?144平方公里，?0.34%），(i) 总生物结皮覆盖率变化（?369平方公里，?0.87%），(j) 苔藓厚度变化（+0.01毫米），(k) 蓝细菌厚度变化（+0.23毫米），(l) 总生物结皮体积变化（?252 × 10^4立方米）。绿色=减少；红色=增加。所有12个情景的完整预测见补充表S9和S10，相应的空间地图见补充图S11–S20。预测不确定性通过10折交叉验证预测的迭代标准差来量化，为特征分布提供了空间上明确的置信估计（见补充图S7–S10）。所有分析均在R v4.4.2（奥地利维也纳的R统计计算基金会；[48]）中进行，使用了terra [49]、randomForest [46]、VSURF [50]和rBayesianOptimization [47]包。

3. 结果 3.1. 生物结皮栖息地分布及其环境驱动因素 在当前条件下，MaxEnt预测的适宜栖息地覆盖了15,290平方公里（占Mu Us沙地的36.09%），在降水梯度上分布不均：高降雨区域（362–434毫米）包含7478平方公里（占区域面积的43.80%），低降雨区域（237–318毫米）包含4425平方公里（占36.54%），中等降雨区域（318–362毫米）包含3383平方公里（占25.69%）（见图2b；补充表S8）。在未来的气候预测下，适宜栖息地面积表现出情景依赖的轨迹，在所检验的12个情景中变化范围从?16.01%（SSP126，2070年代）到+1.37%（SSP126，2090年代）（见补充结果1.1）。通过Jackknife分析确定了栖息地适宜性的前五个预测因子：总土壤钾含量（Tk；AUC = 0.64）、年降水量（BIO12；AUC = 0.62）、总土壤氮含量（Tn；AUC = 0.61）、最干燥月份的降水量（BIO14；AUC = 0.60）和城市土地覆盖（AUC = 0.59）（见图2a）。环境预测因子与栖息地适宜性之间的非线性关系通过响应曲线揭示（见图2c；补充表S6）。总土壤氮含量显示了最强的正响应，当氮含量超过0.79克/千克时，适宜性达到峰值0.93。总土壤钾含量表现出相反的关系，适宜性在低钾浓度（6.8–9.8克/千克）时达到最大值0.84，并且随着钾含量的增加而单调下降。年降水量显示出正的单调响应，在高降雨区域（388–456毫米）时适宜性达到0.70。最干燥月份的降水量表现出单峰响应，在一个狭窄的最佳范围内（1.1–1.6毫米）达到峰值，然后向两个极端方向减少。城市土地覆盖显示出负面影响，适宜性在接近零覆盖时最高（0.68），并且在超过2%时急剧下降，表明生物结皮对人为干扰高度敏感。通过对适宜栖息地像素（n = 244,544）的环境条件进行主成分分析，识别出两个主要的环境梯度（见图2b）。PC1（解释了24.79%的方差）代表了一个降水-生产力-营养梯度，其特征是年降水量（BIO12；载荷=0.85）、粘土含量（Btcly；0.73）、土壤有机碳（Soc；0.71）、植被覆盖比例（FVC；0.71）、总土壤氮含量（Tn；0.80）和次生管理森林（secma；0.76）有很强的正载荷，而海拔（?0.71）和总土壤钾含量（Tk；?0.71）有负载荷。PC2（解释了20.35%的方差）捕捉到了温度-季节性梯度，其中降水季节性（BIO15；0.84）、阳离子交换能力（Cec；0.76）和总磷（Tp；0.75）有正载荷，而等温性（BIO3；?0.84）、土壤pH值（?0.79）和年平均温度（BIO1；?0.67）有负载荷。适宜栖息地像素被聚类为三个不同的组，对应于Jenks分类的降雨梯度：低降雨梯度组占据了负PC1空间，其特征是更高的海拔和钾含量；高降雨梯度组延伸到正PC1空间，具有更高的降水量和植被生产力；中等降雨梯度组位于中间位置。这种清晰的分离证实了降水是结构化Mu Us沙地生物结皮栖息地分布的主要环境因素。

3.2. 未来情景下生物结皮功能特征的预测变化 在当前条件下，观察到生物结皮功能类型的明显空间划分：苔藓主导的生物结皮占据了3231平方公里（占研究区域的7.63%），平均厚度为10.56毫米，而蓝细菌主导的生物结皮占据了2493平方公里（占5.88%），平均厚度为4.88毫米。这两类功能类型合计占据了5724平方公里，总体积为4629 × 10^4立方米（见图3；补充表S9）。值得注意的是，生物结皮的出现集中在高降雨区域（2824平方公里），大约是中等降雨区域（1183平方公里）的两倍。苔藓（10.54–10.64毫米）和蓝细菌（4.83–4.93毫米）的平均厚度在各个降水梯度上相对一致（见补充表S10）。在预测期间，观察到生物结皮覆盖率的显著情景依赖轨迹。到2050年，SSP126下的覆盖率下降最为明显（苔藓：?1.89%；蓝细菌：?1.66%），大约是SSP245下的30–80倍，而总生物结皮覆盖率在SSP126下减少了1502平方公里，但在SSP245下几乎保持稳定（见补充表S9）。到2070年，覆盖率的变化加剧，SSP126下的减少最为显著（苔藓：?3.32%；蓝细菌：?2.80%），大约是SSP585下的5–8倍。独特的是，在此期间SSP245预测了苔藓（+0.37%）和蓝细菌（+0.27%）覆盖率的轻微增加。到2090年，在SSP126情景下，苔藓覆盖情况出现了逆转趋势，从下降转为轻微恢复（+0.26%），而SSP245情景则持续呈现下降趋势（苔藓：?0.95%）（见补充图表S11–S20）。与覆盖面积的变化相比，所有情景和时间段内的苔藓厚度变化都很小。观察到一个显著的脱钩现象：虽然大多数情景下的苔藓覆盖面积减少，但蓝细菌的厚度在整个情景和时间段内持续增加（+0.02至+0.23毫米）。苔藓厚度的变化范围较小（?0.13至+0.08毫米），没有明显的方向性趋势。这种覆盖面积和厚度变化之间的脱钩表明，蓝细菌可能通过垂直生长策略来补偿栖息地的缩减（见补充表S9）。总生物壳体积表现出依赖于情景的变化趋势，结合了覆盖面积和厚度的变化（见补充表S9）。到2070年，SSP126情景下的生物壳体积大幅下降（?2064.76 × 10^4立方米），而SSP245情景下则有所增加（+243.88 × 10^4立方米）。区域间差异明显：到2070年，SSP245情景在低降雨量（+187.34 × 10^4立方米）和中等降雨量（+86.70 × 10^4立方米）地区实现了体积恢复，而高降雨量地区则继续下降（?97.24 × 10^4立方米）。到2090年，SSP126情景下的低降雨量和中等降雨量地区开始恢复，而SSP245情景下的高降雨量地区则持续下降（?302.17 × 10^4立方米）（见补充表S10）。不同地区的响应在降雨梯度上存在显著差异。高降雨量地区对环境变化的敏感度最高，2070年代SSP126情景下的苔藓覆盖面积减少了约3倍（716平方公里），而SSP585情景下仅减少了2倍（228平方公里）。相比之下，中等降雨量地区在不同情景下相对稳定。这些模式表明，接近其降水生态位上限的群落可能最容易受到气候扰动的影响。

4. 讨论
4.1. 控制生物壳栖息地分布的环境因素
Jackknife分析确定土壤氮、年降水量和土壤钾是决定栖息地适宜性的主要因素，揭示了干旱生态系统中的营养-水分耦合控制机制[21,51]。土壤氮对栖息地适宜性的显著正相关响应（当TN > 0.79克/千克时适宜性达到峰值0.93）反映了在养分限制条件下生物固氮作用带来的竞争优势，这与压力梯度假说一致[21,52]。这一现象在Mu Us沙地尤为明显，因为沙质土壤本身氮含量较低，生物壳介导的氮输入成为关键的营养补充[51]。年降水量的正单调响应（最佳范围388–456毫米）证实了水分可用性是半干旱系统中生物壳建立的主要限制因素[21,53]。最干燥月份的降水量（最佳1.1–1.6毫米）的单峰响应表明，生物壳栖息地的适宜性需要最低水分阈值，低于这个阈值会导致干燥超出生理耐受限度[53]，而过多的干旱季节降水量可能促进维管植物的竞争[6]。土壤钾的负相关关系（最佳值6.8–9.8克/千克）与中国西北部高钾土壤常与盐碱条件相关联的观察结果一致[54]，这些条件通过渗透压限制了生物壳的定殖[55,56]。城市土地覆盖率超过2%时栖息地适宜性的急剧下降表明生物壳对人为干扰非常敏感，这与表面踩踏和土壤压实即使在低强度下也会破坏壳结构的证据一致[7,57]。尽管Jackknife分析量化了各个预测因子的贡献，但主成分分析揭示了这些因素如何共同作用以形成环境梯度上的栖息地分布。PC1表示降水量-生产力-养分梯度，适宜栖息地的像素明显分为三个簇，对应不同的降雨区域。这表明降水量作为主要环境过滤器，同时影响土壤水分可用性、植被生产力和养分循环过程，共同决定了栖息地适宜性[51,58]。低降雨量地区（特征为较高海拔和钾含量）与高降雨量地区（植被覆盖度更高且土壤有机碳含量更多）之间的明显分离表明，该地区的生物壳栖息地沿着压力梯度分布，不同端点的环境限制不同[52]。PC2捕捉到了次级的温度-季节性梯度，表明温度制度和降水变异性对栖息地适宜性有额外限制[59]。这两个轴解释了适宜栖息地环境变异的45%，表明虽然降水量相关因素占主导，但仍有很多栖息地异质性可以归因于未测量的局部条件，如微地形和生物相互作用。

4.2. 生物壳特征对未来情景的非线性和不同响应
对未来气候情景的功能特征响应预测显示了明显的非线性动态。最值得注意的是，SSP126情景下的适宜栖息地在2050年代（?10.03%）和2070年代（?16.01%）下降，然后在2090年代（+1.37%）恢复。这种模式可能反映了气候稳定和群落恢复之间的滞后：即使在低排放路径下，中世纪过渡期间的累积升温也超过了生物壳的耐热阈值，而群落的重新组装需要时间[60]。这种依赖于阈值的响应与证据一致，即生物壳对气候变化的敏感性表现出非线性特征[18]，其中适度的环境变化仍在其生理耐受范围内，而更极端的变化则触发群落重组[61]。情景依赖的响应反映了不同的机制路径。实验研究表明，温度升高2–4°C会导致苔藓覆盖面积显著下降（通过改变碳平衡和干燥压力[7]），这解释了在高排放情景下的脆弱性。此外，降水模式的改变——特别是小事件频率的增加——可以在一个生长季节内通过持续的负碳平衡消除以苔藓为主的生物壳[62]。在SSP370和SSP585情景下，温度升高和降水模式的变化共同导致了持续的栖息地丧失，这与预测一致，即增加的干旱减少了全球干旱地区的土壤微生物多样性和丰度[63]。一个值得注意的发现是覆盖面积和厚度响应之间的脱钩：即使空间覆盖面积减少，蓝细菌壳的厚度在整个情景中仍持续增加（+0.02至+0.23毫米）。这种模式可能反映了蓝细菌通过增强细胞外聚合物（EPS）分泌来构建更密集的垂直结构以改善水分保持和耐热性的适应性策略[64]。具体来说，EPS形成了一种多糖基质，增强了持水能力和土壤团聚体的稳定性[65,66]，使蓝细菌能够在干燥条件下巩固垂直结构。此外，蓝细菌和苔藓表现出根本不同的碳吸收策略——苔藓在罕见的大降雨事件中增加光合作用同化，而蓝细菌在变化不定的降水量条件下保持稳定的碳固定[67]——可能将代谢资源重新定向用于垂直强化而非横向扩展。生物相互作用可能进一步加剧这种脱钩：湿润地区的维管植物冠层扩张会竞争性地排除表面的生物壳[68,69]，限制了水平覆盖面积，而蓝细菌则通过剩余空间的垂直生长进行补偿。尽管我们的研究条件下缺乏对这些机制的直接实验验证，但这一模式在多个情景中的一致性及其与记录的生理响应的一致性支持了这种解释[17,70]。然而，我们的相关框架无法完全区分这些过程，省略明确的生物相互作用项可能会高估非生物变量对特征动态的独立影响。这种脱钩表明，仅基于覆盖面积变化的评估可能无法完全捕捉气候对生物壳生态系统功能的影响。区域间气候敏感性的差异也很明显：高降雨量地区（362–434毫米）的脆弱性最大（2070年代SSP126情景下苔藓覆盖面积减少了716平方公里），而中等降雨量地区则相对稳定。这种模式与接近生态位边界的群落更容易受到环境扰动一致[71]。可能有多种共同作用机制导致了这种增加的脆弱性。首先，362–434毫米区域的生物壳可能处于其降水生态位的上限，进一步的升温会增加蒸散作用并减少有效土壤水分，可能触发干旱生态系统特有的非线性阈值响应[61,72]。其次，较高的降水量支持更大的维管植物生产力，通过冠层遮荫和凋落物积累加剧了对生物壳的竞争性排除[68,69]——这种压力可能随着气候驱动的植被变化而加剧。第三，我们的情景中纳入的LUH2土地利用预测部分捕捉了人为压力和降水梯度之间的协变（例如，在生产更高地区的更高放牧强度），尽管这些预测的空间分辨率限制了干扰强度的亚区域变异性检测。

4.3. 限制、保护意义和未来方向
有几个限制需要考虑。首先，由于预测的栖息地缩减，方法学上的一个问题是跨情景解释模型性能指标的问题。未来情景模型显示出更高的区分度（测试AUC：0.568 → 0.743–0.795；TSS：0.236 → 0.443–0.533；见补充表S5），但这并不表明预测能力的改善本身，而是反映了物种分布建模中普遍存在的现象[38,39]。随着适宜栖息地从研究面积的36%减少到20–30%，真正不适宜的背景比例增加，使得分类任务变得更容易。因此，这种改善主要是由特异性提高（0.618 → 0.768–0.921）驱动的——即模型正确识别不适宜区域的能力提高了——而敏感性保持相当，确认模型检测适宜栖息地的能力在不同情景下没有根本变化。类似的验证问题也适用于随机森林特征模型，其中70/30的随机划分不能完全排除训练和测试站点之间的空间自相关。三个设计特征部分缓解了这一风险：稀疏的采样网络（约42,000平方公里内的84个站点）限制了地理上接近的训练-测试对；在三个降雨区进行分层采样（见补充图S1）确保评估涵盖了主要的环境梯度；10折交叉验证和独立保留测试结果之间的一致性（R2 = 0.45–0.54；见补充表S7）确认估计值不是由单一有利划分驱动的。虽然考虑了空间块交叉验证，但由于样本量（n = 84）较小，这会导致测试集过小而无法实现。观察到的R2值低于通常报告的自相关膨胀模型的范围，这与区域尺度上相关特征建模的内在预测限制一致。其次，我们的快照调查无法捕捉年际动态，84个特征采样点虽然分布在所有主要降雨区，但可能无法完全代表区域内的异质性。虽然通过逐层标准差图量化了空间明确的预测不确定性（见补充图表S7–S10），但当前的样本量限制了按地形或土壤类型分层的稳健残差分析；随着扩展采样网络的可用，应优先进行此类诊断。第三，物种分布模型中的空间-时间替代假设可能无法完全捕捉环境变化对动态群落的响应[40,73]。作为相关模型，MaxEnt和随机森林捕捉的是统计关联而非因果机制，无法在新未来条件下明确模拟非线性相互作用或反馈循环[73]；能够捕捉这些动态的过程基础替代方法需要参数化数据（例如，物种特定的生理阈值、种间竞争系数），这些数据在区域尺度上尚未获得[73]。因此，预测的特征变化应解释为在当前特征-环境关系条件下的环境合理轨迹，而不是确定性预测。此外，气候预测中的固有不确定性也会影响我们的预测[41]。为了部分解决模型不确定性，我们使用了14个GCM的集合，并比较了四种SSP情景下的响应；尽管如此，仍需要对独立时间数据进行验证。尽管存在这些限制，我们的框架将生物壳研究从分布映射扩展到了功能特征的量化，使得能够估计生态系统服务提供能力。生物壳对全球生物地球化学循环有重要贡献，约占陆地净初级生产力的7%，以及高达46%的生物固氮[4]。它们还在土壤稳定、水文调节和碳封存中发挥关键作用[6,56]。对于Mu Us沙地来说，作为黄河流域的重要生态屏障，这里生成的空间预测可以用于通过识别高功能价值和高气候脆弱性的区域来指导保护优先级。我们的结果对保护规划有几个启示。首先，SSP126情景下2070年代被确定为关键脆弱窗口，表明如果在这一时期之前实施保护干预措施，可能会最有效，使生物壳群落在压力达到峰值之前建立韧性。我们的模型确定表面干扰是一个关键威胁——栖息地适宜性在城市覆盖率超过2%时急剧下降（第3.1节）——这表明在2060年代之前建立禁牧区和限制基础设施扩张是优先干预措施[74,75]。其次，气候敏感性的空间异质性——高降雨量地区显示出最大的脆弱性——表明保护工作应优先考虑这些生物壳接近生态极限的过渡区。在这些区域内，以 Artemisia ordosica 和 Caragana korshinskii 为主的灌木丛应当成为被动恢复的目标区域——在这些地方，下层植被的生物结皮发育最为旺盛[10,11]——通过禁止牲畜进入和保持植被覆盖来实现恢复。VSURF 分析表明，NDVI（归一化差异植被指数）和土壤有机碳是预测覆盖特征最可靠的因素（见补充图 S6）[76]。第三，SSP126 模型预测到到 2090 年，即使不采取积极的恢复措施，自然群落也有望实现恢复[77]。在需要积极恢复的退化区域，可以使用本地来源的蓝细菌（例如研究区域中主要形成的生物结皮物种 Microcoleus vaginatus）进行接种，并结合微地形改造以增强水分保持能力，这可以将生物结皮的恢复时间从几十年缩短到几年[74,78]。我们框架生成的空间明确的特征预测和脆弱性地图（图 3 和图 4；补充图 S11–S20）为制定区域特定的恢复计划提供了定量基础，使干预措施的时间和强度与当地的生物结皮状况和预测的气候趋势相匹配。未来的研究应朝四个互补的方向发展：首先，应建立长期监测网络，以验证 SSP126 模型预测的非单调变化趋势，并识别生态系统状态转变的早期预警指标；其次，应整合高光谱遥感技术进行动态特征监测，从而在区域尺度上验证空间预测结果并检测时间变化[79]；第三，应开展控制实验，弄清不同情景下响应差异和覆盖度与厚度脱钩的生理机制，特别是EPS（氨甲酰磷酸）在耐受压力中的作用；第四，应开发空间明确的保护框架，将我们的功能特征预测与土地利用规划相结合，确定 Mu Us 沙漠及类似干旱生态系统中的重点保护区域和最佳干预时机。

5. 结论
利用 Mu Us 沙漠地区的分布数据和特征数据，我们结合 MaxEnt 生境模型和随机森林回归方法，预测了当前和未来情景下的生物结皮功能特征。研究得出三个主要结论：首先，土壤氮含量、年降水量和土壤钾含量是影响生物结皮分布的主要环境因素；其次，功能特征表现出不同的响应：在 SSP126 情景下，覆盖度呈下降趋势（苔藓减少 3.32%，蓝细菌减少 2.80%，但到 2090 年苔藓覆盖度有所回升至+0.26%），而总体积在 2070 年减少了 2064.76 × 10^4 立方米；相比之下，生物结皮的厚度保持稳定，即使覆盖度下降，蓝细菌的厚度仍有所增加（从+0.02毫米增加到+0.23毫米）；第三，高降水量过渡带（362–434毫米）的脆弱性最高，在 SSP126 情景下，到 2070 年苔藓覆盖度下降幅度是 SSP585 情景的三倍，而且体积损失一直持续到 2090 年。这些结果表明应在预计的中期压力高峰期（SSP126 情景下的 2070 年）之前优先在这些过渡带进行保护干预，并建立长期监测网络以验证预测趋势。

补充材料
以下支持信息可下载：
https://www.mdpi.com/article/10.3390/land15030436/s1
图 S1：Mu Us 沙漠地区的降水量梯度带
图 S2 和 S3：MaxEnt 模型中使用的环境变量地图
图 S4：MaxEnt 模型的性能诊断指标（AUC 和 CBI）
图 S5：所有环境变量的响应曲线
图 S6：不同情景下的 VSURF 变量选择频率热图
图 S7–S10：四种功能特征的预测不确定性地图
图 S11–S20：10种未来情景下的特征变化预测（SSP126 和 SSP585，见图 4）
表 S1：栖息地适宜性和特征建模所使用的环境变量
表 S2：数据来源和预处理方法
表 S3：208个候选模型的 MaxEnt 模型调优结果
表 S4：适宜栖息地划分的阈值比较
表 S5：未来情景下的 MaxEnt 模型性能指标
表 S6：环境响应曲线分析和最佳范围
表 S7：随机森林模型的超参数和性能指标
表 S8：不同降水量区域的适宜栖息地面积
表 S9：当前和未来情景下生物结皮功能特征的统计摘要
表 S10：不同降水量区域的区域尺度特征预测

补充结果 1.1：未来情景下适宜栖息地的时间动态