基于MaxEnt模擬歐亞大陸氣候變化下叉子圓柏的潛在分布*

王愛君 路東曄 張國盛 黃海廣,2 王 穎 呼斯楞 敖 民

(1.內(nèi)蒙古農(nóng)業(yè)大學(xué)林學(xué)院 呼和浩特 010019； 2.內(nèi)蒙古自治區(qū)林業(yè)科學(xué)研究院 呼和浩特 010010；3.內(nèi)蒙古鄂爾多斯市烏審旗林業(yè)局 鄂爾多斯 017300)

氣候變化對物種分布有至關(guān)重要的影響，小到生物個(gè)體生理代謝、繁育，大到生態(tài)系統(tǒng)結(jié)構(gòu)組成和功能發(fā)揮。近年來城市化進(jìn)程加快，自然生態(tài)環(huán)境面臨加速破碎化，聯(lián)合國政府間氣候變化專門委員會(huì)在第5次評(píng)估報(bào)告(AR5)中指出全球氣候?qū)⒗^續(xù)變暖，至21世紀(jì)末地球的平均溫度將增加0.3～4.5 ℃(相對于1986—2005的基準(zhǔn)氣候條件)，相應(yīng)的降水、輻射、潛在蒸散以及其他氣候變量的季節(jié)性變化也將隨氣候變暖而發(fā)生改變(Stockeretal.， 2014)。生態(tài)環(huán)境變化會(huì)直接影響到物種的地理分布格局及生態(tài)系統(tǒng)結(jié)構(gòu)、功能與穩(wěn)定性(Chenetal.， 2011； Gibsonetal.， 2011； Dielemanetal.， 2015)。應(yīng)用特定的算法分析物種潛在地理分布，已成為區(qū)域生態(tài)學(xué)和生物地理學(xué)研究的熱點(diǎn)之一(Bellardetal.， 2012)。

為了保護(hù)和加強(qiáng)植物資源管理，抑制生物多樣性的減少，人們依據(jù)客觀存在的種群地理位置及其環(huán)境變量信息構(gòu)建物種分布模型模擬不同氣候下植物潛在地理分布變化范圍。所謂生態(tài)位模型(Ecological Niche Modelling，ENM)就是利用物種現(xiàn)有分布點(diǎn)與環(huán)境因子的關(guān)聯(lián)性，推測物種生存所需因子的有效范圍，并將運(yùn)算結(jié)果投射到不同時(shí)間或地理空間中，以概率的形式來預(yù)測物種潛在分布區(qū)域(朱耿平等， 2013)。隨著科學(xué)技術(shù)的發(fā)展，開發(fā)了基于不同算法的多種生態(tài)位模型，如BioClim、DoMain、GARP、GLM、MaxEnt(Elithetal.， 2006； Qiaoetal.， 2012)等，每種方法各有優(yōu)劣。最大熵模型(Maximum Entropy，MaxEnt)是目前使用最為廣泛的生態(tài)位模型(Ahmedetal.， 2015； Barbosaetal.， 2015； Vazetal.， 2015)，在預(yù)測史前地質(zhì)時(shí)期與未來氣候情景方面都有良好的應(yīng)用(Bertrandetal.， 2012)，也是物種保育(Kumaretal.， 2009)、譜系地理(Wangetal.， 2015)、物種潛在分布區(qū)模擬(李璇等， 2018)等研究領(lǐng)域的首選模型。MaxEnt模型的優(yōu)勢在于其特征類的選擇可以運(yùn)行各種數(shù)量的樣本，處理變量間復(fù)雜的相互作用(Wiszetal.， 2008； Van Proosdijetal.， 2016)。在現(xiàn)有的物種分布建模算法中，MaxEnt被證明表現(xiàn)最好，同時(shí)預(yù)測結(jié)果受個(gè)別錯(cuò)誤數(shù)據(jù)的影響最低(Grahametal.， 2008)。管畢財(cái)?shù)?2016)利用中國91個(gè)分布位點(diǎn)信息及氣候數(shù)據(jù)，基于MaxEnt和ArcGIS軟件模擬預(yù)測四照花(Cornuskousasubsp.chinensis)分布格局，推測其冰期避難所。郭飛龍等(2020)基于MaxEnt模型預(yù)測了胡楊(Populuseuphratica)世界范圍內(nèi)潛在適宜區(qū)，認(rèn)為胡楊地理分布受多種環(huán)境變量影響，僅以氣候變量模擬的胡楊適宜分布區(qū)與實(shí)際分布范圍存在較大偏差，而多類環(huán)境變量綜合模擬更能反映其現(xiàn)實(shí)種群的分布特征。生態(tài)位模型不僅在植物方面有廣泛的應(yīng)用，同時(shí)也適用于動(dòng)物。朱井麗等(2018)基于MaxEnt模型對松嫩平原丹頂鶴(Grusjaponensis)秋遷期生境適宜性進(jìn)行了評(píng)價(jià)。最大熵模型與地理信息科學(xué)的有機(jī)結(jié)合，可以預(yù)測不同氣候情景下的物種潛在分布區(qū)，深入理解物種進(jìn)化歷程； 估算分布面積變化情況，探究物種遷移路線，預(yù)估未來種群可能會(huì)受到全球變化的威脅程度，也可推測擁有豐富遺傳多樣性的避難所。

叉子圓柏(Juniperussabina)為柏科(Cupressaceae)刺柏屬(Juniperus)匍匐常綠灌木，因其發(fā)達(dá)的根系網(wǎng)與強(qiáng)大的萌蘗力，使它既可以生長在石質(zhì)山坡、河谷上，也可以在固定、半固定的沙丘及灘地邊扎根，是優(yōu)良的水土保持與固沙造林資源樹種，對維持生態(tài)環(huán)境穩(wěn)定安全具有重要作用。叉子圓柏還可入藥，作飼料和薪炭林，也是城市綠化、盆景園藝的珍貴樹種。叉子圓柏自然地理分布區(qū)集中在歐洲南部、亞洲的中部及東部地區(qū)，是旱中生常綠灌木，屬泛北極植物區(qū)系，具有極強(qiáng)的適應(yīng)性，是重要建群種。由于氣候趨干旱化、人為毀林開荒、不合理放牧等現(xiàn)象頻繁發(fā)生，導(dǎo)致叉子圓柏種群退化嚴(yán)重，大多數(shù)為長勢差、質(zhì)量低的塊狀殘林，因此叉子圓柏灌叢林的保護(hù)與修復(fù)刻不容緩。叉子圓柏在生理生化、栽培技術(shù)、生長與環(huán)境關(guān)系(孟和等， 2010)、生態(tài)學(xué)特性、遺傳多樣性、基因組大小(路東曄等， 2018)、譜系地理(路東曄等， 2020)等方面已開展了廣泛研究。叉子圓柏在歐亞大陸廣大區(qū)域內(nèi)呈斷續(xù)分布且各分布區(qū)生境因子差異大，存在局部地區(qū)資源分布信息不詳，缺乏對叉子圓柏歐亞大陸分布區(qū)生態(tài)環(huán)境變量特征和生態(tài)過程的全面解析(郭飛龍等， 2020)。García-Cervign等(2018)選取西班牙東南部安達(dá)盧西亞自治區(qū)的4個(gè)種群，研究表明冬季降水延遲會(huì)影響叉子圓柏在夏季干旱條件下的生長。本文應(yīng)用MaxEnt、BioClim、DoMain 3種預(yù)測模型和ArcGIS軟件，收集已知的歐亞大陸叉子圓柏種群地理位置信息，基于氣候環(huán)境變量和地形(海拔)數(shù)據(jù)，模擬其潛在地理分布區(qū)，使用AUC值(受試者工作特征(ROC)曲線下方面積)和Kappa值對3種模型進(jìn)行比較分析與篩選。采用MaxEnt模型預(yù)測不同氣候下適宜區(qū)的動(dòng)態(tài)變化，運(yùn)用ArcGIS 10.2軟件“SDM Tools”估算各個(gè)時(shí)期分布區(qū)面積變化，質(zhì)心總體遷移趨勢。旨在探討制約叉子圓柏地理分布的主導(dǎo)環(huán)境變量，明確叉子圓柏對生態(tài)環(huán)境變量的需求，為科學(xué)合理有效開展叉子圓柏種群保護(hù)修復(fù)、種質(zhì)資源保育及人工栽培提供參考，對促進(jìn)歐亞大陸乃至全球叉子圓柏種群和種質(zhì)保育建設(shè)具有現(xiàn)實(shí)意義。

1 材料與方法

1.1 數(shù)據(jù)來源

叉子圓柏地理數(shù)據(jù)： 國家植物標(biāo)本資源庫(http:∥www.cvh.org.cn)、全球生物多樣性信息平臺(tái)(http:∥www.gbif.org)、中國國家標(biāo)本資源平臺(tái)(http:∥www.nsii.org.cn)，查閱收集植物的相關(guān)文獻(xiàn)資料(Adamsetal.， 2017； 2018； 路東曄等， 2020)和野外采點(diǎn)數(shù)據(jù)。經(jīng)檢索，共獲得標(biāo)本采集記錄種群信息和野外采點(diǎn)數(shù)據(jù)2 950條記錄。

環(huán)境變量信息： 末次盛冰期(Last Glacial Maximum，LGM，距今約22 000年)、全新世中期(Middle Holocene，MID，約6 000年前)、當(dāng)前(Current， 1970—2000年)、未來(Future， 2070：2061—2080年)的氣候變量和地形(海拔)數(shù)據(jù)均來源于世界氣候數(shù)據(jù)庫(WorldClim，https:∥www.worldclim.org)。未來氣候情景采用美國大氣研究中心(NCAR)開發(fā)的CCSM4模型(Otto-Bliesneretal.， 2006)的推算數(shù)據(jù)，該模型對氣候變量模擬與預(yù)測的準(zhǔn)確性都比較高(李文慶等， 2019)?；谡g氣候變化專門委員會(huì)發(fā)布的第5次評(píng)估報(bào)告，選取典型濃度路徑(RCPs)，RCP2.6和RCP8.5分別代表未來溫室氣體低和高2種排放情形(Van Vuurenetal.， 2011)，建立預(yù)測模型。即在RCP2.6的情景下，預(yù)測全球平均氣溫上升0.3～1.7 ℃，全球平均海平面上升0.26～0.54 m； 在RCP8.5的情景下，預(yù)測全球平均氣溫上升2.6～4.8 ℃,全球平均海平面上升0.45～0.81 m(Alexanderetal.， 2013)。下載的19個(gè)氣候變量和地形(海拔)數(shù)據(jù)，其空間分辨率均為2.5′。

1.2 數(shù)據(jù)處理

剔除標(biāo)本無詳細(xì)地理位置、信息重復(fù)及人工林的數(shù)據(jù)，利用Google Earth軟件(http:∥earth.google.com)核查刪除部分分布位置不正確的點(diǎn)，確定種群位置經(jīng)緯度。為降低樣點(diǎn)的空間自相關(guān)且與環(huán)境數(shù)據(jù)范圍保持一致，將世界地圖分成5 km×5 km(空間分辨率為2.5′)為1個(gè)樣本單元格網(wǎng)，刪除格網(wǎng)內(nèi)冗余點(diǎn)，共篩選確定267個(gè)有效樣點(diǎn)用于建模。

應(yīng)用ArcGIS10.2軟件(https:∥www.Esri.com)空間分析功能，對19個(gè)氣候變量進(jìn)行相關(guān)性檢驗(yàn)。以變量之間相關(guān)性值> 0.90為標(biāo)準(zhǔn)，保留有生物學(xué)意義且在初始模型中貢獻(xiàn)率大的變量，減少由于各變量間高度相關(guān)而導(dǎo)致模型過度擬合，最終確定Bio1、Bio2、Bio4、Bio7、Bio8、Bio9、Bio12、Bio14、Bio15、Bio17、Bio18、Bio19和Elev(海拔)，共計(jì)13個(gè)環(huán)境變量(表1)。

表1 預(yù)測叉子圓柏地理分布的環(huán)境變量Tab.1 The environmental variables for prediction of geographical distribution of Juniperus sabina

1.3 模型篩選

將267個(gè)叉子圓柏分布數(shù)據(jù)和環(huán)境變量共同導(dǎo)入MaxEnt3.4.1軟件中。特征類選擇(feature class selection)是MaxEnt重要的參數(shù)設(shè)置(邱浩杰等， 2020)，本研究選擇自動(dòng)特征。設(shè)置25%為驗(yàn)證集(testing data)，75%為訓(xùn)練集(training data)，最大迭代次數(shù)為10 000，重復(fù)訓(xùn)練(replicates)為3，同時(shí)交叉驗(yàn)證(crossvalidation)(Phillips， 2010)防止人為統(tǒng)計(jì)化。勾選刀切法(Jackknife)檢驗(yàn)各環(huán)境變量對叉子圓柏分布的貢獻(xiàn)率，Logistic輸出.asc。

同樣將267個(gè)叉子圓柏分布數(shù)據(jù)和環(huán)境變量導(dǎo)入DIVA-GIS7.5軟件(http:∥www.diva-gis.org)(楊芙蓉等， 2019)，設(shè)置25%為驗(yàn)證集(testing data)，75%為訓(xùn)練集(training data)，重復(fù)訓(xùn)練(replicates)為3，將13個(gè)環(huán)境數(shù)據(jù)轉(zhuǎn)換成.grd格式并生成stack數(shù)據(jù)集，在Modeling-BioClim/DoMain模塊中依次進(jìn)行BioClim和DoMain模型預(yù)測。

使用受試者工作特征(receiver operating characteristic，ROC)曲線下方面積(area under curve，AUC值)和Kappa值對模型預(yù)測結(jié)果進(jìn)行檢驗(yàn)。二者取值范圍均為0～1，值越接近1說明模型模擬精確性越高(Swets， 1988)。AUC值大于0.90，預(yù)測效果極好； Kappa值大于0.80，預(yù)測效果極佳(楊芙蓉等， 2019)。AUC值及其標(biāo)準(zhǔn)差共同反映模型的總體水平，AUC值越大，標(biāo)準(zhǔn)差越小，模型的穩(wěn)定性越高(Elithetal.， 2011； 陳新美等， 2012)，且隨機(jī)取樣對模型影響小。

采用ArcGIS10.2軟件，對分布區(qū)結(jié)果進(jìn)行適宜性劃分和可視化?；诓孀訄A柏靈敏度-特異度和最大化(sensitivity-specificity sum maximization approach)閾值0.20，設(shè)定物種出現(xiàn)概率小于0.20的地區(qū)為叉子圓柏非適生區(qū)，出現(xiàn)概率0.20～0.40為叉子圓柏低適生區(qū)，概率0.40～0.60為叉子圓柏中適生區(qū)，概率大于0.60為高適生區(qū)。

將叉子圓柏歐亞大陸分布區(qū)內(nèi)267個(gè)有效樣點(diǎn)地理數(shù)據(jù)與氣候、地形環(huán)境變量進(jìn)行MaxEnt、BioClim、DoMain模型模擬，3種預(yù)測模型的AUC與Kappa值表明，MaxEnt和DoMain模型均可用于叉子圓柏潛在分布區(qū)模擬，且預(yù)測效果好，穩(wěn)定性高(表2)。

表2 3個(gè)模型的AUC和Kappa值Tab.2 AUC and Kappa values of three models

叉子圓柏潛在適宜分布區(qū)橫跨歐亞大陸(圖1)，基本反映出叉子圓柏在歐亞大陸的分布特征，但存在3種模型的原理與算法不同，預(yù)測的適宜分布區(qū)域略有差異。MaxEnt模型預(yù)測亞洲的主要適宜分布區(qū)為亞洲西部(土耳其、格魯吉亞、亞美尼亞)、亞洲中部(哈薩克斯坦、吉爾吉斯斯坦)、亞洲東部(中國)的干旱和半干旱地區(qū)； 歐洲為歐洲南部(西班牙、法國、德國、意大利、捷克、瑞士、奧地利、波蘭、俄羅斯、保加利亞、克羅地亞、塞爾維亞)地區(qū)。DoMain模型預(yù)測的分布結(jié)果整體向北偏移，亞洲的主要適宜分布區(qū)為亞洲中部(哈薩克斯坦)、亞洲東部(中國、蒙古)地區(qū)； 歐洲為歐洲北部(西班牙、法國、德國、捷克、奧地利、波蘭、俄羅斯、烏克蘭、立陶宛、白俄羅斯、羅馬尼亞)地區(qū)。BioClim模型預(yù)測的分布區(qū)進(jìn)一步向北偏移，亞洲的主要適宜分布區(qū)為亞洲西部(土耳其、格魯吉亞、亞美尼亞)、亞洲中部(哈薩克斯坦、吉爾吉斯斯坦)、亞洲東部(中國)地區(qū)； 歐洲為歐洲北部(德國、捷克、奧地利、波蘭、俄羅斯、烏克蘭、白俄羅斯、羅馬尼亞、保加利亞、塞爾維亞、斯洛伐克)地區(qū)。MaxEnt、BioClim、DoMain模型預(yù)測高適生區(qū)面積分別為118.28×103、157.57×103、195.64×103km2，呈遞增規(guī)律，這與3種模型的運(yùn)算原理有關(guān)： MaxEnt模型是基于最大熵理論(馬松梅等， 2010)，在已知的約束條件下，對未知事物不做任何假設(shè)，受取樣點(diǎn)的影響最小，熵最大時(shí)事物最接近真實(shí)狀態(tài)，其預(yù)測結(jié)果比較保守； BioClim模型是基于環(huán)境包絡(luò)原理，假定此環(huán)境氣候因素落在環(huán)境包絡(luò)內(nèi)，則物種能在該環(huán)境中正常生長、繁殖后代，故預(yù)測結(jié)果也相對保守； DoMain模型是根據(jù)Gower距離算法得到點(diǎn)-點(diǎn)的相似矩陣，受取樣點(diǎn)影響較大，其預(yù)測結(jié)果與分布點(diǎn)關(guān)系緊密，分布非連續(xù)且均以樣本點(diǎn)為中心向周圍分散(楊芙蓉等， 2019)?；诹私猸h(huán)境變量對叉子圓柏地理分布貢獻(xiàn)率為研究目標(biāo)，鑒于3種模型的AUC、Kappa值、預(yù)測分布圖比較相近的情況下，MaxEnt模型更具優(yōu)勢。

圖1 MaxEnt、DoMain、BioClim模型模擬叉子圓柏潛在分布Fig. 1 Simulation of potential distribution of J. sabina by MaxEnt, DoMain, BioClim models

1.4 環(huán)境變量重要性評(píng)估

在MaxEnt3.4.1軟件中基于刀切法，計(jì)算不同環(huán)境變量測試增益(training gain)，結(jié)合自動(dòng)生成的環(huán)境變量貢獻(xiàn)率與置換重要值，檢驗(yàn)不同環(huán)境變量限制叉子圓柏地理分布的重要性大小。

1.5 適生區(qū)動(dòng)態(tài)變化及遷移趨勢

基于MaxEnt模型繪制末次盛冰期、全新世中期、當(dāng)前、未來情景下叉子圓柏潛在地理分布區(qū)。運(yùn)用Reclassfy分類工具，統(tǒng)計(jì)各類適宜區(qū)面積。利用ArcGIS 10.2軟件中“SDM Tools”中“Distribution Changes between Binary SDMs”工具，獲得末次盛冰期(LGM)—全新世中期(MID)、全新世中期(MID)—當(dāng)前(Current)、當(dāng)前(Current)—未來(Future)2070年叉子圓柏適生地?cái)U(kuò)張、收縮和穩(wěn)定區(qū)域面積變化值。采用“SDM Tools”中“Centroid Changes(Lines)”工具，將各個(gè)時(shí)期的分布區(qū)域集合為一個(gè)中心點(diǎn)，創(chuàng)建矢量文件標(biāo)注叉子圓柏分布方向及位置變化，探討歐亞大陸叉子圓柏分布區(qū)遷移趨勢。

2 結(jié)果與分析

2.1 制約叉子圓柏分布的主導(dǎo)環(huán)境變量

MaxEnt模型運(yùn)行過程中自動(dòng)生成的環(huán)境變量貢獻(xiàn)率及置換重要值表明，基于氣候和地形2類環(huán)境變量模擬，貢獻(xiàn)率大于5%的分別為年均氣溫Bio1(28.80%±0.86%)、海拔Elev(29.00%±0.36%)、溫度季節(jié)性變化Bio4(13.00%±0.98%)、最干月降水量Bio14(7.10%±2.14%)和最熱季降水量Bio18(5.70%±0.49%)，累積貢獻(xiàn)率為83.6%； 置換重要值排名前5的分別是海拔Elev(45.70%±5.67%)、年均氣溫Bio1(18.70%±2.12%)、溫度季節(jié)性變化Bio4(8.00%±7.08%)、降水量季節(jié)性變化Bio15(7.70%±1.50%)和最干季平均溫度Bio9(6.20%±1.08%)，累計(jì)置換重要值為86.3%(表3)。刀切法檢驗(yàn)不同環(huán)境變量的測試增益結(jié)果表明，年均氣溫Bio1、溫度季節(jié)性變化Bio4、海拔Elev、最干季平均溫度Bio9、年溫度變化范圍Bio7排在前5位(圖2)。統(tǒng)一分析環(huán)境變量貢獻(xiàn)率、置換重要值及測試增益，采用2類環(huán)境變量模擬，地形(海拔)、大氣溫度(年均氣溫Bio1、溫度季節(jié)性變化Bio4)是限制叉子圓柏分布的主導(dǎo)環(huán)境變量，同時(shí)也受大氣降水(降水量季節(jié)性變化Bio15、最干月降水量Bio14)的調(diào)控。

表3 各環(huán)境變量貢獻(xiàn)率Tab.3 Percent contribution of each environmental variable

圖2 刀切法檢驗(yàn)環(huán)境變量的測試增益結(jié)果Fig. 2 Jackknife test to examine training gain results of environmental variables

2.2 叉子圓柏歐亞大陸潛在適生區(qū)的變化特征

利用ArcGIS軟件對MaxEnt模型模擬結(jié)果重新分類，分別統(tǒng)計(jì)叉子圓柏各適生區(qū)面積，即LGM時(shí)期為72.894×103km2、MID時(shí)期為72.245×103km2、當(dāng)前時(shí)期為663.115×103km2、未來RCP2.6和RCP8.5時(shí)期分別為67.116×103km2和 55.355×103km2，呈現(xiàn)先增加后減少的變化趨勢(表4)。MaxEnt 模型模擬結(jié)果表明，亞洲是叉子圓柏的主分布區(qū)，亞洲的適生區(qū)面積在LGM時(shí)期占86.9%、MID時(shí)期占87.0%、當(dāng)前時(shí)期占57.8%、未來RCP2.6和RCP8.5時(shí)期分別占84.1%和79.2%。

表4 各時(shí)期叉子圓柏適宜分布區(qū)面積①Tab.4 Suitable distribution area of J. sabina in each period km2

基于MaxEnt模型模擬預(yù)測可知，叉子圓柏適生區(qū)范圍及面積隨著氣候環(huán)境變化而變化，并呈現(xiàn)先擴(kuò)張后退縮的特征。

從LGM—MID時(shí)期叉子圓柏適生區(qū)面積雖然下降了0.64×103km2，但在歐洲的意大利、瑞士、奧地利境內(nèi)的阿爾卑斯2 000～3 000 m的山脈、俄羅斯西南部與格魯吉亞交界的大高加索山脈，亞洲西部的格魯吉亞與亞美尼亞交界的小高加索山脈、土耳其東北部的大阿勒山，亞洲中部的吉爾吉斯斯坦境內(nèi)的帕米爾-阿賴山脈、哈薩克斯坦境內(nèi)的塔爾巴哈臺(tái)山，亞洲東部的中國境內(nèi)的阿爾泰山、天山、祁連山脈、青海日月山和陰山的亞高山地帶、青海湖、渾善達(dá)克沙地及蒙古的杭愛山和草原植被帶，發(fā)生擴(kuò)張，分布區(qū)面積擴(kuò)增27.32×103km2。在歐洲的法國東南部、意大利中部、保加利亞、塞爾維亞、克羅地亞部分地帶、俄羅斯西南部，亞洲西部的格魯吉亞、亞美尼亞與土耳其東北部的低海拔山脈(1 000～1 500 m)，亞洲中部的吉爾吉斯斯坦境內(nèi)的帕米爾-阿賴山脈、天山山脈，亞洲東部的中國境內(nèi)的天山、祁連山脈、呂梁山、賀蘭山、六盤山、柴達(dá)木盆地東南區(qū)、青海南山，以及歐亞大陸其他低緯度、低海拔、盆地等生境，為主要退縮區(qū)，退縮面積約為31.08×103km2。相對穩(wěn)定區(qū)域的面積達(dá)39.99×103km2，多集中分布于中亞的天山和阿爾泰山山脈(表5，圖3A)。

MID—當(dāng)前時(shí)期，是歐亞大陸叉子圓柏種群擴(kuò)張時(shí)期，適生區(qū)面積從72.245×103km2快速增加到663.115×103km2，分布區(qū)增加面積為559.71×103km2。此期間叉子圓柏分布遍及歐洲南部、亞洲中部及東部(圖3B)。僅在蒙古北部高山草地區(qū)、吉爾吉斯斯坦境內(nèi)的天山、中國青海湖及附近山脈有少部分的收縮，收縮面積約為17.39×103km2。

在RCP2.6的氣候情景下，到2070年，亞洲中部吉爾吉斯斯坦中部天山區(qū)域及亞洲東部蒙古北部高山區(qū)域以及中國祁連山脈和青海日月山是叉子圓柏潛在適生區(qū)面積增加的主要區(qū)域(圖3C)，新增面積約為22.80×103km2(表5)。除歐洲的意大利、瑞士、奧地利境內(nèi)的阿爾卑斯山脈、俄羅斯西南部與格魯吉亞交界的大高加索山脈、格魯吉亞與亞美尼亞交界的小高加索山脈、土耳其東北部的大阿勒山，亞洲中部吉爾吉斯斯坦境內(nèi)的天山山脈，亞洲東部中國境內(nèi)的阿爾泰山和天山山脈、蒙古草原地帶外，其余均處于收縮狀態(tài)。

在RCP8.5氣候情景下，叉子圓柏歐亞大陸潛在分布范圍保持與在RCP2.6氣候情景下相同的縮減趨勢，縮減582.81×103km2(表5)。新增的范圍主要分布于天山、蒙古高山區(qū)和中國境內(nèi)的阿爾泰山、柴達(dá)木盆地，新增22.75×103km2。與RCP2.6氣候情景下增加的范圍相比較，RCP8.5氣候情景下新增的區(qū)域更近一步向高緯度擴(kuò)張(圖3D)。

表5 不同氣候時(shí)期叉子圓柏潛在分布區(qū)的動(dòng)態(tài)變化Tab.5 Dynamic changes of potential distribution areas of J. sabina in different climatic periods

圖3 不同氣候情景下叉子圓柏分布格局動(dòng)態(tài)變化Fig. 3 Distribution pattern dynamic changes of J. sabina under different climate scenariosA： 末次盛冰期(LGM)— 全新世中期(MID)的叉子圓柏適生地?cái)U(kuò)張、收縮和穩(wěn)定區(qū)域； B： 全新世中期(MID)—當(dāng)前(Current)的叉子圓柏適生地?cái)U(kuò)張、收縮和穩(wěn)定區(qū)域； C： 當(dāng)前(Current)—未來(Future)2070年RCP2.6的叉子圓柏適生地?cái)U(kuò)張、收縮和穩(wěn)定區(qū)域； D： 當(dāng)前(Current)—未來(Future)2070年RCP8.5的叉子圓柏適生地?cái)U(kuò)張、收縮和穩(wěn)定區(qū)域。A： The expansion, contraction and unchange regions of J. sabina from LGM to MID; B： The expansion, contraction and unchange regions of J. sabina from MID to Current; C： The expansion, contraction and unchange regions of J. sabina from Current to Future 2070(RCP2.6); D： The expansion, contraction and unchange regions of J. sabina from Current to Future 2070(RCP8.5).

2.3 叉子圓柏分布質(zhì)心的變化

叉子圓柏在末次盛冰期的分布中心位于中國新疆西部(25° N，100° E)，全新世中期位于中國新疆西北部(30° N，105° E)，當(dāng)前的分布中心位于哈薩克斯坦西南部(28° N，65° E)，在RCP2.6(2070)氣候情景下的分布中心位于中國新疆北部(33° N，110° E)，RCP8.5(2070)氣候情境下的分布中心位于蒙古西部(35° N，113° E)。叉子圓柏的分布中心從末次盛冰期至全新世中期，整體向東北遷移； 從全新世中期到當(dāng)前，出現(xiàn)一個(gè)回折現(xiàn)象； 但從當(dāng)前至未來，又繼續(xù)呈現(xiàn)向北延伸的趨勢。

2.4 叉子圓柏適宜生境因子的閾值范圍

環(huán)境變量貢獻(xiàn)率、置換重要值表明，年均氣溫和海拔的貢獻(xiàn)率相近(29%)(表3)，但是在置換重要值方面海拔明顯大于年均氣溫，對叉子圓柏歐亞大陸潛在地理分布格局起著主導(dǎo)作用。將年均氣溫(Bio1)、海拔(Elev)、溫度季節(jié)性變化(Bio4)、年降水量(Bio12)依次導(dǎo)入MaxEnt模型中進(jìn)行單因子建模，以存在概率大于0.40為閾值計(jì)算4個(gè)主要?dú)夂颦h(huán)境因子的適宜范圍。歐亞大陸叉子圓柏適生區(qū)年均氣溫適宜范圍約為0～7.5 ℃，海拔適宜范圍為900～2 600 m，溫度季節(jié)性變化適宜范圍為5.5～6.5 ℃，年均降水量的適宜范圍為200～800 mm。這些數(shù)據(jù)總體可以定量地描述叉子圓柏適宜生境的氣候環(huán)境特征，與孟和等(2010)的報(bào)道基本一致。

3 討論

3.1 模型模擬結(jié)果的可靠性

樣本采集范圍及樣本量是決定物種分布模型模擬結(jié)果可靠性的重要因素。一般而言，收集到的樣點(diǎn)覆蓋度均勻、區(qū)域廣、樣本量大，獲取到的環(huán)境信息則越豐富，模型建立約束條件越多(莊鴻飛等， 2018)，預(yù)測物種潛在分布圖就越精準(zhǔn)，模型估計(jì)精度越高(趙曉冏等， 2018)。AUC值被公認(rèn)為是模型預(yù)測準(zhǔn)確度的最佳衡量指標(biāo)(Hughesetal.， 2003)。Kappa統(tǒng)計(jì)量是一種閾值依賴評(píng)價(jià)指標(biāo)，也被大量用作模型評(píng)價(jià)參數(shù)(Maneletal.， 2010； Seguradoetal.， 2004)。模擬過程中篩選的環(huán)境變量是否具有代表性、完整性以及種群的分布位置、環(huán)境變量的相關(guān)性，同樣制約著評(píng)價(jià)結(jié)果是否可靠。本研究基于環(huán)境變量相關(guān)性分析，剔除相關(guān)系數(shù)大于0.9的變量，MaxEnt模型模擬結(jié)果較能體現(xiàn)叉子圓柏符合柏科分布帶的特征，接近前人的描述(毛康珊， 2010)。叉子圓柏歐亞分布區(qū)的樣本能代表其分布區(qū)的生境，MaxEnt模型模擬叉子圓柏在中國適宜分布區(qū)的面積(129.643×103km2)大約是王林和等(2014)報(bào)道的基于衛(wèi)星遙感影像面積(0.032 7×103km2)的3 000倍以上，與郭飛龍等(2020)、劉超等(2018)和王茹琳等(2017)的研究結(jié)果類似，都存在模擬面積大于實(shí)際分布面積的情況。這可能是僅使用單一非生物環(huán)境變量，沒有考慮限制植物分布的其他生物因素，如： 繁殖能力、種間競爭、人類干擾等(Dülgeroluetal.， 2019)； 也有可能是由于物種自身較弱的遷移能力和存在遷移屏障造成物種沒有到達(dá)其適宜生境(Sobernetal.， 2005； 2009)。在某些地區(qū)叉子圓柏更新庇護(hù)植物烏柳(Salixcheilophila)、沙棘(Hippophaerhamnoides)等被砍伐后無法更新擴(kuò)展形成叉子圓柏灌叢林(山中典和等， 2000)，此外城市擴(kuò)張、土地用途變更、樵采、地表淺水過度利用等人類社會(huì)活動(dòng)也是造成原始種群消失的重要因素，可見生物因子嚴(yán)重影響著叉子圓柏的實(shí)際分布狀況。

3.2 制約叉子圓柏分布的環(huán)境變量

植物種群的分布特征是生物和非生物因素共同作用的綜合結(jié)果，也是植物對環(huán)境響應(yīng)的最終體現(xiàn)。叉子圓柏為北溫帶分布類型中的古北極植物種，至少起源于第三紀(jì)，曾是北半球森林植被的重要組成成分。MaxEnt模擬預(yù)測顯示叉子圓柏適宜分布區(qū)與實(shí)際分布擬合程度大體一致。環(huán)境變量貢獻(xiàn)率、置換重要值及測試增益顯示，地形(海拔)、氣候(溫度和降水)對叉子圓柏地理分布格局貢獻(xiàn)最大，這與石質(zhì)山地是叉子圓柏的主要生境相關(guān)。溫度和水分的水平分布特征決定著植物生理活動(dòng)和生化過程(Wangetal.， 1998)，而隨著海拔的變化對氣候特征進(jìn)行了重塑，在隱域植被形成中起著關(guān)鍵作用，水熱因子在塑造植物生態(tài)適應(yīng)方面發(fā)揮著重要作用(李文慶等， 2019)，此外土壤性質(zhì)、水文地質(zhì)等立地相關(guān)的環(huán)境變量同樣制約著叉子圓柏的分布。

第四紀(jì)以來，地球氣候的劇烈變化對當(dāng)今物種的分布格局影響重大(Hewitt， 2004)。末次盛冰期(LGM)是距今最近的極寒冷時(shí)期，全球陸地24%被冰雪覆蓋，兩極冰蓋充分?jǐn)U張，至少覆蓋至30° N中低緯度地區(qū)(王紹武， 2011)。本研究結(jié)果發(fā)現(xiàn)，在冰期叉子圓柏的分布區(qū)多位于30° N左右，分布中心均位于低緯度、低海拔山地和低盆地，較當(dāng)前分布偏南，很可能是由于高海拔和高緯度地區(qū)冰川形成，迫使物種向低緯度、低海拔地區(qū)(1 000～1 500 m)遷移，同時(shí)受到山地地形的影響，形成特殊的微生境，為叉子圓柏及其他植物種提供了冰期避難所得以存活。

全新世中期(MID)全球呈現(xiàn)溫暖濕潤且較穩(wěn)定的特征，平均年均氣溫比LGM時(shí)期高4～5 ℃，海平面上升到與現(xiàn)代相差無幾的高度，降雨量增加(表6)。雖然MID時(shí)期氣溫和降雨量對植物生存有利，但與LGM時(shí)期相比叉子圓柏并未迅速擴(kuò)張其分布面積，而分布區(qū)位置發(fā)生了遷移。通過模型預(yù)測MID時(shí)期溫度升高，冰川開始融化，山地雪線普遍上升了約300～500 m，低海拔區(qū)域由陸地變成湖泊，物種開始向高海拔山地遷移(海拔1 500～2 500 m)，叉子圓柏的分布區(qū)也向北推進(jìn)。MID時(shí)期氣候的另一個(gè)特點(diǎn)是早期氣候濕潤，中期以后較為干旱，全球氣候還處在一個(gè)波動(dòng)激烈的時(shí)期，有B/A暖期，也有YD事件。在YD時(shí)期南歐比北歐更為干旱，孢粉記錄表明希臘甚至比LGM時(shí)期還要干旱，年降水量可能低于150 mm(王紹武， 2011)，導(dǎo)致遷移過程中叉子圓柏分布區(qū)面積微縮(表4)。

表6 各時(shí)期叉子圓柏適生區(qū)的年均氣溫(Bio1)、海拔(Elev)、溫度季節(jié)性變化(Bio4)、年降水量(Bio12)變化范圍及均值Tab.6 The variation range and mean value of annual average temperature(Bio1), elevation(Elev), seasonal variation of temperature(Bio4) and annual precipitation(Bio12) in suitable growing areas of J. sabina in each period

MaxEnt模型模擬結(jié)果表明，當(dāng)前叉子圓柏呈條帶狀分布在歐亞大陸30°～50° N，分布面積明顯擴(kuò)增。從全新世后期到當(dāng)前，環(huán)境因素也發(fā)生了復(fù)雜的變化，基于沉錯(cuò)鉆孔孢粉記錄重建的全新世溫度發(fā)現(xiàn)，溫度最高值出現(xiàn)在初期(中國第四紀(jì)孢粉數(shù)據(jù)庫小組， 2000)，后期(亞大西洋時(shí)期)氣候趨涼爽潮濕，喜冷濕植物群落擴(kuò)張，出現(xiàn)水青岡(Faguslongipetiolata)、云杉(Piceaasperata)、冷杉(Abiesfabri)林。叉子圓柏又作為云杉林下伴生樹種，也應(yīng)大量繁殖擴(kuò)展其分布區(qū)。質(zhì)心分析也表明當(dāng)前相較于MID時(shí)期，叉子圓柏分布區(qū)向西偏南方向偏移擴(kuò)大。

未來由大氣CO2增加引起的全球增溫與MID時(shí)期的增溫效應(yīng)類似(鄭益群等， 2004)。隨著全球氣候變暖，許多動(dòng)植物都有向高緯度和高海拔地區(qū)遷移的趨勢(Sekerciogluetal.， 2008； 方精云等， 2018)?；赗CP2.6和RCP8.5的CO2排放量，測算到2070年叉子圓柏適宜生存區(qū)的海拔在2 500～3 000 m，低緯度、低海拔區(qū)域的叉子圓柏分布面積將大量減少，面臨著消失的風(fēng)險(xiǎn)。MaxEnt模型模擬預(yù)測也表明，未來2070年(RCP2.6和RCP8.5)全球叉子圓柏主要分布在蒙古北部高山草原區(qū)，意大利、瑞士、奧地利境內(nèi)的阿爾卑斯山，俄羅斯西南部與格魯吉亞交界的大、小高加索山，土耳其東北部的大阿勒山，吉爾吉斯斯坦境內(nèi)的天山，中國境內(nèi)的阿爾泰山、天山、祁連山這些高山地帶。分布區(qū)總面積與當(dāng)前相比銳減，不及MID時(shí)期，潛在適生生境喪失或更加破碎化。胡菀等(2020)在研究觀光木(Tsoongiodendronodorum)的潛在地理分布變遷中也得到了相同的結(jié)果，即未來將導(dǎo)致觀光木的生境減少遠(yuǎn)甚于歷史時(shí)期。

3.3 叉子圓柏資源保護(hù)

叉子圓柏是歐亞大陸亞高山地、石質(zhì)山地及沙地的重要生態(tài)樹種。20世紀(jì)以來，雖然人們的生態(tài)保護(hù)意識(shí)日益增強(qiáng)，如： 中國西北部地區(qū)對叉子圓柏的保護(hù)、恢復(fù)、重建工作取得了重大進(jìn)展，部分地區(qū)已建立起省市級(jí)叉子圓柏自然保護(hù)區(qū)，有效地改善了當(dāng)?shù)厣鷳B(tài)環(huán)境，防止了過度放牧的破壞，但自然叉子圓柏灌叢林生境資源與發(fā)展?fàn)顩r仍不容樂觀，天然林分結(jié)構(gòu)退化沒有得到根本遏制。天然叉子圓柏更新雖然擁有匍匐莖萌蘗和種子繁殖2種模式，但其種子庫小、質(zhì)量差、更新環(huán)境喪失等原因限制了許多種群的有性更新(王哲等， 2005； 張國盛等， 2006)。自然更新能力弱，匍匐莖萌蘗成為維持叉子圓柏種群生存的主要繁殖方式，進(jìn)化潛力低，一些種群的遺傳多樣性下降或偏低(路東曄等， 2020)。叉子圓柏為雌雄異株異花授粉樹種，自然分布區(qū)廣，在歷史進(jìn)化過程中不同地理種群、不同發(fā)育階段，都經(jīng)受著環(huán)境變量選擇壓力的不同(劉建全， 2016)，在形態(tài)特征、生理特性、生態(tài)習(xí)性等方面產(chǎn)生了遺傳差異，豐富了種內(nèi)遺傳變異，也形成許多變種(Adamsetal.， 2017； 2018； 路東曄等， 2018； 2020)。因此，在叉子圓柏分布區(qū)范圍內(nèi)，開展種質(zhì)遺傳資源調(diào)查、收集、管理與評(píng)價(jià)等一系列工作，并在此基礎(chǔ)上建立原地和異地遺傳資源保存林，來防止叉子圓柏多樣性的丟失是十分必要的。

基于本研究的結(jié)果，意大利、奧地利、瑞士境內(nèi)的阿爾卑斯山，俄羅斯與格魯吉亞接壤的大、小高加索山脈，中亞的天山、阿爾泰山在4個(gè)時(shí)期內(nèi)均存在穩(wěn)定分布的叉子圓柏種群，側(cè)面說明這些地方極有可能是叉子圓柏的冰期避難所。此外在4個(gè)時(shí)期內(nèi)，叉子圓柏歐亞大陸潛在適生區(qū)的變化特征顯示叉子圓柏適生分布區(qū)為山區(qū)，且在未來分布預(yù)測中，呈現(xiàn)向高緯度高海拔地區(qū)遷移的趨勢，因此要高度重視山區(qū)保護(hù)，確定優(yōu)先保護(hù)單元； 相應(yīng)的沙區(qū)生境破碎化嚴(yán)重，面臨種群消失的風(fēng)險(xiǎn)，要加強(qiáng)沙區(qū)人工林水分條件管理，滿足其天然更新的需求。叉子圓柏在歐亞大陸分布整體而言，西部變化要大于東部，青藏高原的隆起阻擋了西風(fēng)帶，從而導(dǎo)致西部整體受季風(fēng)氣候影響明顯。叉子圓柏多分布于山地，多為社會(huì)經(jīng)濟(jì)條件相對落后的牧區(qū)草場，為促進(jìn)當(dāng)?shù)亟?jīng)濟(jì)發(fā)展，過度放牧、不合理樵采、濫墾、水資源浪費(fèi)等，造成天然叉子圓柏大面積銳減。因此，改善當(dāng)?shù)孛裆鷹l件，創(chuàng)新傳統(tǒng)農(nóng)耕模式，加大科普宣傳實(shí)踐力度，提高保護(hù)意識(shí)，是叉子圓柏保護(hù)的基本措施。

4 結(jié)論

叉子圓柏地理分布受多種環(huán)境變量共同影響，采用氣候和海拔變量進(jìn)行MaxEnt模型模擬，能較完整地?cái)M合出叉子圓柏的實(shí)際分布狀況，反映出叉子圓柏真實(shí)的分布特征?；谀┐问⒈凇⑷率乐衅?、當(dāng)前和未來(2070年)4個(gè)時(shí)期叉子圓柏生境適宜性評(píng)價(jià)，年均氣溫、海拔、溫度季節(jié)性變化是限制叉子圓柏分布的主導(dǎo)環(huán)境變量，山地是其主要生境及避難所，亞洲是其主分布區(qū)。本研究結(jié)果可為叉子圓柏資源管理與發(fā)展、種群恢復(fù)與重建等提供理論依據(jù)。