氣候變化背景下我國長柄扁桃潛在適生區預測*

張淑寧 陳俊興 敖 敦 紅 梅 張雅茜 包福海 王 淋烏云塔娜 白玉娥 包文泉

（1.內蒙古農業大學 呼和浩特 010018；2.中國林業科學研究院經濟林研究所 鄭州 450003）

植被分布格局是植物對包括氣候在內的環境條件長期適應的結果（Fanget al.，2006）。據IPCC 第五次評估報告，過去100 年全球氣溫上升近1 ℃，預測到2100 年仍將升溫1.4～5.8 ℃，這將對陸地生態系統產生巨大影響（Stockeret al.，2014）。研究表明，全球氣候變化對植物的物候、生長和生殖產生顯著影響，進而引起植物分布范圍的遷移，使植物分布由低緯度和低海拔向高緯度和高海拔遷移，或造成其分布范圍縮?。═homaset al.，2004；Matiaset al.，2017）。尤其是瀕危植物，因其遺傳多樣性較低，對氣候變化更敏感，若不采取有力保護措施，將大規模滅絕，全球生態系統會遭到嚴重破壞（Zhuet al.，2012）。研究瀕危植物與氣候因子的關系和預測其適生區受氣候變化的影響，可為其保護利用提供重要理論依據。

生態位模型可基于物種分布數據及其對應的環境參數來獲取物種的生態需求，進而預測適生區域（Schorret al.，2012）。近年來，生態位模型在植物適生區預測、植物保護和入侵植物管理以及探究氣候變化對物種分布格局的影響方面得到了廣泛應用（Elithet al.，2009） 。 目前常用的生態位模型有MaxEnt、CLIMEX、GLM、HABITAT、GARP 和BIOCLIM 模型，其中，MaxEnt 模型預測準確、可靠，穩定性強（Elithet al.，2009），已在互花米草（Spartina alterniflora）（張靜涵等，2020）和朝鮮崖柏（Thuja koraiensis）（蘭雪涵等，2021）等眾多物種中應用，是目前預測物種分布區及其變遷的首選模型。

長柄扁桃（Amygdalus pedunculata）屬薔薇科（Rosaceae）落葉灌木，自然分布于內蒙古、河北、陜西等干旱半干旱區的沙地或山地，是我國特有的木本油料樹種（蘇貴興，1987）。其核仁營養豐富，富含油酸、亞油酸和蛋白質等成分，是食品、化妝品和藥品等的重要原料，具有極高的開發利用價值（郭春會等，2005）。但氣候變化、基礎建設及亂砍濫伐等使長柄扁桃生境嚴重破壞，分布范圍和數量逐年減少，現已被列入我國重點保護的瀕危植物，長柄扁桃野生資源的保護與修復工作已刻不容緩。然而，目前對長柄扁桃的研究主要集中于人工栽植、成分提取、生理生化、遺傳多樣性評價和良種選育等方面（慕曉煒等，2014；孫樹臣等，2020）；而關于長柄扁桃適生區預測及受氣候變化影響的研究較少。本文將優化MaxEnt 模型的參數，確定長柄扁桃當前適生區域，推測末次間冰期、末次盛冰期和全新世中期3 個歷史時期以及未來氣候情景SSP126、SSP245 和SSP585 下的長柄扁桃適生區變化趨勢，以期為長柄扁桃資源的保護與人工栽植提供科學依據。

1 研究方法

1.1 地理分布數據的收集

在2015—2019 年，對長柄扁桃自然分布區進行了系統調查，從內蒙古、陜西、河北、山西等地共獲取61 個分布數據；結合相關研究文獻、國家標本館、全球生物多樣性信息網絡數據庫，獲取長柄扁桃187 個分布數據，對其刪除錯誤及重復分布點，最終獲得148 個分布點（圖1）。

圖1 我國長柄扁桃地理分布點信息Fig.1 Distribution information of A.pedunculata in China

1.2 環境因子的獲取及篩選

通過分析長柄扁桃分布區的氣候、地形和土壤因子，初步篩選出34 個環境因子作為建模變量，包括19個生物氣候因子、14 個土壤因子和1 個地形因子（海拔）；其中，氣候因子數據取自世界氣候數據庫（www.worldclim.org）的末次間冰期（last inter glacial）、末次盛冰期（last glacial maximum）、全新世中期（mid holocene）、當前（1960—1990 年）、2050s（2041—2060年）和2090s（2081—2100 年）6 個時期，包括最暖季平均氣溫、月平均晝夜氣溫差、最熱月最高氣溫、最濕月降水量、年均氣溫、最濕季平均氣溫、最干季平均氣溫、氣溫變異系數、氣溫年較差、年降水量、最暖季降水量、等溫性、最冷季平均氣溫、最干季降水量、最冷季降水量、最冷月最低氣溫、降水量變異系數、最濕季降水量和最干月降水量。土壤因子和地形因子數據來源于世界土壤數據庫的中國土壤數據集（http://www.fao.org/faostat/en/#data.）；土壤因子包括0～30 cm 土層中的礫石體積百分比、陽離子交換量、土壤質地類型、砂礫百分比、黏粒百分比、粉砂百分比、黏性成分陽離子交換量、土壤密度、pH 值、交換性鹽基總和、鹽基飽和度、可交換鈉鹽百分比、電導率和有機碳含量。

將氣候因子和海拔因子數據按1∶100 萬中國行政地圖矢量圖進行剪裁，并利用ArcGis 將其轉換為ASC 格式；運用ArcGis 提取MU_GLOBAL 圖層的表層土壤因子柵格圖層，將其轉為ASC 格式；最后將各環境因子的空間分辨率統一設置為30″，使環境因子數據柵格圖像范圍和像元大小完全一致。

本研究涉及的地圖數據下載于自然資源部官網（http://www.mnr.gov.cn/）。以社會經濟共享路徑（shared socio-economic pathways, SSP）的低濃度（SSP126）、中濃度（SSP245）和高濃度（SSP585）排放情景代表未來氣候系統模式，假設未來50 年土壤和地形因子不發生明顯變化（曹麗格等，2012；翁宇威等，2020），對長柄扁桃進行適生區預測。為避免各環境因子間的共線性而導致模型預測結果過度擬合，本研究利用方差膨脹因子VIF 和Pearson 相關性檢驗，篩選相關性小于0.8 且VIF 值小于10 的因子（張天蛟等，2017）；同時，對相關系數大于0.8 的因子保留其生態學意義更重要的因子（Yeet al.，2020）。最終，篩選出8 個氣候因子、4 個土壤因子和1 個地形因子（表1）。

表1 預測長柄扁桃地理分布的環境因子Tab.1 The environmental variables for prediction of geographical distribution of A.pedunculata

1.3 物種分布模型的構建與評價

將長柄扁桃148 個分布數據和所篩選的13 個環境因子數據導入MaxEnt 模型，隨機選取75%的樣點作為訓練數據集構建模型，將剩余25%的樣點作為測試數據集驗證模型，設置10 次重復；并用Jackknife 評估各環境因子的權重，結合其貢獻率確定主導環境因子。

為準確評價模型的可靠性，采用刀切法分析環境因子的貢獻率，利用受試者工作特征曲線ROC 下的面積AUC 值評估所建模型精度（Fandet al.，2014）。通常認為AUC 值越大預測越精確，其中，AUC 值在0.5～0.7 時表示預測較差，0.8～0.9 時表示較好，而0.9～1.0 時則表示非常好（Freitaset al.，2019）。

1.4 模型的優化

應用ENMeval 程序包優化MaxEnt 模型，將調控倍頻設置為0.5～8.0，每次間隔0.5，共設置16 種調控倍頻；并采用9 個特征組合，即：L、LQ、H、LQH、LQHP、LQHPT、QHP、QHPT 和HPT。ENMeval 程序包將上述144 種參數組合進行測試，最終依據Akaike信息量準則（AICc），以delta.AICc 作為測試模型的復雜度和擬合度，選取參數組合，其中delta.AICc 值越低，表示模型預測越準確（李安等，2020；趙光華等，2021）。

1.5 數據處理

根據MaxEnt 模型預測的長柄扁桃適宜性閾值，結合自然斷點法劃分長柄扁桃生境適宜性指數，并運用ArcGis 軟件將其可視化。將長柄扁桃生境適宜性劃分為4 個等級，即不適生區（0～0.1）、較不適生區（0.1～0.3）、一般適生區（0.3～0.5）和高度適生區（0.5～1.0）。比較不同時期的適生區變化，獲得氣候變化背景下長柄扁桃未來空間分布格局變化；調用SDMToolbox工具，計算長柄扁桃當前和未來氣候情景下的適生區質心位置，使用dists 函數計算其質心遷移距離，以質心位置的變化表示其適生區空間分布的遷移方向。用ArcGis 將當前和未來適生區數據轉換為柵格數據并對其進行模糊疊加，得到長柄扁桃潛在地理分布圖層，再進行重新劃分適生等級，獲取長柄扁桃適生區劃圖。

2 結果與分析

2.1 MaxEnt 模型的優化與評價

將148 個分布點和13 個環境因子導入MaxEnt 模型，預測長柄扁桃的潛在分布區。在MaxEnt 默認參數下（RM=1，FC=LQHPT），delta.AICc 值為225.87；而當RM=1.5、FC=LQHPT 時，delta.AICc 值為0；依據AIC 信息準則，此時模型最優（圖2），模型的10 次預測訓練集和測試集AUC 值也均大于0.90，均值為0.967，其訓練遺漏率較低（圖3），說明模型預測較好，結果可靠。

圖2 不同參數設置下MaxEnt 模型的評估結果Fig.2 Evaluation results of MaxEnt model under different settings

圖3 最優參數設置下MaxEnt 模型ROC 曲線和遺漏率曲線Fig.3 ROC curve and Omission curve of MaxEnt model under optimal parameter setting

2.2 長柄扁桃在我國當前潛在地理分布

基于MaxEnt 模型，構建長柄扁桃當前適生分布區（圖4），其當前總適生區（一般適生區和高度適生區）面積為447 274.3 hm2，主要分布于內蒙古高原和黃土高原地區，其中高度適生區為219 166.7 hm2，占我國國土總面積的2.3%（表2）。當前我國長柄扁桃高度適生區集中于內蒙古的“呼包鄂”地區，此外，赤峰、烏蘭察布和陜西榆林也少有分布。一般適生區則以高度適生區為中心向四周分布，從西至東依次分布于甘肅、寧夏、內蒙古、陜西、山西和河北等省區。

表2 長柄扁桃不同時期的適生區面積①Tab.2 The suitable living area for A.pedunculata under different climates scenarios/years

圖4 基于MaxEnt 模型預測的我國長柄扁桃當前適生區分布Fig.4 Distribution of current suitable areas of A.pedunculata in China based on MaxEnt

2.3 長柄扁桃在我國歷史潛在地理分布

通過對我國長柄扁桃歷史潛在地理分布的推測發現，在末次間冰期，長柄扁桃適生分布區較廣，西可達西藏，東可達山東，該時期長柄扁桃的高度適生區面積為620 431.9 hm2，為當前高度適生區的2.83 倍（圖5，表2）；而在末次盛冰期，長柄扁桃的適生區嚴重收縮，空間喪失率高達81.44%（圖6）；并且這種收縮趨勢延續至全新世中期，此時高度適生區面積僅為末次盛冰期的44.8%（圖6，表3）；但全新世中期有內蒙古東部的通遼，山西太原以北，以及北京、河北等較高緯度地區已成為新增適生區。因此，我國長柄扁桃從末次間冰期至末次盛冰期，其空間分布格局變化顯著，主要由高緯度和低緯度向中高緯度收縮；而從末次盛冰期至全新世中期，適生區則由東向西進一步收縮（圖5）。

表3 長柄扁桃不同時期適生區的變化①Tab.3 Change of A.pedunculata suitable area in different periods

圖5 歷史、當前及未來不同氣候情景下長柄扁桃在中國的適生區分布Fig.5 Distribution of suitable areas of A.pedunculata under historical, current and future climate scenarios in China SSP：社會經濟共享路徑Shared socio-economic pathways.

圖6 歷史、當前及未來不同氣候情景下長柄扁桃在中國的空間分布格局變化Fig.6 The spatial pattern of A.pedunculata in China under historical, current and future climate scenarios SSP：社會經濟共享路徑Shared socio-economic pathways.

2.4 未來氣候情景下的長柄扁桃適生區預測

在未來不同氣候情景下長柄扁桃適生區除在2050s-SSP126 情景下較當前有增加外（增加42.18%），其余未來氣候情景下均有不同程度的縮減（圖6）。在2050s-SSP126 情景下，適生區喪失面積最少（92 100 hm2），新增率最大（63.32%），新增區可分為東、西部兩區，東部區為內蒙古錫林郭勒、赤峰、通遼以及內蒙古與吉林省交界區域，而西部區為甘肅、寧夏與內蒙古交界帶（表3）。該情景下，長柄扁桃適生區面積約為635 954.8 hm2，是當前適生區的1.42 倍，說明此時將有大面積非適生區轉變為適生區。

2090s-SSP245 情景下，長柄扁桃適生區面積約為345 833.3 hm2，較當前適生區收縮了近1/3，適生區新增率最?。?.33%），新增區域分布于天津、河北、山西等地。在2090s-SSP126 情景下，長柄扁桃適生區喪失面積最大（241 200 hm2），喪失率高達53.83%，而新增率僅為9.15%；此時，吉林、遼寧、通遼及赤峰地區的適生區將大面積喪失，僅有內蒙古錫林郭勒以及陜西延安的部分地區新增為適生區。

通過預測未來適生區變化可知，長柄扁桃對不同氣候情景下氣候變化的響應存在一定差異（圖5）。在SSP126 情景下，高度適生區面積具有先增后減的趨勢；而在SSP245 和SSP585 情景下，高度適生區面積均隨時間呈縮減趨勢，說明長柄扁桃對氣候變化的響應較敏感。

2.5 長柄扁桃適生區質心遷移路線

通過對空間分布格局變化的分析發現，在未來不同氣候情景下，長柄扁桃適生區的遷移距離存在一定差異，但其遷移趨勢基本一致（圖7，表4）。當前長柄扁桃適生區質心位于內蒙古烏蘭察布市涼城縣（112.52°E，40.55°N）；而當氣候情景為SSP126-2090s時，適生區質心向西北遷移78 638 m，此時適生區質心位于呼和浩特市土左旗（110.91°E，40.95°N）；當氣候情景為SSP585-2090s 時，適生區質心向西遷移340 468 m，遷移至內蒙古自治區烏拉特前旗（108.50°E，40.49°N）。在未來氣候變化背景下，適生區質心整體往西北遷移，并且隨著溫室氣體排放濃度的加大，遷移距離具有向西擴大的趨勢。

表4 歷史、當前和未來不同氣候情景下長柄扁桃適生區質心的遷移距離①Tab.4 Centroid migration distance of suitable areas of A.pedunculata under historical, current and future climate scenarios/m

2.6 環境因子分析

刀切法分析結果顯示（表5），最暖季度降水量、降水量變異系數、氣溫變異系數、年均溫、表層土壤鹽基飽和度、年降水量、最濕季度平均氣溫等7 個環境因子對長柄扁桃地理分布的貢獻率大于4%，其累計貢獻率和重要置換值分別為87.9%和87.2，其中，最暖季度降水量的貢獻率最大（22.7%），而降水量變異系數的貢獻率次之（20.3%），表明上述7 個因子是影響長柄扁桃適生分布的主要環境因子，是構建模型的最重要環境變量。

表5 參與建模的環境因子貢獻率及置換重要值Tab.5 Percentage contribution and permutation importance of environment variables for A.pedunculata in the MaxEnt model

依據環境因子響應曲線可知（圖8），長柄扁桃對6 個主要環境因子的響應曲線均呈單峰（圖8），表明所建模型的預測較為準確，其中，年均溫、氣溫變異系數、降水量變異系數、年降水量、最暖季度降水量和表層土壤鹽基飽和度對長柄扁桃分布的高度適宜范圍分別為2.3～9.8 ℃、1 090～1 380、97%～119%/130%～150%、180～500 mm、148～296 mm 和92.9%～99.4%。

圖8 長柄扁桃對6 個主要環境因子的響應曲線Fig.8 Response curves of A.pedunculata to five major environmental factors

3 討論

3.1 優化模型的長柄扁桃適生區預測

本文基于氣候、土壤和地形因子，參照最新的Checkerboard2 法，調用R 語言ENMeval 程序包，優化了MaxEnt 模型的參數。模型參數的優化可限制背景數據與校準區域，使預測分布區覆蓋當前的分布點，改變正則化水平，提高模型預測與實際分布區的擬合度，并通過檢查地理預測圖對模型準確性進行衡量（Phillipset al.，2017）。因此，優化后的模型可有效降低其復雜性，能夠顯著提高模型預測的準確性（Phillipset al.，2017；趙光華等，2021）。屈振江等（2016）和葉學敏等（2019）曾采用MaxEnt 模型預測我國紅富士蘋果（Malus pumila）和南酸棗（Choerospondias axillaris）的適生區研究，但均未對MaxEnt 模型進行參數優化，導致預測誤差較大，較難解釋其研究結果。本文采用優化參數后的MaxEnt 模型，將默認參數下MaxEnt 模型的delta.AICc 值由225.87 降低至0，形成最優模型，降低了模型的復雜度，顯著提高了模型預測準確性，其預測結果也與課題組前期調查結果相似（包文泉等，2018）。這與Rong 等（2019）預測青海云杉（Picea crassifolia）適生區時提出的MaxEnt 模型參數優化可使預測結果更具說服力的結論一致。可見，采用MaxEnt 模型預測物種適生區，其參數優化至關重要，是決定模型預測準確的關鍵因素之一。

長柄扁桃根系發達，小灌木，具有極強的抗旱、抗寒和耐貧瘠性，喜生于干旱及荒漠區的固定沙地、石礫質陽坡、山麓等地（Heet al.，2021）。通過本研究團隊前期對長柄扁桃的野外調查可知，我國長柄扁桃當前實際分布區主要在內蒙古、陜西、山西、河北、寧夏和甘肅（包文泉等，2018），而本研究應用優化后的MaxEnt 模型所預測長柄扁桃適生區也恰好主要位于內蒙古、陜西、河北、甘肅、山西等地，預測與實際分布區基本一致，表明MaxEnt 模型可用于長柄扁桃分布區預測，并且優化后的MaxEnt 模型預測可靠、準確。但由于本研究僅采用了MaxEnt 模型，預測結果可能仍存在精度不夠高或一些不確定性，因此后續還應采用Bioclim、Domain、Biomod2 和Garp 等模型同MaxEnt相結合，以組合模型預測加以驗證本研究結果，所得結果可能會更精確。

3.2 環境因子對長柄扁桃適生分布的影響

分析表明，最暖季度降水量、降水量變異系數、氣溫變異系數、年均溫、表層土壤鹽基飽和度、年降水量和最濕季度平均氣溫是影響長柄扁桃分布的主要環境因子。環境因子的響應曲線顯示，長柄扁桃適生區的年均溫為2.3～9.8 ℃，降水量變異系數為97%～119%/130%～150%，最暖季度降水量為148～296 mm，這與內蒙古地區長柄扁桃適生區的年均溫與降水特征相符（馬梓策等，2019），也與長柄扁桃分布廣泛的陜西榆林地區的年均溫與降水量變異系數特點基本相符（劉曉瓊等，2017）。前期調研發現，長柄扁桃年周期中不同物候期對土壤水分含量及環境濕度的要求存在較大差異，在果實快速生長期與花芽分化期要求降雨量較高，而花期及營養生長期對降水要求相對較低（包文泉等，2016），這與本研究中降水量變異系數和最暖季度降水量等是影響長柄扁桃分布的主要環境因子這一結果一致。此外，本研究篩選出4 個土壤因子進行建模，但僅有表層土壤鹽基飽和度的貢獻率大于4%（8.4%），這與長柄扁桃根系發達、抗旱、耐瘠薄，對土壤條件要求較低，能在荒山、荒地、沙地、巖石等生態脆弱區正常生長的生物學特性相一致（烏云塔娜等，2014）。

本研究表明，長柄扁桃主要適生區分布在我國北方地區，溫度和降水是影響長柄扁桃適生的主要因子，其貢獻值顯著大于土壤和地形因子。但由于土壤和地形因子是較復雜的環境因子，并且二者還可對溫度和降水進行再分配，因此，溫度、降水、土壤、地形因素對長柄扁桃適生分布的影響均不能忽視，針對四者間存在何種關系以及具體影響機制是什么，后期還需進一步研究和分析。

3.3 氣候變暖對長柄扁桃適生分布的影響

通過對長柄扁桃歷史、當前和未來潛在分布區的模擬表明，氣候變暖對長柄扁桃適生區分布的影響顯著。從長柄扁桃歷史分布區的推測可知，在末次間冰期，長柄扁桃分布于我國北方大部分地區，面積約為當前適生區的1.5 倍；而在末次盛冰期，長柄扁桃適生區從東部向西部嚴重收縮，集中于黃土高原和華北平原地區，表明末次盛冰期的氣候變化嚴重縮小了長柄扁桃的分布范圍，使其僅能在中高緯度地區的適生區得以留存；此后，全新世中期是距今最近的一次大暖期，氣候變暖使長柄扁桃適生區較末次盛冰期縮減近一半，并向高緯度區遷移，新增區集中于內蒙古東北部。這一結果與眾多植物歷史分布區研究相符，物種的適生區范圍與氣溫升高呈負相關（鄭卓等，2013），在氣溫升高時向高緯度地區遷移（劉金陵等，2004；白偉寧等，2014）。本研究中長柄扁桃在其歷史上呈動態分布，從末次間冰期至末次盛冰期，長柄扁桃的分布格局變化最顯著，由西南往東北收縮；而從末次盛冰期至全新世中期，則由東向西收縮；可見，冰期與間冰期交替的歷史氣候事件促進了長柄扁桃當前分布格局的形成，但其影響程度和作用機制還尚未清楚，后期還需更深入探究。

與當前相比，未來高溫高濕環境使長柄扁桃適生分布區不斷向西北擴張，但擴張面積遠低于喪失面積，表明長柄扁桃對氣候變化響應較敏感，尤其在高濃度溫室氣體情景下更敏感。這一結果與絕大數北方落葉林木對氣候變化的響應一致（Zhaoet al.，2021）。

長柄扁桃空間分布格局的變化分析也表明，隨著氣候變暖，適生空間有向西北遷移的趨勢；并且在溫室氣體的低、中排放濃度下，適生區質心的遷移距離較?。欢跍厥覛怏w高排放濃度下，適生區幾何中心的遷移距離較大。這與溫帶樹種適生分布區在氣候變暖下往高緯度、高海拔區遷移的趨勢相符（Yeet al.，2020；張殷波等，2019）。因此，未來由大氣CO2濃度增加引起的全球氣候變化同樣是影響長柄扁桃分布的主要因素，尤其是低緯度、低海拔區的長柄扁桃分布面積較當前銳減，具有適生生境喪失的可能。

4 結論

本研究對MaxEnt 模型默認參數進行優化，得到最優模型的特征組合為線性、二次型、片段化、乘積型與閾值性特征，其調控倍率為1.5。優化后模型的訓練遺漏率和復雜度較低，預測結果與實際分布一致，準確可靠。當前我國長柄扁桃適生分布區主要在內蒙古高原和黃土高原地區，屬于中高緯度區。在歷史、當前和未來氣候情景下，長柄扁桃對氣候變化的響應較為敏感；在未來不同氣候變化情景下，長柄扁桃適生區有縮減的趨勢，由低海拔和低緯度向中高緯度和高海拔地區遷移，尤其在溫室氣體高排放濃度下遷移距離更大。

在全球氣候增溫增濕背景下，長柄扁桃在吉林、遼寧、赤峰、通遼等高緯度地區的適生區將大面積消失，而在中高緯度和高海拔（內蒙古呼和浩特、鄂爾多斯、錫林郭勒，陜西延安）有新增適生區出現。隨著未來溫室氣體排放濃度增加，長柄扁桃當前適生區的東南部區可能不再適宜。因此，長柄扁桃產業的今后發展，應選擇當前和未來均適宜的區域，如重點選擇內蒙古的‘呼包鄂’，以及內蒙古、陜西、山西、河北與寧夏交接區域，可避免未來氣候變暖帶來的損失。