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

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

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

全球變暖對溫度和降水格局的改變影響著植物的分布格局，因為它們的分布很大程度上取決于非生物因素，如降水、溫度、土壤和海拔等[1]。因此，氣候變化對植物分布的影響是近年來研究所關(guān)注的焦點。對氣候變化如何影響植物分布的研究有助于農(nóng)業(yè)生產(chǎn)和生物多樣性保護，并促進生態(tài)系統(tǒng)的可持續(xù)性。

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

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

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

1 材料與方法

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

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

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

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

研究考慮了60個可能影響醉馬芨芨草分布的環(huán)境變量，包括19個生物氣候變量、36個氣候變量、5個土壤變量以及海拔(表1)。生物氣候變量雖然由氣候變量衍生，但是二者并不相同[14]。當前(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)濟、技術(shù)進步和資源利用等因素，提出了共享社會經(jīng)濟路徑(Shared socioeconomic pathways，SSPs)來代表全球可持續(xù)發(fā)展路徑、中等可持續(xù)發(fā)展路徑、地方可持續(xù)發(fā)展路徑、不平衡發(fā)展路徑和常規(guī)發(fā)展路徑[17]。研究選擇4種SSPs(SSP126，SSP245，SSP370和SSP585)進行預測，它們分別代表碳排放由低到高的未來氣候情景。在SSP126情景中，全球CO2排放量大幅減少，到本世紀末全球氣溫升高約1.8℃；在SSP245情景中，CO2排放量在本世紀中葉開始下降之前徘徊在當前水平，到本世紀末氣溫將上升2.7℃；在SSP370和SSP585情景中，CO2排放量將大幅度增加，到本世紀末平均氣溫分別上升3.6℃和4.4℃[17]。研究使用北京氣候中心氣候系統(tǒng)模式(BCC-CSM2-MR)作為全球氣候模式(Global climate model，GCM)，因為已有研究證實該模式非常適合模擬中國的氣候變化[18-19]。

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

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

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

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

之前的研究表明特征類(Feature class,FC)和正則化乘數(shù)(Regularization multiplier,RM)的設(shè)置對MaxEnt模型預測的穩(wě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包通過計算候選模型的偏ROC(Partial ROC,pROC)、遺漏率(Omission rates,OR)和AICc值綜合考慮模型的統(tǒng)計顯著性、預測能力和復雜性[26]。最終選擇具有統(tǒng)計學顯著性(pROC=0)，小于遺漏率閾值(OR<5%)，且AICc值最小的模型參數(shù)[26]。

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

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

研究使用最大訓練靈敏度加特異性邏輯(Maximum training sensitivity plus specificity logistic,MTSS)閾值對物種的“存在-不存在”進行劃分。研究證實該方法在模型輸出基于僅存在(presence-only)數(shù)據(jù)時，受“出現(xiàn)點-背景點”比例以及物種流行度影響較小，可降低流行度低的物種的遺漏誤差和流行度高的物種的錯判誤差[27]，且具有客觀性、等效性、判別力強等特點[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]來可視化當前和未來氣候條件下醉馬芨芨草適生區(qū)分布的差異。即根據(jù)MTSS閾值制作“存在-不存在”的二進制圖層，并通過未來各氣候情景適生區(qū)預測的二進制圖層中減去當前適生區(qū)模擬的二進制圖層觀察變化。使用ArcGis 10.7進行地理數(shù)據(jù)分析，并使用Microsoft Excel和GraphPad Prism 8軟件統(tǒng)計數(shù)據(jù)與繪制統(tǒng)計圖。

2 結(jié)果與分析

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

經(jīng)過篩選，最終使用229個醉馬芨芨草發(fā)生點數(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)化與評價

使用Kuenm包一共創(chuàng)建了1 160個候選模型，所有候選模型都具有統(tǒng)計學意義(pROC=0)，其中139個候選模型符合5%遺漏率標準(圖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ù)運行MaxEnt，并進行10次交叉驗證，平均AUC為0.955，標準偏差為0.011 7。這表明模型預測結(jié)果具有很好的準確性和可信度，可用于模擬醉馬芨芨草的適生區(qū)分布預測。

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

利用刀切法檢驗環(huán)境變量對醉馬芨芨草適生區(qū)分布預測收益的重要性(圖4)。僅使用12月平均最高氣溫(tmax12)時，所獲得的增益值最高，其次是海拔(elev)、最冷季降水量(bio19)和降水變異系數(shù)(bio15)。這表明這些環(huán)境變量中含有很多有效信息。除去tmax12后，增益值減少最多，其次是bio19，elev和bio15。這說明這些環(huán)境變量中含有許多其他變量中不可替代的重要信息。同時，上述4個環(huán)境變量對模型的累積貢獻率為96.8%(表2)。因此，tamx12,elev,bio19和bio15是影響醉馬芨芨草分布預測的關(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)境變量的響應曲線。4個關(guān)鍵環(huán)境變量與邏輯值之間都為單峰響應關(guān)系。tmax12，elev，bio19和bio15分別在—1.83℃，3 082.69 m，4.47 mm和97.65時，模型輸出邏輯值達到峰值。根據(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 刀切法檢驗環(huán)境變量對分布預測收益的重要性Fig.4 Effects of environmental variables on the gain of distribution prediction using Jackknife test

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

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

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

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

圖6 當前氣候背景下醉馬芨芨草適生區(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)濟途徑氣候背景下醉馬芨芨草適生區(qū)分布預測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)；各時期總適生區(qū)面積占比(E)；各時期總適生區(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 當前與未來不同氣候條件下醉馬芨芨草適生區(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)境變量，模擬和預測了在當前和未來四種氣候情景下醉馬芨芨草適生區(qū)的分布。最終模型的遺漏率、AICc和AUC分別為3.51%，5 335.42和0.955，環(huán)境變量的響應曲線較為平滑。這表明我們的模型輸出結(jié)果是穩(wěn)健可靠的。

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

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

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

4 結(jié)論

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