當(dāng)前(1970—2000)與未來氣候變化情境下中國醉馬芨芨草潛在分布預(yù)測

黃睿杰， 張春艷， 溫雨婷， 吳晨晨， 路 浩， 趙寶玉

(西北農(nóng)林科技大學(xué)動物醫(yī)學(xué)院， 陜西 楊凌 712100)

全球變暖對溫度和降水格局的改變影響著植物的分布格局，因?yàn)樗鼈兊姆植己艽蟪潭壬先Q于非生物因素，如降水、溫度、土壤和海拔等[1]。因此，氣候變化對植物分布的影響是近年來研究所關(guān)注的焦點(diǎn)。對氣候變化如何影響植物分布的研究有助于農(nóng)業(yè)生產(chǎn)和生物多樣性保護(hù)，并促進(jìn)生態(tài)系統(tǒng)的可持續(xù)性。

物種分布模型(Species distribution models，SDMs)可以將物種發(fā)生點(diǎn)數(shù)據(jù)與相關(guān)環(huán)境變量相結(jié)合，預(yù)測物種適宜生存的環(huán)境[2]。最大熵(Max entropy model，MaxEnt)模型是利用最大熵原理預(yù)測物種生境的SDMs，其假設(shè)物種發(fā)生點(diǎn)由物種分布概率確定，而此概率是由環(huán)境因子對發(fā)生點(diǎn)的約束性計(jì)算出來的[3]。過去的研究已證實(shí)MaxEnt相比其他SDMs具有更強(qiáng)的預(yù)測能力和準(zhǔn)確性[3-4]。由于其具有方便易用、對樣本數(shù)量需求小和預(yù)測表現(xiàn)優(yōu)越等優(yōu)點(diǎn)，MaxEnt已廣泛應(yīng)用于生物多樣性保護(hù)、外來入侵物種防治、病蟲害防治和作物栽培等方面。但是目前在毒害草防控領(lǐng)域的運(yùn)用還較少。

醉馬芨芨草(Achnatheruminebrians)，俗稱醉馬草、醉針草、醉針茅等(圖1)，屬于禾本科芨芨草屬多年叢生草本植物，是分布于我國西北、華北天然草地的主要毒害草之一[5]。由于長期的過度放牧和干旱氣候影響，醉馬芨芨草在退化草地競爭中占據(jù)優(yōu)勢，種群數(shù)量和生長面積不斷增加，降低草原生產(chǎn)力，放牧動物誤食或因饑餓被迫采食醉馬芨芨草易引起急性中毒，從而嚴(yán)重影響著草地畜牧業(yè)發(fā)展和草地生態(tài)安全[6]。然而醉馬芨芨草在藥用、造紙、飼草化利用、機(jī)場驅(qū)鳥、生態(tài)等方面又存在其利用價(jià)值[6-7]。目前關(guān)于醉馬芨芨草的研究多為對其化學(xué)成分、生物特性、內(nèi)生真菌、毒理學(xué)等方面研究[8]，而關(guān)于其不同氣候背景下適生區(qū)分布格局的研究尚未報(bào)道。本研究將基于優(yōu)化的MaxEnt模型對醉馬芨芨草當(dāng)前和未來(2050s和2070s)4種氣候情景下適生區(qū)分布進(jìn)行模擬和預(yù)測，并探究影響其分布的關(guān)鍵環(huán)境變量。研究將為醉馬芨芨草的防治和資源利用、草地生態(tài)學(xué)研究和草地放牧管理提供理論依據(jù)和參考。

圖1 醉馬芨芨草單株及其生境Fig.1 Individual plant and habitat of A. inebrians

1 材料與方法

1.1 醉馬芨芨草發(fā)生點(diǎn)數(shù)據(jù)

研究所用醉馬芨芨草發(fā)生點(diǎn)數(shù)據(jù)共515個，來源于：(1)本課題組2014—2021年在中國西北地區(qū)的野外調(diào)查，使用手持式GPS(GARMIN GPSMAP 621 sc)記錄植物群落野外生長的經(jīng)緯度；(2)中國數(shù)字植物標(biāo)本館(https://www.cvh.ac.cn/)和全球生物多樣性信息機(jī)構(gòu)(Global Biodiversity Information Facility，https://www.gbif.org/)的網(wǎng)絡(luò)數(shù)據(jù)庫。對于一些沒有地理坐標(biāo)但有準(zhǔn)確的地理描述的標(biāo)本，我們使用百度地圖坐標(biāo)拾取系統(tǒng)(https://api.map.baidu.com/lbsapi/getpoint/index.html)查詢坐標(biāo)。

MaxEnt模型的一個假設(shè)是物種的發(fā)生點(diǎn)在所有區(qū)域都是系統(tǒng)或隨機(jī)抽樣的[3]。但是在實(shí)際調(diào)查中，往往會受到各種原因?qū)е虏蓸悠?，從而使模型過擬合影響預(yù)測能力[9-11]。因此，在刪除錯誤點(diǎn)后，采用ENMTools軟件(https://github.com/danlwarren/ENMTools)對發(fā)生點(diǎn)數(shù)據(jù)進(jìn)行篩選[12-13]，最終每個5 km×5 km的柵格內(nèi)只含1個發(fā)生點(diǎn)數(shù)據(jù)。

1.2 環(huán)境變量數(shù)據(jù)

研究考慮了60個可能影響醉馬芨芨草分布的環(huán)境變量，包括19個生物氣候變量、36個氣候變量、5個土壤變量以及海拔(表1)。生物氣候變量雖然由氣候變量衍生，但是二者并不相同[14]。當(dāng)前(1970—2000)與未來(2050s與2070s)的生物氣候、氣候變量以及海拔數(shù)據(jù)來源于世界氣候網(wǎng)站(Worldclim,https://www.worldclim.org/)[15]，而土壤數(shù)據(jù)來自于世界土壤數(shù)據(jù)庫v1.2(https://daac.ornl.gov/)[16]。政府間氣候變化專門委員會(Intergovernmental Panelon Climate Change，IPCC)綜合考慮人口、經(jīng)濟(jì)、技術(shù)進(jìn)步和資源利用等因素，提出了共享社會經(jīng)濟(jì)路徑(Shared socioeconomic pathways，SSPs)來代表全球可持續(xù)發(fā)展路徑、中等可持續(xù)發(fā)展路徑、地方可持續(xù)發(fā)展路徑、不平衡發(fā)展路徑和常規(guī)發(fā)展路徑[17]。研究選擇4種SSPs(SSP126，SSP245，SSP370和SSP585)進(jìn)行預(yù)測，它們分別代表碳排放由低到高的未來氣候情景。在SSP126情景中，全球CO2排放量大幅減少，到本世紀(jì)末全球氣溫升高約1.8℃；在SSP245情景中，CO2排放量在本世紀(jì)中葉開始下降之前徘徊在當(dāng)前水平，到本世紀(jì)末氣溫將上升2.7℃；在SSP370和SSP585情景中，CO2排放量將大幅度增加，到本世紀(jì)末平均氣溫分別上升3.6℃和4.4℃[17]。研究使用北京氣候中心氣候系統(tǒng)模式(BCC-CSM2-MR)作為全球氣候模式(Global climate model，GCM)，因?yàn)橐延醒芯孔C實(shí)該模式非常適合模擬中國的氣候變化[18-19]。

由于環(huán)境變量之間的多重共線性會導(dǎo)致模型過擬合，影響模型的輸出結(jié)果，所以使用ENMTools對環(huán)境變量進(jìn)行相關(guān)性分析[13，20]。如果一組環(huán)境變量的相關(guān)性系數(shù)絕對值大于0.7，則根據(jù)刀切法(Jackknife)檢驗(yàn)結(jié)果除去貢獻(xiàn)較低的那一個[21-22]。

表1 研究使用的環(huán)境變量Table 1 List of environmental variables used in this study

1.3 模型優(yōu)化及評價(jià)

研究使用MaxEnt 3.4.4軟件預(yù)測醉馬芨芨草的適生區(qū)分布[3]，隨機(jī)選取75%的發(fā)生點(diǎn)數(shù)據(jù)用于訓(xùn)練模型，而剩下25%的數(shù)據(jù)用于檢驗(yàn)。采用以往研究中常用的10倍交叉驗(yàn)證(Ten-fold cross-validation)確保模型的準(zhǔn)確性[11]。

之前的研究表明特征類(Feature class,FC)和正則化乘數(shù)(Regularization multiplier,RM)的設(shè)置對MaxEnt模型預(yù)測的穩(wěn)健性影響較大，且默認(rèn)設(shè)置的參數(shù)通常不是最優(yōu)的[23-24]。因此，研究采用R 3.6.3軟件(https://www.r-project.org/)和Kuenm包(https://github.com/marlonecobos/kuenm)優(yōu)化MaxEnt模型。Kumen包通過計(jì)算候選模型的偏ROC(Partial ROC,pROC)、遺漏率(Omission rates,OR)和AICc值綜合考慮模型的統(tǒng)計(jì)顯著性、預(yù)測能力和復(fù)雜性[26]。最終選擇具有統(tǒng)計(jì)學(xué)顯著性(pROC=0)，小于遺漏率閾值(OR<5%)，且AICc值最小的模型參數(shù)[26]。

AUC即受試者工作曲線(ROC)下的面積，是一種不依賴閾值的模型表現(xiàn)評估，常用于評價(jià)MaxEnt模型[11]。AUC值越接近1，則表示模型越可靠。通常0.6

1.4 適生區(qū)等級劃分與變化比較

研究使用最大訓(xùn)練靈敏度加特異性邏輯(Maximum training sensitivity plus specificity logistic,MTSS)閾值對物種的“存在-不存在”進(jìn)行劃分。研究證實(shí)該方法在模型輸出基于僅存在(presence-only)數(shù)據(jù)時(shí)，受“出現(xiàn)點(diǎn)-背景點(diǎn)”比例以及物種流行度影響較小，可降低流行度低的物種的遺漏誤差和流行度高的物種的錯判誤差[27]，且具有客觀性、等效性、判別力強(qiáng)等特點(diǎn)[11]。根據(jù)模型輸出的邏輯值(logistic value)，即適生指數(shù)，將醉馬芨芨草的生境劃分為：非適生區(qū)(0～MTSS)、低適生區(qū)(MTSS～0.3)、中適生區(qū)(0.3～0.5)、高適生區(qū)(0.5～1)。

使用ArcGIS 10.7軟件與SDMTools 2.5 (http://www.sdmtoolbox.org/)[28]來可視化當(dāng)前和未來氣候條件下醉馬芨芨草適生區(qū)分布的差異。即根據(jù)MTSS閾值制作“存在-不存在”的二進(jìn)制圖層，并通過未來各氣候情景適生區(qū)預(yù)測的二進(jìn)制圖層中減去當(dāng)前適生區(qū)模擬的二進(jìn)制圖層觀察變化。使用ArcGis 10.7進(jìn)行地理數(shù)據(jù)分析，并使用Microsoft Excel和GraphPad Prism 8軟件統(tǒng)計(jì)數(shù)據(jù)與繪制統(tǒng)計(jì)圖。

2 結(jié)果與分析

2.1 發(fā)生點(diǎn)與環(huán)境變量篩選

經(jīng)過篩選，最終使用229個醉馬芨芨草發(fā)生點(diǎn)數(shù)據(jù)以及8個環(huán)境變量(bio2，bio15，bio19，elev，t_ece，t_oc，t_ph，tamx12)用于建模。這些環(huán)境變量之間的相關(guān)性系數(shù)絕對值均小于0.7(圖2)。

圖2 用于建模的環(huán)境變量相關(guān)性Fig.2 Correlation of environmental variables used to modeling

2.2 模型優(yōu)化與評價(jià)

使用Kuenm包一共創(chuàng)建了1 160個候選模型，所有候選模型都具有統(tǒng)計(jì)學(xué)意義(pROC=0)，其中139個候選模型符合5%遺漏率標(biāo)準(zhǔn)(圖3)。最終選擇符合5%遺漏率并且AICc值最小的模型。最優(yōu)模型pROC=0、遺漏率=3.51%、AICc=5 335.42，其參數(shù)為：RM=3.3，F(xiàn)C=TH。然后用該RM和FC參數(shù)運(yùn)行MaxEnt，并進(jìn)行10次交叉驗(yàn)證，平均AUC為0.955，標(biāo)準(zhǔn)偏差為0.011 7。這表明模型預(yù)測結(jié)果具有很好的準(zhǔn)確性和可信度，可用于模擬醉馬芨芨草的適生區(qū)分布預(yù)測。

2.3 影響醉馬芨芨草適生區(qū)分布的關(guān)鍵環(huán)境變量

利用刀切法檢驗(yàn)環(huán)境變量對醉馬芨芨草適生區(qū)分布預(yù)測收益的重要性(圖4)。僅使用12月平均最高氣溫(tmax12)時(shí)，所獲得的增益值最高，其次是海拔(elev)、最冷季降水量(bio19)和降水變異系數(shù)(bio15)。這表明這些環(huán)境變量中含有很多有效信息。除去tmax12后，增益值減少最多，其次是bio19，elev和bio15。這說明這些環(huán)境變量中含有許多其他變量中不可替代的重要信息。同時(shí)，上述4個環(huán)境變量對模型的累積貢獻(xiàn)率為96.8%(表2)。因此，tamx12,elev,bio19和bio15是影響醉馬芨芨草分布預(yù)測的關(guān)鍵環(huán)境變量。

圖3 所有候選模型的遺漏率與AICc值自然對數(shù)Fig.3 Omission rate VS the natural logarithm of AICc of all candidate models

圖5為模型輸出邏輯值與4個關(guān)鍵環(huán)境變量的響應(yīng)曲線。4個關(guān)鍵環(huán)境變量與邏輯值之間都為單峰響應(yīng)關(guān)系。tmax12，elev，bio19和bio15分別在—1.83℃，3 082.69 m，4.47 mm和97.65時(shí)，模型輸出邏輯值達(dá)到峰值。根據(jù)MTSS來劃分環(huán)境變量的適生閾值(表2)。tamx12，elev,bio19和bio15的適生閾值分別為—7.5～7.9℃，1 053.8～4 590.9 m，1.6～24.6 mm和25.5～111.5。其他環(huán)境變量的適生閾值如下：bio2(>10.09℃)，t_ece(>0 ds(m-1)，t_oc(>0%)，t_ph(>5.06)。

圖4 刀切法檢驗(yàn)環(huán)境變量對分布預(yù)測收益的重要性Fig.4 Effects of environmental variables on the gain of distribution prediction using Jackknife test

圖5 模型輸出邏輯值與關(guān)鍵環(huán)境變量的響應(yīng)曲線Fig.5 Response curves between logistics output and key environmental variables

表2 關(guān)鍵環(huán)境變量的貢獻(xiàn)率與適生閾值Table 2 Percent contribution and adaptive threshold of key environmental variables

2.4 當(dāng)前醉馬芨芨草適生區(qū)分布

在當(dāng)前氣候條件下，醉馬芨芨草適生區(qū)分布模擬如圖6所示。結(jié)果表明：醉馬芨芨草總適生區(qū)面積約為2.33×106km2，約占中國大陸面積的24.24%，主要分布于新疆、西藏、四川、青海、甘肅、內(nèi)蒙古、寧夏、陜西和山西等省區(qū)。高、中、低適生區(qū)面積分別為5.48×105km2，9.70×105km2和7.87×105km2，分別占中國大陸面積約5.71%，10.10%和8.42%。

2.5 未來醉馬芨芨草適生區(qū)分布

根據(jù)未來氣候數(shù)據(jù)預(yù)測了4種共享社會經(jīng)濟(jì)途徑情景下2050s和2070s醉馬芨芨草適生區(qū)分布(圖7)，并對不同時(shí)期和氣候情景各等級適生區(qū)占中國大陸面積比例與總適生區(qū)平均海拔進(jìn)行統(tǒng)計(jì)分析(圖8)。圖8E顯示，在SSP245，SSP370和SSP585情景下，適生區(qū)面積逐漸增加，其中SSP370情景增加最明顯；在SSP126情景下，2050s增加而2070s減少，但是2070s總適生區(qū)面積相對于當(dāng)前是增加的。圖8 A,B,C和D提示，在4種未來氣候情景下，高適生區(qū)面積變化不大，中適生區(qū)面積有小幅度減小，低適生區(qū)面積有較大幅度的增加。在4種未來氣候情景下，總適生區(qū)平均海拔皆高于當(dāng)前，至少增加118.75 m(圖8 F)。研究將當(dāng)前與未來不同氣候條件下醉馬芨芨草適生區(qū)空間分布進(jìn)行對比(圖9)。結(jié)果顯示，在未來不同情景下山西地區(qū)適生區(qū)有很明顯的縮小，其次是新疆西南部；在新疆西北部、西藏中部、青海中部、內(nèi)蒙古東北部以及黑龍江和吉林的部分地區(qū)有明顯的擴(kuò)張。總體上在東南和西南方向縮減，向西北和東北方向擴(kuò)張。

圖6 當(dāng)前氣候背景下醉馬芨芨草適生區(qū)分布模擬Fig.6 Simulation of distribution of A. inebrians suitable area under the current climate注：審圖號GS(2022)1873，下同Note:On drawing No. GS(2022)1873,the same as below

圖7 未來(2050s,2070s)4種共享社會經(jīng)濟(jì)途徑氣候背景下醉馬芨芨草適生區(qū)分布預(yù)測Fig.7 Prediction of suitable area of A. inebrians under four shared socioeconomic pathways scenarios in the future (2050s,2070s)

圖8 未來氣候背景下各等級適生區(qū)面積占比(A,B,C和D)；各時(shí)期總適生區(qū)面積占比(E)；各時(shí)期總適生區(qū)平均海拔變化(F)Fig.8 Proportion of suitable areas at different levels under future climate background (A,B,C and D)；proportion of total suitable area in each period (E);average elevation change of total suitable area in different periods (F)

圖9 當(dāng)前與未來不同氣候條件下醉馬芨芨草適生區(qū)變化對比Fig.9 Changes in the area of suitable habitat of A. inebrians under different future climate conditions relative to that under current climate conditions

3 討論

本研究利用優(yōu)化后的MaxEnt模型和土壤、地形、氣候和生物氣候等環(huán)境變量，模擬和預(yù)測了在當(dāng)前和未來四種氣候情景下醉馬芨芨草適生區(qū)的分布。最終模型的遺漏率、AICc和AUC分別為3.51%，5 335.42和0.955，環(huán)境變量的響應(yīng)曲線較為平滑。這表明我們的模型輸出結(jié)果是穩(wěn)健可靠的。

溫度、降水、海拔等環(huán)境因素影響植物的生長發(fā)育，同時(shí)也制約著植物物種的分布[29]。刀切法檢驗(yàn)結(jié)果表明，影響醉馬芨芨草分布的關(guān)鍵環(huán)境變量是12月平均最高氣溫(tmax12)、海拔(elev)、最冷季降水量(bio19)和降水變異系數(shù)(bio15)。這些環(huán)境變量對最終模型的累積貢獻(xiàn)達(dá)到了96.8%。利用MTSS對這些環(huán)境變量的適宜性閾值進(jìn)行了劃分。各關(guān)鍵變量的適宜范圍是：tmax12(—7.5～7.9℃)，elev(1 053.8～4 590.9 m)，bio19(1.6～24.6 mm)和bio15(25.5～111.5)。這與之前報(bào)道醉馬芨芨草主要生長在高山及亞高山草原，海拔1 700～4 200 m的山坡、草地、田邊等，具有抗寒、耐旱等生物特性基本相符[6,30]。此外，陳雅琦等[31]證實(shí)醉馬芨芨草中的內(nèi)生真菌使其具有良好的耐鹽堿性。岳永寰等[32]發(fā)現(xiàn)醉馬芨芨草種子在pH 6～10范圍中有良好的發(fā)芽率(>90%)，提示其對極端環(huán)境有極強(qiáng)耐受性。在本研究中，表層土壤鹽度(t_ece)適生閾值>0 ds·m-1，在鹽度持續(xù)升高的情況下依然保持很高的適生指數(shù)(0.66)；表層土壤PH(t_ph)適生閾值>5.06，pH 6～10范圍時(shí)有較高的適生指數(shù)(0.38～0.65)。即醉馬芨芨草有很強(qiáng)的的鹽堿耐受性，與前面的研究一致。

在目前的氣候條件下，醉馬芨芨草適生區(qū)主要分布在新疆、西藏、四川、青海、甘肅、內(nèi)蒙古、寧夏、陜西和山西等省區(qū)，約占中國大陸面積24.24%。其中高適生區(qū)主要分布在新疆、青海、甘肅、寧夏、內(nèi)蒙古、西藏以及四川西部。模擬結(jié)果與已知醉馬芨芨草地理分布基本一致[6]。在未來4種氣候變化情景下(2050s，2070s)，醉馬芨芨草適宜生境的分布面積有所增加。這可能與其強(qiáng)大的環(huán)境適應(yīng)能力有關(guān)。在SSP245，SSP370，SSP585情境下，醉馬芨芨草適生區(qū)呈增加趨勢，SSP370最為明顯，SSP245和SSP585次之。在SSP126情景下，醉馬芨芨草適生區(qū)先增加后減少，但仍然比當(dāng)前適生區(qū)面積要多。這可能是因?yàn)檫m當(dāng)增加溫度可以促進(jìn)植物的擴(kuò)張，而持續(xù)升溫則會帶來一定負(fù)面影響[33]。預(yù)測未來醉馬芨芨草適生區(qū)平均海拔將增加，并且向西北和東北方向擴(kuò)張。這與之前的研究觀點(diǎn)類似。例如，Lenoir等[34]發(fā)現(xiàn)隨著全球變暖氣候發(fā)生變化，大多數(shù)植物的地理分布將會向高緯度和高海拔地區(qū)轉(zhuǎn)移。李曉辰[35]預(yù)測大白刺(Nitrariaroborowskii)、霸王(Zygophyllumxanthoxylon)、沙拐棗(Calligonummongolicum)等多數(shù)荒漠植物的適生區(qū)在未來有向北增加趨勢，且分布重心也將向高緯度地區(qū)遷移。段義忠等[36]預(yù)測四合木(Tetraenamongolica)未來適生區(qū)呈現(xiàn)破碎化分布，并有向東北、西北高緯度地區(qū)遷移的趨勢。

但值得注意的是，本研究存在一定局限性。因?yàn)樯镆蛩亍⑷藶橐蛩?、植被覆蓋度等同樣會影響植物生境適宜性與傳播[37-39]。例如，家畜可促進(jìn)醉馬芨芨草的傳播，尤其是羊更容易增加其種子的擴(kuò)散強(qiáng)度[40]。此外，風(fēng)媒介和水媒介也是醉馬芨芨草種子擴(kuò)散的途徑[41]。人為的過度放牧使草地退化，使醉馬芨芨草在競爭中占有較強(qiáng)優(yōu)勢，從而加速其擴(kuò)增[42]。因此，在實(shí)際應(yīng)用中應(yīng)綜合考慮模擬結(jié)果的不確定性，才能更有效地預(yù)測未來氣候變化條件下醉馬芨芨草適生區(qū)分布格局，從而對這種毒害草的防控與資源化利用做出合理的決策。

4 結(jié)論

綜上所述，利用優(yōu)化后的MaxEnt模型能較好的預(yù)測醉馬芨芨草分布。12月平均最高氣溫(tmax12)、海拔(elev)、最冷季降水量(bio19)和降水量變異系數(shù)(bio15)是主導(dǎo)其分布的關(guān)鍵環(huán)境變量。隨著未來氣候變化，醉馬芨芨草適生區(qū)面積將增加，并向高緯度和高海拔地區(qū)遷移。本研究為醉馬芨芨草的防治和資源利用、草地生態(tài)學(xué)研究和草地放牧管理提供理論依據(jù)和參考。