中药麻黄为麻黄科植物草麻黄、中麻黄或木贼麻黄的干燥草质茎^[1]。其中, 木贼麻黄(Ephedra equisetina)是麻黄科麻黄属直立小灌木, 具有重要的药用、经济和生态价值^[2—7], 其野生资源濒危, 被列为国家二级保护植物^[8—9]。主要分布于河北、山西、内蒙古、陕西、甘肃、宁夏及新疆等省区, 生长在干旱山地、荒漠、戈壁等处。

木贼麻黄是千百年来抗击瘟疫的一味关键药, 如今在我国抗击新型冠状病毒肺炎疫情中也发挥了积极的作用^[10]。木贼麻黄全株可入药, 其麻黄碱的含量比其他种类都高, 是提制麻黄碱的重要原料^[11]。木贼麻黄具有抗旱、耐盐碱等优势, 对西北地区的水土保持、生态平衡和荒漠化防治具有一定的社会效益和生态效益^[12—13]。在过去20年中, 由于气候巨大差异导致木贼麻黄栖息地的变化、商业市场对麻黄的需求增加、人类不加保护的乱采乱挖及连年的干旱及病虫害等原因, 使麻黄资源面临枯竭^[14—18]。为了持续获得木贼麻黄资源, 部分地区开始人工种植木贼麻黄, 但因生境限制导致人工种植和野生麻黄的药用成分差异很大^[19—20]。此外, 木贼麻黄种植点和品种选择缺乏科学指导以致资源浪费。因此, 木贼麻黄确定最适栖息地和优先保护区对于促进麻黄资源保护与利用至关重要。

本研究通过最大熵(MaxEnt)和ArcGIS软件预测气候变化下木贼麻黄在中国范围内潜在适生区分布, 探索人类活动对木贼麻黄分布的影响, 并结合土地利用类型科学选择木贼麻黄最适栖息地和优先保护区, 为保护和持续利用木贼麻黄资源提供理论基础。

1 材料与方法 1.1 木贼麻黄分布数据

本研究所用木贼麻黄的分布数据有以下3个来源：(1)通过查阅中国国家标本资源平台(http://www.nsii.org.cn/)、中国数字植物标本馆(http://www.cvh.org.cn)、全球生物多样性信息平台(http://www.gbif.org/)获取木贼麻黄在中国的分布点经纬坐标信息。(2)通过2021年6月—2022年4月期间陆续对新疆、甘肃、宁夏和陕西4省部分木贼麻黄天然种群进行了实地勘查, 并整理了58个分布点(附表 1)。(3)通过查阅并整理国内外相关文献中的报道整理出木贼麻黄在中国的地理分布区数据。木贼麻黄在中国分布点最终整理得到269个分布点, 为避免集群效应所致模型过拟合, 利用ArcGIS 10.4.1将所收集数据以10km设置缓冲区, 去除叠置数据, 最终得到201个有效分布点(图 1)。

1.2 自然环境数据

本文所用的自然环境数据包括气候数据、土壤数据和地形数据, 其来源为：(1)从全球气候数据库WordClim(https://www.worldclim.org)下载海拔数据和气候数据, 分别选取当前时期(1970—2000年)气候数据和未来两种排放情景(RCP2.6和RCP8.5)下2050年、2070年的气候数据, 每个时期19个气候变量。(2)土壤数据来源于“国家青藏高原科学数据中心”(http://data.tpdc.ac.cn), 其中包括12个土壤数据。(3)地形数据是通过ArcGIS 10.4.1空间分析模块从海拔中提取坡度和坡向。共选取了35个独立的自然环境变量(表 1)进行研究。

1.3 数据处理

选择的35个自然环境变量间存在一定的相关性, 应避免过度拟合。将自然环境数据导入ArcGIS 10.4.1中, 坐标统一为WGS—1984。利用数据管理工具, 执行批量重采样命令, 对各个来源的数据统一分辨率。通过转换工具, 将35个自然环境变量全部转为“asc”格式, 添加筛选后木贼麻黄的地理分布数据导入MaxEnt 3.4.1软件进行建模运算, 随机选取75%的数据作为训练集用于模型预测, 其余25%数据用于验证模型精度, 模型重复运行5次, 得到自然环境变量贡献度。利用ArcGIS 10.4.1软件对各变量图层进行多重共线性分析, 检验各变量之间的相关性。结合各自然环境变量的贡献率和相关性, 确定主要环境变量。

1.4 MaxEnt模型构建与评估

本研究使用MaxEnt 3.4.1模型预测木贼麻黄在中国的潜在分布。将筛选后的木贼麻黄分布数据和自然环境变量数据导入MaxEnt软件中, 分布数据中随机选取75%的数据作为训练集, 剩余25%的数据作为测试集做模型预测^[21—23], 设置迭代数为10000, 模型重复10次, 通过刀切法对自然环境变量数据进行重检验, 其他参数默认。

计算受试者工作特征曲线(ROC)下面积(AUC)量化MaxEnt模型预测准确性, AUC的取值范围为[0—1], AUC值越接近1, 模型预测结果越精准, 可信度越好^[24]。一般而言, 当AUC≤0.6时认为模型预测结果失败; 0.6＜AUC≤0.7时, 模型预测结果较差; 0.7＜AUC≤0.8时, 模型的预测结果一般; 0.8＜AUC≤0.9时, 模型预测结果良好; 0.9＜AUC≤1.0时, 模型预测结果优秀^[25—26]。

1.5 木贼麻黄适生等级划分

MaxEnt模型输出木贼麻黄在中国范围内的分布概率(P), 使用ArcGIS 10.4.1软件ArcToolbox的重分类(Reclassify)命令, 根据分布概率P值将木贼麻黄的生境适宜性分布区划分为四个适生区等级(表 2)。

1.6 人类活动干扰下木贼麻黄适宜生境分析

人类活动数据来源于国际地球科学信息网络中心(https://sedac.ciesin.columbia.edu)的人类足迹(HF)数据层^[27], 此数据层是由人类土地利用、基础设施建设、人口密度、公路、铁路等8个全球数据层生成的人类影响指数通过归一化得到的, 能客观且全面的反应人类活动的强弱及空间分布状态。对HF的原始数据进行坐标与格式转换等处理^[28]。将处理后的HF、木贼麻黄分布数据和主要自然环境数据导入MaxEnt模型中运行10次, 得到人类活动干扰下木贼麻黄适宜生境, 与自然环境(无人类活动)影响下的适宜生境作比较。

1.7 木贼麻黄潜在适生区与土地利用遥感信息叠加分析

从地理国情监测云平台(http://www.dsac.cn)下载2020年土地利用遥感数据, 根据一级分类对土地进行划分, 一级分类将土地划分为耕地、林地、草地、水域、建设用地、未利用土地和海洋。将木贼麻黄高适生区与土地利用遥感数据进行叠加分析。

2 结果与分析 2.1 自然环境变量的筛选

本研究选取的自然环境变量较多, 而各个自然环境变量之间存在一定的相关性, 一定程度上会出现模型过度拟合, 预测不准确的现象。利用ArcGIS 10.4.1软件提取35个自然环境变量的相关性(图 2), 再将变量导入MaxEnt 3.4.1软件, 模型运算5次, 得到自然环境变量的贡献率。综合考虑相关性较低和贡献率相对较高的自然环境变量, 最终确定21个自然环境变量参与模型建模(表 3)。

2.2 模型预测结果

本研究的ROC曲线验证结果显示：自然环境影响下训练集AUC值为0.919, 测试集AUC值为0.854;人类活动干扰下训练集AUC值为0.948, 测试集AUC值为0.925。其中, 人类活动干扰下AUC值均大于自然环境影响下的AUC值, 添加HF对木贼麻黄地理分布具有一定的影响。

2.3 潜在地理分布及适生区评价

将MaxEnt模型预测结果添加至ArcGIS 10.4.1中, 选择浮点型转换为栅格数据, 对栅格数据进行重分类, 按照分布概率P值从低到高将木贼麻黄的生境适宜区划分为非适生区、低适生区、中适生区和高适生区四个等级, 分别比较自然环境影响和人类活动干扰下木贼麻黄适生区面积变化。

自然环境影响下(图 3), 在RCP2.6情境中, 2050年相对当前时期适生区总面积约增加20.5×10⁴km², 2070年相对当前时期适生区总面积约增加54.6×10⁴km²; 在RCP8.5情景中, 2050年相对当前时期适生区总面积约增加105.9×10⁴km², 2070年相对当前时期适生区总面积约增加162.4×10⁴km²。人类活动干扰下(图 4), 在RCP2.6情境中, 2050年相对当前时期适生区总面积约增加173.5×10⁴km², 2070年相对当前时期适生区总面积约增加40.8×10⁴km²; 在RCP8.5情景中, 2050年相对当前时期适生区总面积约增加39.5×10⁴km², 2070年相对当前时期适生区总面积约增加2.2×10⁴km²。

其中, 总适生区面积最大是RCP8.5-2070年(自然环境影响下), 面积约530.7×10⁴km²; 总适生区面积最小的是当前时期(人类活动干扰下), 面积约207.3×10⁴km²; 高适生区面积最大的是RCP2.6-2050年(自然环境影响下), 面积约196.9×10⁴km²; 高适生区面积最小的是当前时期(人类活动干扰下), 面积约12.0×10⁴km²。自然环境影响下的总适生区面积均多于人类干扰下的总适生区面积, 自然环境影响下的高适生区面积同样多于人类干扰下的高适生区面积。

2.3.1 当前时期木贼麻黄潜在分布

在自然环境影响下(图 5), 木贼麻黄高适生区主要分布于甘肃中南部(兰州、白银、临夏、武威中部、张掖零星分布)、宁夏西北部(中卫、银川、吴忠北部)、陕西中北部(延安、榆林东南角、渭南中部、西安北部、咸阳南缘)、新疆西北部(阿勒泰中北部、塔城西北缘和东南缘、博乐中部、昌吉南部、石河子、伊犁中北部、喀什中西部)、山西西部边缘有零星分布、河南北部有零星分布, 高适生区总面积约38.5×10⁴km², 占国土面积的4.01%。中适生区主要以高适生区为核心向外发散分布, 主要分布于甘肃东部、陕西边缘、山西、河南、河北西缘、内蒙古中部、新疆(天山山脉附近、阿尔泰山脉附近、西缘)、青海零星分布、四川零星分布、山东零星分布, 中适生区总面积约138.8×10⁴km², 占国土面积的14.46%。低适生区以中适生区为核心向外发散分布, 主要分布于甘肃西缘、内蒙古锡林郭勒盟、湖北中部、山东北缘、河北东部、辽宁北缘、吉林西缘、新疆西北缘, 低适生区总面积约191.0×10⁴km², 占国土面积的19.90%。除了上述地区外, 研究区域的其他地区均为非适生区, 这部分地区主要包括新疆南部、西藏、青海西部、云南、四川大部分地区、贵州、重庆、广西、湖南、江西、广东、海南、福建、浙江、安徽、江苏、内蒙古东北部、黑龙江、吉林部分地区、台湾、香港等地区, 非适生区面积约591.7×10⁴km², 占国土面积的61.63%。

在人类活动干扰下(图 5), 木贼麻黄总适生区面积较自然环境影响下的总适生区面积大幅度减少, 分布稀疏零散。其主要零星分布于甘肃、宁夏、陕西、山西、新疆等地。其中, 高适生区面积约12.04×10⁴km², 占国土面积的1.25%, 较自然环境影响下高适生区面积减少2.76%;中适生区面积约45.3×10⁴km², 占国土面积的4.72%, 较自然环境影响下中适生区面积减少9.74%;低适生区面积约149.9×10⁴km², 占国土面积的15.62%, 较自然环境影响下低适生区面积减少4.28%;总适生区较自然环境影响下适生区面积减少16.78%。

2.3.2 未来时期木贼麻黄的潜在分布

选取RCP2.6、RCP8.5情景下2050年和2070年的气候变量, 导入MaxEnt模型中对木贼麻黄的分布进行模拟预测。在这两种情景下, 未来时期总适生区面积与当前时期相比都增加了。

RCP2.6情景下(图 6), 2050年自然环境影响下木贼麻黄总适生区面积约388.9×10⁴km², 占国土面积的40.51%, 较当前时期总适生区分布面积增加了2.14%。其中高适生区面积大幅度增加, 中适生区、低适生区面积均减少，甘肃中北部、宁夏、陕西中北部、山西、河南西部边缘、河南北部边缘、内蒙古、吉林西部边缘、新疆中部及边缘分布、青海中部等地区由中适生区升为高适生区; 新疆中部偏南部由非适生区升为中适生区; 青海西北部、内蒙古西部边缘由非适生区升为中、高适生区; 四川中西部零星、西藏东部边缘、云南北部边缘、湖北北部由低适生区降为非适生区。

2050年人类活动干扰下木贼麻黄总适生区面积约380.8×10⁴km², 占国土面积的39.67%, 较当前时期人类活动干扰下总适生面积增加了18.08%。但较2050年自然环境影响下高适生区面积约减少145.5×10⁴km², 中适生区约增加19.3×10⁴km², 低适生区总适生区面积约增加11.8×10⁴km², 总适生区约减少8.0×10⁴km²。

2070年自然环境影响下木贼麻黄总适生区面积约422.9×10⁴km², 占国土面积的44.06%, 较当前时期分布面积增加了5.70%。其中高适生区面积大幅增加、中适生区面积略增加、低适生区面积减少。与当前时期相比, 高、中、低适生区均向外微微扩散。

2070年人类活动干扰下木贼麻黄总适生面积约248.1×10⁴km², 占国土面积的25.85%, 较当前时期人类活动干扰下总适生区面积增加了4.25%。但较2070年自然环境影响下高适生区面积减少64.7×10⁴km², 中适生区减少96.1×10⁴km², 低适生区面积减少14.0×10⁴km², 总适生区面积减少174.8×10⁴km²。

RCP8.5情景下(图 7), 2050年自然环境影响下木贼麻黄总适生面积约474.2×10⁴km², 占国土面积的49.40%, 较当前分布面积增加了11.03%。其中高适生区大幅度增加, 中、低适生区面积也略增加。高适生区以陕西中部为核心向东扩散至河北西部, 向西扩散至青海东部, 新疆高适生区微微向四周扩散, 内蒙古略有零星分布。四川东部、重庆、贵州、云南中南部由非适生区升为低适生区。

2050年人类活动干扰下木贼麻黄总适生面积约246.8×10⁴km², 占国土面积的25.71%, 较当前时期分布面积增加了4.12%。其中高适生区、中适生区、低适生区面积均增加。但较2050年自然环境影响下高适生区面积减少81.9×10⁴km², 中适生区减少96.4×10⁴km², 低适生区面积减少49.0×10⁴km², 总适生区面积减少227.4×10⁴km²。

2070年自然环境影响下木贼麻黄总适生面积约530.7×10⁴km², 占国土面积的55.29%, 较当前时期分布面积增加了16.92%。其中高适生区大幅度增加, 中、低适生区面积也略增加。甘肃中北部、宁夏、陕西北部、山西中南部、河南、新疆西北部、河北西部边缘、山东西部边缘由中适生区升为高适生区, 贵州、云南、湖南、江西、江苏、上海、浙江、福建等零星地段由非适生区升为低适生区。

2070年人类活动干扰下木贼麻黄总适生面积约209.5×10⁴km², 占国土面积的21.82%。较当前时期分布面积增加了0.23%。其中高适生区和中适生区面积增加, 低适生区面积减少。但较2070年自然环境影响下高适生区面积减少14.4×10⁴km², 中适生区减少116.3×10⁴km², 低适生区面积减少61.3×10⁴km², 总适生区面积减少321.2×10⁴km²。

2.3.3 高适生区质心分布变化

从预测结果可知, 2050年和2070年我国木贼麻黄的适宜区面积均有所增加, 在人类活动干扰下总适生区面积少于自然环境影响下的总适生区面积, 可见人类活动因子对木贼麻黄的生长分布具有一定的影响, 预测得到的适生区的形状并不规则, 在人类活动干扰下适生区较自然环境影响下适生区破碎化。采用高适生区质心位置变化来表现木贼麻黄在不同情景下地理迁移变化(图 8)。自然环境影响下, RCP2.6-2050年高适生区质心相对当前时期有一定程度的东移, 从肃北蒙古自治县(39.39N, 96.78E)移至内蒙古阿拉善右旗(40.59N, 102.95E); RCP2.6-2070高适生区质心从肃北蒙古自治县移至肃北裕固族自治县(39.14N, 98.44E); RCP8.5-2050年和RCP8.5-2070年高适生区质心相对当前时期均有一定程度的东移, 从肃北蒙古自治县移至甘州区(39.16N, 100.71E)(39.28N, 100.49E)。在人类活动干扰下, RCP2.6-2050年高适生区质心相对当前时期有一定程度的东移, 从肃南裕固族自治县(38.49N, 100.44E)移至民勤县(38.75N, 103.65E); RCP2.6-2070年高适生区质心轻微向西移; RCP8.5-2050年高适生区质心轻微向南移, 由肃南裕固族自治县移至青海祁连县(38.12N, 100.53E); RCP8.5-2070年轻微向东北方向移动, 从肃南裕固族自治县移质山丹县(38.79N, 101.04E)。

2.4 木贼麻黄潜在高适生区与土地利用遥感信息叠加分析

通过木贼麻黄高适生区与2020年土地利用类型叠加分析(图 9), 木贼麻黄的高适生区多数已被开发利用, 还有一部分处在未利用的土地上, 主要分布在新疆、甘肃、内蒙古和宁夏地区。在自然环境影响下, 高适生区处于未利用土地上的平均面积约为23.5×10⁴km², RCP2.6-2050年的高适生区未利用土地面积最多, 约为66.3×10⁴km²; RCP2.6-2070年和RCP8.5-2070年的高适生区未利用土地面积较少, 约为15.8×10⁴km²。在人类活动干扰下, 高适生区处于未利用土地上的平均面积约为1.6×10⁴km², RCP2.6-2050年的高适生区未利用土地面积最多, 约为4.6×10⁴km²; 当前时期的高适生区未利用土地面积较少, 约为0.5×10⁴km²。总体而言, 预测结果中显示自然环境影响下未利用土地均多于人类活动干扰下的未利用土地, 未来时期下木贼麻黄高适生区处于未利用土地上的面积均多于当前时期。

3 讨论

在中国麻黄有上千年的用药历史, 时至今日仍是中医的一味常用中药。本研究使用MaxEnt和ArcGIS模型模拟了中国范围内木贼麻黄的潜在分布, 揭示了人类活动对木贼麻黄未来生长分布的影响, 明确了木贼麻黄最适栖息地和优先保护区。这些对木贼麻黄的发现将提高其药用价值、保护木贼麻黄资源和调节生态平衡等做出重要贡献^[29]。

MaxEnt模型已广泛用于物种潜在分布区研究, 玉竹等中药材^[30—32]、梭梭等荒漠植物^{[22, 33]}也应用该模型进行了适宜性区划研究, 并取得了较好的预测结果。前人研究表明, AUC值作为模型预测衡量标准, 区间为[0—1], 值越大表明该模型模拟效果良好^[34]。本研究中, 自然环境下训练集AUC值为0.919, 人类活动干扰下训练集AUC值上升至0.948, 其AUC值均在0.9—1之间, 表明本研究预测结果可靠度高, 模型模拟效果极好^{[24, 34]}。同时, 添加人类活动因子提高了该模型预测精度, 说明人类活动因子和木贼麻黄分布之间的相关性很大。研究结果显示, 人类活动干扰下木贼麻黄的总适生区面积明显小于自然环境影响下总适生区面积, 总适生区面积严重减少且呈破碎化状态, 表明人类活动因子对木贼麻黄适生区的影响主要体现在适生区面积减少和适生区等级退化, 即高适生区大面积向中、低适生区过渡, 多沦为低、非适生区, 人类活动威胁到木贼麻黄资源可持续发展, 并进一步影响其地理分布^[35]。盐碱、温度、降水等自然因子对麻黄属植物生长发育有关键影响^[36—43]。本研究中发现, 酸碱度(pH)、海拔(Altitude)和年平均温度(Bio1)对模型预测结果贡献度较高, 表明盐碱、海拔和温度是木贼麻黄适生区主要影响因素。此外, 在牛蒡子、黄芪、中麻黄等物种的研究中也发现这些因子对其地理分布分布具有重要影响^[44—47], 这与药用、荒漠植物的生物学特性相关。

当前时期木贼麻黄在我国的适生区分布较为广泛, 主要集中在甘肃、宁夏、陕西、山西、四川、河北、河南、内蒙古、新疆等地, 主要以高适生区为核心向外发散分布, 其中高适生区多分布于甘肃、陕西、山西、宁夏和新疆, 这与《中国植物志》记载和多数学者研究的麻黄分布结果基本一致^{[2, 8, 9, 48—51]}。未来时期, 木贼麻黄适生区呈逐渐扩张趋势, 不同干扰、不同情景下在2050年、2070年其适生区面积均有所增加, 表明在全球气候变暖背景下木贼麻黄生存空间较20世纪末期有良好回转的趋势。多数物种在适应气候变化过程中, 核心分布区会向高纬度迁移^[52]。本研究中, 在RCP2.6情境下, 木贼麻黄高适生区质心总体向高纬度迁移(仅自然环境影响下2070年轻微向低纬度迁移); RCP8.5情境下, 质心总体轻微向低纬度迁移(仅人类活动干扰下2070年向高纬度迁移), 总体质心有“东迁”的趋势, 其中迁移方向、迁移距离在不同情景下存在差异, 可能是由于人类活动和气候变暖程度不一、该物种对海拔和温度条件的要求较宽泛有关^[53]。本研究中参与建模的环境指标包括了气候因子、人类活动、土壤和地形等因素, 但仍存在一定的局限性, 同时, 研究尺度的选择也是物种分布预测的不确定性因素之一^[54], 又因气候变暖会增加极端天气出现的频率以致栖息地的地退化和丧失, 所以木贼麻黄实际适生区分布可能要小于此次预测结果^[55]。

近年来气候环境恶化和人类活动严重威胁到木贼麻黄资源可持续利用, 其栖息地限制了人工种植和野生麻黄的主要药用成分, 特别是麻黄碱和伪麻黄碱的含量存在很大差异^{[19—20, 56—57]}。基于前期调查和本研究结果, 建议在新疆北部靠近塔克拉玛干沙漠附近、甘肃武威和张掖等高适生区, 制定木贼麻黄资源人工种植激励机制, 大力发展种植基地; 在甘肃西部、古尔班通古特沙漠、库木塔格沙漠、内蒙古中部等低适生区, 建立自然保护区, 引导和鼓励当地群众参与恢复保护工作。未来木贼麻黄高适生区多以被开发利用, 部分分布在耕地、草地和未利用土地, 建议在绿洲边缘的荒漠区采取封禁保护措施快速恢复栖息地, 在未利用土地因地制宜发展木贼麻黄种植业。

4 结论

迄今气候变化和人类活动干扰对木贼麻黄适生区的预测研究鲜有报道。本研究使用MaxEnt模型评估和预测木贼麻黄在中国的适生区, 结果显示：(1)模型AUC值均在0.9以上, 表明模拟效果好; (2)未来木贼麻黄适生区呈扩张趋势, 影响其地理分布的主要因子为酸碱度、海拔、年平均温度; (3)未来木贼麻黄高适生区多数已被开发, 剩余部分多分布于耕地、草地和未利用土地上, 几何质心有向“东迁”的趋势; (4)人类活动干扰下适生区面积严重减少且呈破碎化状态, 今后应密切关注人类活动对木贼麻黄地理分布的影响; 建议日后在新疆、甘肃、内蒙古等地区进行自然恢复和人工恢复相结合的治理方式, 推动麻黄资源和天然药物产业的可持续发展。

附表: