末次間冰期以來及未來氣候情景下紅皮云杉適生分布區模擬*

王 亞 王軍輝 王福德 劉逸夫 譚燦燦 袁艷超 聶 穩劉建鋒 常二梅 賈子瑞

（1.林木遺傳育種國家重點實驗室 中國林業科學研究院林業研究所 北京 100091；2.中國林業科學研究院森林生態環境與自然保護研究所 北京 100091；3.黑龍江省林業科學研究所 哈爾濱 150081）

植物與氣候的相互關系一直是植物學和生態學研究的熱點（安芷生等，2005），氣候是決定物種分布的主要因素之一，物種分布格局變化能直觀反映氣候在不同時期的變化?；谖锓N分布模型（species distribution models, SDMs）和物種分布點的經緯度數據與環境因子，通過數學方法可得到各環境因子的貢獻率，構建物種潛在分布范圍，并以此推測未來和過去時期的適生分布區。物種分布預測在生態學、生物地理學等領域十分重要，已在土地利用變化、未來氣候變化情景（Faleiroet al.，2013）中得到廣泛應用，用來研究生物多樣性損失（Bertrandet al.，2012）、瀕危物種管理和保護（Guoet al.，2016）、生物熱點區域確定（Walterset al., 2019）等。物種分布預測的可靠性受環境變量選取和模型參數的影響（Fourcadeet al.，2018;Sillero, 2011），環境變量之間的相關性會引入冗余信息并導致模擬結果偏差（Verbruggenet al.，2013）。因此，去除多余變量可消除變量間重線性關系，利于精確模擬物種的可能適生分布，并為生態保護等理論問題提供重要參考數據（Ashcroftet al.，2011）。

全球氣候變化背景下的物種適生分布區不斷變化。第四紀末次間冰期以來，全球氣候頻繁振蕩，極大影響物種種群動態變化（張愛平等，2018）。物種經歷盛冰期和間冰期的不斷回旋，逐步形成當前的地理分布格局（鄒旭等，2018）。物種分布模型研究表明，一些物種的分化和形成與歷史事件相關。對瀕危物種奇瓦瓦云杉（Picea chihuahuana）的最大熵模型分析發現，自末次盛冰期以來其適生區一直減少（Pinedo-Alvarezet al.，2019）；由于森林碎片化和氣候變化，其孤立種群在未來面臨消失的危險（Aguilar-Sotoet al.，2015）。對紫果云杉（Picea purpurea）、麗江云杉（Picea likiangensis）、青杄（Picea wilsonii）的潛在分布區預測發現，等溫性、最暖季降水量和年均溫是這3 個物種生態位分化的主要因子（張愛平等，2018），而且這3個物種在末次盛冰期發生擴張。對云杉屬34 個物種的生態位分析表明，氣候對云杉形態和物種分化的影響是由地理隔離、生態位保守性和生態適應等許多生物和非生物因素介導的，溫度和降水分別是云杉系統發育和形態初級分化的主要驅動因素（Wanget al.，2017）。對東亞孑遺植物中全新世、末次盛冰期和未來的分布區預測結果發現，孑遺植物物種豐度和氣候避難所相對應，中國西南部和越南北部是許多第三紀孑遺植物的長期避難所（Tanget al.，2018）。

紅皮云杉（Picea koraiensis）又名臭松，為松科（Pinaceae）云杉屬（Picea）植物，是中國東北重要的針葉樹種和古老的孑遺樹種，主要分布于大興安嶺北部、小興安嶺、長白山、完達山、張廣才嶺、老爺嶺等地區，地理坐標范圍在41°40'—53°15'N、122°—131°E 之間，垂直分布范圍為250～1 100 m，屬當地主要造林與用材樹種。紅皮云杉分布于針闊混交林中，并處于針葉林到闊葉林的過渡帶，是研究植物分布響應氣候變化的理想植物之一。紅皮云杉具有涵養水源、防風固沙、水土保持等生態效益以及較高的經濟價值，在木材生產和生態服務方面均有重要作用。目前，關于紅皮云杉的研究主要涉及生態學系統特征和生物量（張含國等，2003）、遺傳多樣性（王秋玉等， 2004）、礦質元素含量（姜興林等，2002）等，不同時期氣候變化下的紅皮云杉地理分布格局變化尚未見報道。鑒于此，本研究基于最大熵模型（Merowet al.，2013），結合GIS 技術，利用末次間冰期（last interglacial, LIG, 120～140 ka BP）、末次盛冰期（last glacial maximum, LGM,21 ka BP）、全新世中期（mid-holocene, MH, ca.6 000 years ago）、現代（current, 1960—1990 年）和未來2 個時期（2050s 和2070s）的氣候數據和物種分布點數據，推測紅皮云杉過去、現代和未來的適生分布區，識別主要氣候影響因子，揭示氣候對紅皮云杉地理分布格局形成的影響，推測種群遷移的歷史動態和生物避難所，并預測未來分布變化趨勢，以期為該物種的資源管理和多樣性保護提供科學依據。

1 材料與方法

1.1 物種分布數據收集

通過3 種方法收集紅皮云杉地理分布數據：第一，實地調查收集紅皮云杉的分布地點記錄，每個調查點之間的距離不小于20 km，記錄分布點的位置名稱、經緯度和海拔；第二，檢索查詢相關文獻并記錄紅皮云杉分布的地理信息；第三，搜索中國虛擬植物標本館（http://www.cvh.ac.cn/）和全球生物多樣性信息網絡（https://www.gbif.org/）等數據庫，收集紅皮云杉的分布點數據。匯總分布點信息，刪除地理信息不完整和重復的分布記錄，最終確定59 個紅皮云杉分布點。根據MaxEnt 模型的輸入文件格式，將物種分布點數據按照物種名、分布點經度和緯度的順序進行整理。

基于ENMTools v1.4 軟件（Warrenet al.，2010），使用Trim duplicate occurrences 函數對整理好的物種分布數據進行冗余篩選，以避免地理聚集導致的空間自相關。最后共獲得55 個紅皮云杉有效分布點數據，用來構建模型。從自然資源部網站（http://www.mnr.gov.cn/）下載中國地圖，使用ArcGIS v10.4.1 軟件將紅皮云杉分布點數據在地圖上展示出來（圖1）。

圖1 當前紅皮云杉在中國的分布記錄Fig.1 Current distribution of P.koraiensis based on occurrence records

1.2 環境變量

云杉屬植物生長發育和繁殖受多種氣候因子影響，如溫度、降水等。本研究選取與溫度和降水相關的19 個環境變量進行最大熵模型分析,其中，現代19個與溫度和降水有關的生物氣候變量（Bio1—Bio19）來自WoldClim version 1.4（Ficket al.，2017），精度為2.5 arc-minutes（弧分）。為避免因高度相關的氣候變量之間的相互影響而導致模擬結果過度擬合，利用MaxEnt v3.3.4 的刀切法（jackknife）將19 個氣候變量按貢獻率大小排序，篩選出貢獻率大于1 的變量。利用ENMTools 1.4 在19 個氣候變量之間進行相關系數分析，一對相關系數大于|0.9|的氣候變量被視為非獨立變量。對每對顯著相關的變量，只選1 個有較大貢獻的變量用于分析。

1.3 過去和未來氣候數據

過去和未來時期的氣候變量數據從Worldclim version 1.4（Ficket al.，2017）下載，其中LGM 和MH 時期的數據精度為2.5 arc-minutes（弧分），LIG 時期的數據精度為30 arc-seconds（弧秒）。大氣環流模型選擇CCSM4（community climate system model four），CCSM4 是能夠模擬過去和未來氣候變化對植物分布影響的最有效的氣候模型之一（Abdelaalet al.，2019），已成功用于類似研究（Araret al.，2020; Liet al.，2020）。利用ArcGIS v10.4.1 將所有柵格文件轉換成ASCII 格式，將LIG 數據重抽樣為2.5 arc-minutes。下載未來2050s 和2070s 兩個時間段的氣候數據，選擇RCP4.5（representative concentration pathway）和RCP8.5 兩種排放氣候場景，用于模擬紅皮云杉未來地理分布格局。

1.4 最優模型選擇

由于要素類型（feature class，FC）和正則化乘數（regularization multiplier，RM）2 個變量對模型結果影響較大，為獲得更精準的適生區預測結果，本研究利用MaxEnt v3.3.4 軟件、R 3.6.2 程序和kuenm 程序包（Coboset al.，2019）對這2 個變量進行優化。首先利用R 軟件的kuenm 程序包將每個物種csv 格式的數據點分布文件拆分，70%的地理分布點用于訓練集構建、30%用于測試集驗證?；诓鸱趾蟮牡乩矸植键c數據和篩選后的環境變量，利用R 語言的kuenm 軟件包從MaxEnt 模型的40 個正則化常數（0.1～4.0）和10種要素類型（L、Q、P、H 的2 種組合和3 種組合）中選擇最優模型組合，要素類型含義如下：L（linear features）為線性，Q（quadratic features）為二次型，P（product features）為乘積型，H（hinge features）為片段化。利用kuenm 程序包逐一檢驗上述400 組參數，根據最佳模型選擇標準（Coboset al.，2019），選擇統計學上顯著、遺漏率（omission rate）低于閾值（0.05）且ΔAICc（delta Akaike information criterion）值不高于2 的模型，具有ΔAICc 最小值的模型被作為最優模型。

1.5 紅皮云杉潛在分布預測

應用MaxEnt 最大熵模型預測紅皮云杉在過去、當前和未來氣候情景下的潛在分布。將選定的物種分布點數據和環境變量導入MaxEnt v3.3.4 軟件。采用10 倍交叉驗證法對模型進行訓練和驗證，輸出格式為Logistic，輸出文件類型為asc，最大迭代次數設置為5 000，線程設置為10，重復次數為10 次。MaxEnt建模完成后，采用受試者工作特征曲線（receiver operating characteristic，ROC）評估模型，并使用曲線下面積（area under the ROC curve，AUC）測量模型可信度（Alloucheet al.，2006），AUC＞0.8 時表明物種適生區預測精度較好，AUC＞0.9 時表明預測精度非常好。本研究中，使用ArcGIS v10.4.1 軟件創建當前、過去（LIG、LGM 和MH）和未來（2050s 和2070s）的適生區預測圖。MaxEnt 生成結果導入ArcGIS v10.4.1 軟件，并以中國地圖為底層進行分析。使用自然斷點法劃分紅皮云杉的適生區等級：非適生區（0～0.1）、低適生區（0.1～0.3）、中適生區（0.3～0.5）和高適生區（0.5～1），計算各時期的適生面積。

1.6 適生分布區質心轉移

為明確紅皮云杉適生分布區的變化趨勢，在ArcGIS v10.4.1 軟件中加載SDMToolbox （Brown，2014）工具箱，利用Centroid changes command 命令，計算不同時期的分布質心。將這些質心在ArcGIS v10.4.1 軟件中利用矢量箭頭投影，表示適生區變化的遷移距離和方向。

2 結果與分析

2.1 模型選擇和精度評價

為精準預測適生區，本研究篩選出6 個環境變量，并選擇最優的要素類型和正則化乘數，用于最大熵模型分析。6 個環境變量為年均氣溫（Bio1）、晝夜溫差與年溫差比值（Bio3）、氣溫季節性變動系數（Bio4）、降水量季節性變化（Bio15）、最干季度降水量（Bio17）、最暖季度降水量（Bio18），最優模型要素類型為PH（product features, hinge features），最優正則化乘數為2.3。10 次重復訓練數據集和測試數據集的平均AUC 分別為為0.995 8 和0.995 7，均接近于1，說明模型有較高穩定性和精度，模擬結果可信度高（圖2）。

圖2 基于MaxEnt 模型預測紅皮云杉分布的受試者工作特征曲線Fig.2 The ROC curve predicting the distribution of P.koraiensis based on MaxEnt

2.2 環境變量篩選

根據MaxEnt 模型的輸出環境因子貢獻率和刀切法分析結果，比較各環境變量對紅皮云杉潛在分布區變化的貢獻率（表1、圖3），最暖季度降水量（Bio18）、氣溫季節性變動系數（Bio4）和降水量季節性變化（Bio15）為3 個主要變量，貢獻率分別達40.6%、28.9%和24.6%，年均氣溫（Bio1）、最干季度降水量（Bio17）和等溫性（晝夜溫差與年溫差比值）（Bio3）的貢獻率分別為4%、1.5%和0.4%。

表1 影響紅皮云杉適生分布的環境變量貢獻百分率Tab.1 Contribution percentage of environmental variables affecting the suitable distribution of P.koraiensis

圖3 環境變量刀切法檢驗Fig.3 Jackknife tests for environmental variables

重要變量響應曲線（圖4）顯示出對紅皮云杉分布影響較大的變量（貢獻率＞2）適生范圍：最暖季度降水量（Bio18）為332.9～527.9 mm，氣溫季節性變動系數（Bio4）為129.4～173.5，降水量季節性變化（Bio15）為87.6～111.2，年均氣溫（Bio1）為-1.4～5.6 ℃。

圖4 紅皮云杉對重要環境變量的響應曲線Fig.4 Response curves of P.koraiensis to important environmental variables

2.3 目前、過去和未來適生分布區

紅皮云杉當前適生分布區包括中國的大興安嶺、小興安嶺、長白山脈及其支脈等地區，涉及黑龍江、吉林、遼寧等省，其中高適生區的分布范圍與實際分布相符（圖5）。紅皮云杉當前的總適生區面積達133.58×104km2，其中，低、中、高適生分布區面積分別為40.93×104、42.36×104和50.29×104km2（表2）。

表2 紅皮云杉在多種氣候條件下的適生面積①Tab.2 Prediction of suitable area of P.koraiensis under different climate scenarios

圖5 當代、過去和未來氣候情景下紅皮云杉在中國適生分布區Fig.5 Suitable distribution areas of P.koraiensis in China under contemporary, past and future climate conditions

歷史氣候條件下，紅皮云杉在末次間冰期（LIG）、末次盛冰期（LGM）和中全新世（MH）的適生分布區面積分別為目前的102.59%，105.20% 和100.53%；高適生區和總適生區的變化趨勢一致，末次盛冰期種群擴張，全新世中期及以后種群持續收縮（圖5、表2）。

未來氣候情景下，紅皮云杉的總適生區和高適生區面積均呈減少趨勢，高適生區邊緣呈片段化分布，而小興安嶺腹地、長白山腹地、老爺嶺以及完達山等地的種群較為穩定（圖5、表2）。紅皮云杉在未來氣候情景中2 種溫室氣體排放模式下（RCP4.5-2050s、RCP4.5-2070s、RCP8.5-2050s 和RCP8.5-2070s）的總適生區面積分別為當前的99.06%、96.88%、98.05%和95.08%。在高排放（RCP8.5）模式下，紅皮云杉適生區面積縮減更大。RCP4.5 模型預測2050s 時期高適生區縮減的區域在大興安嶺北部邊緣（塔河縣）到中段（嫩江市）、長白山西部（輝南縣和撫松縣）；2070s 時期這2 個區域持續縮減，大興安嶺北部（蒙克山地區）和中部（高峰地區）、張廣才嶺北部（勃利縣和樺南縣）和東部（海林市）的高適生區消失。RCP8.5 模型預測2050s 時期高適生區縮減的區域在張廣才嶺北部（勃利縣和依蘭縣）和東部（海林市）；2070s 時期，紅皮云杉高適生區邊緣持續碎片化，長白山西部（輝南縣和撫松縣）的高適生區消失。

2.4 適生分布區質心轉移

通過重建現代、過去和未來不同時期氣候條件下紅皮云杉的適生區質心坐標，分析質心變化。紅皮云杉LIG 時期的分布質心坐標在46°14'27.38''N,127°1'46.92''E；LGM 時期分布質心（46°10'54.01''N,127°14'49.34''E）向東南轉移21.62 km；MH 時期分布質心（46°15'0.44''N, 127°18'54.83''E）向東北轉移8.96 km；當前時期分布質心（46°6'14.27''N, 126°14'28.81''E）向西南轉移104.43 km。未來氣候條件下，當氣候情景為中排放濃度（RCP4.5）時，2050s年代分布質心（46°25'2.47''N,127°17'54.03''E）向東北轉移105.72 km；2070s 年代分布質心（46°15'33.81''N, 127°3'15.60''E）向西南轉移27.4 km；當氣候情景為高排放濃度（RCP8.5）時，2050s年代分布質心（45°49'27.43''N, 126°57'0.64''E）向東北轉移105.72 km；2070s 年代分布質心（46°11'31.04''N,127°1'25.55''E）向西南轉移27.4 km（表3、圖6）。

表3 不同時期紅皮云杉適生區質心坐標和遷徙距離①Tab.3 Centroid coordinates and migration distance of suitable growth areas for P.koraiensis in different periods

圖6 不同時期紅皮云杉適生區質心的遷移路線Fig.6 Migration route of centroid in the suitable growth area of P.koraiensis in different periods

3 討論

3.1 影響紅皮云杉分布的環境因子

根據自然選擇學說，物種形成是由于環境不同而引起的適應性改變。環境是物種形成中最重要的因素，而全球氣候變化在環境變化中占有重要地位。本研究采用19 個環境變量，計算Pearson 相關系數，排除高相關性變量中貢獻率較低的變量，并通過AUC評價建模精度。本研究訓練集和測試集的AUC 均高于0.99，說明模擬結果具有較高可靠性，當前適生分布區的預測更接近實際情況。

氣候是決定物種分布變化的關鍵因素，其中溫度因子通過影響植物光合作用、呼吸作用來控制物種分布；而植物生長發育也受降水條件制約（Abeliet al.，2020）。本研究采用MaxEnt 自帶刀切法，不僅可提高模型預測可信度（De Cauweret al.，2014），而且能客觀、精確地評價19 個氣候變量。最暖季度降水量（Bio18）、氣溫季節性變動系數（Bio4）、降水量季節性變化（Bio15）對紅皮云杉的分布影響最大，其次為年均氣溫（Bio1）、最干季度降水量（Bio17）、晝夜溫差與年溫差比值（Bio3）等環境因子。其中，與水分相關的因子（Bio18、Bio15、Bio17）貢獻率為66.7%，說明與溫度因子相比，紅皮云杉適生區分布受降水影響更大，與王改妮（2018）、段國珍（2019）對紅皮云杉及其近緣種分布格局的模擬結果相似。這是因為紅皮云杉在云杉屬物種分布區的最北端，北方河流冬季有封凍期，降水量較少，對于耐陰、喜濕潤的云杉屬植物來說，生長季降水量不足會限制生長，長期缺水則會導致植物的適應性分化。與其近緣種白杄（Picea meyeri）和沙地云杉（Picea mongolica）相比，紅皮云杉適應濕冷環境的習性十分明顯，無論是大氣濕度還是土壤水分供給狀況均高于前二者。與云杉組（Sect.Picea）內其他類群相比，紅皮云杉是對水分需求較高的種群，其分布區內的年均降水量僅次于云杉（Picea asperata）和青扦（Picea wilsonii）；與同域種魚鱗云杉（Picea jezoensis）（云杉屬魚鱗云杉組sect.Omorica）相比，紅皮云杉更適應低濕的土壤環境，在排水良好的山坡，紅皮云杉的競爭力不及魚鱗云杉，分布數量較少。此外，影響擦木（Sassafras tzumu）、柳杉（Cryptomeria japonica）等植物適生區分布的主要限制因子也是降水（劉想等，2018；戴旻峻等，2022），這些植物穩定生長在山地溝谷和水分肥沃的地區，同樣，紅皮云杉主要生長在河谷溪旁、河岸階地和山麓。綜上表明，紅皮云杉對水分的需求較高，降水是影響紅皮云杉分布和擴散潛力的主要限制因子，因此紅皮云杉在推廣種植或異地保護時，可選擇在降水充沛地區。

溫度在影響紅皮云杉分布上也起到較重要作用。根據刀切法結果，年均氣溫（Bio1）對紅皮云杉分布的影響不可替代，Bio1 的適宜范圍在-1.4～5.6 ℃，最適值為4.3 ℃，在-1.4～4.3 ℃時，分布概率隨年均溫升高而增大，在4.3～5.6 ℃時，分布概率隨年均溫升高而降低，該因子可能限制紅皮云杉繼續向南擴張。氣溫季節性變動系數（Bio4）的貢獻率達28.9%，Bio4 的適宜范圍在129.4～173.5，最適值為155.6。Bio4 代表溫度在一年中的變化幅度，數值越大溫度變幅越大，因此均一穩定的生長環境可能對促進紅皮云杉種群擴張有積極作用。

3.2 紅皮云杉在不同氣候情景下的分布格局變化和質心轉移

植物在與氣候的長期交互作用下形成一定的地理分布格局。通常，擁有更廣生態位的物種具有更強的環境適應能力，可充分利用其周邊的自然資源（Yinget al.，2016），因此在預測紅皮云杉林分更新狀況時必須考慮其對氣候變化的適應情況。本研究采用數據庫數據與實地調查數據相結合的方法，利用古氣候數據、現代氣候數據和未來氣候變化情景，構建紅皮云杉物種分布模型，并重建LIG、LGM 和MH 3個歷史時期的地理分布格局。

研究發現，紅皮云杉種群適宜生境曾經歷多次適應性轉換，包括末次盛冰期的種群擴張以及全新世中期及之后的種群持續收縮。末次間冰期的陸地平均氣溫比現代高1～2 ℃（何韜，2019），與現代更相似；末次盛冰期的陸地平均氣溫比現代低5～10 ℃，而中高緯度地區低20 ℃以上（何韜，2019）。本研究結果和Stewart 等（2010）提出的理論是吻合的，即適應寒冷的植物在間冰期的分布范圍比盛冰期小，如白云杉（Picea glauca）（Stewartet al.，2010 ；Bannisteret al.，2001），以往對云杉屬其他物種的生態位模擬分析中也有一致的結果（Shenet al.，2019；王改妮，2018；段國珍，2019）。相對于末次間冰期，末次盛冰期氣候較冷，紅皮云杉的質心向東南方向轉移；全新世中期氣候再次回暖，紅皮云杉的質心向東北方向轉移。紅皮云杉的質心轉移與氣候變暖條件下物種向溫度更低的區域遷移的規律是相符的。

溫度適當升高可以加速植物生長發育，并影響種群數量變化。隨著全球溫度持續上升，在未來氣候變化2 種溫室氣體排放情景下，紅皮云杉的適生區面積均呈現持續減少趨勢，高度適生區呈片段化分布，全球氣候變暖不利于紅皮云杉適生區擴張。對附子（Aconitum carmichaeli）（韋潔敏，2020） 、 金錢松（Pseudolarix amabilis）（張文秀，2020）、桫欏（Alsophila spinulosa）（楊啟杰等，2021）等植物的潛在分布預測分析也發現，隨著未來全球氣候變暖，這些物種適生區面積會持續減少，與本研究結論一致。未來氣候變化情景中，與RCP4.5 排放路徑相比，RCP8.5 排放路徑下的適生區面積縮減速度較大，紅皮云杉在高排放情景下適生區面積變化更顯著，對氣候變化響應明顯。隨著未來氣候變暖，紅皮云杉總適生區質心先向高緯度遷移，后逐漸向低緯度回落，出現這種情況的主要原因是紅皮云杉潛在適生區在向西南方向海拔更高的長白山地區擴張。全球變暖影響植物的多樣性、分布和遷移模式，大多數物種遷移到高海拔和高緯度地區（Parmesan，2006）。木梨（Pyrus xerophila）、四合木（Tetraena mongolica）未來的潛在分布范圍有向西北方向的高緯度地區遷移的趨勢（劉超等，2020；段義忠等，2019）。未來氣候情景下，紅皮云杉的高適生區縮減，邊緣呈片段化分布，縮減區域包括大興安嶺的塔河、嫩江、蒙克山和高峰，長白山西部的輝南和撫松，張廣才嶺的勃利、樺南、依蘭和海林等地，建議對這些區域的種群加強就地保護、建立自然保護區、減少人為破壞與干擾。小興安嶺腹地和長白山腹地的種群生境較為穩定，是建立紅皮云杉異地保護的最佳候選區域。

3.3 紅皮云杉避難所推測

當全球經歷第四紀冰期的漫長考驗時，眾多生物依靠生物避難所的保護而得以幸存，并在冰期后重新擴散。識別生物避難所有3 個主要途徑：一是基于生物地理模式，通過收集目標生物體的古生物學、生態學和遺傳學數據，量化生物避難所的時空緯度；二是調查物種避難所棲息地的產生過程；三是通過模擬物種分布格局對環境變化的響應來識別避難所（胡忠俊等，2013）。冰期避難所理論表明，避難所在環境波動期間的生態條件相對穩定，可促進物種多樣性的積累。因此物種具有較高水平遺傳多樣性的地理分布區域，可以作為尋找避難所的依據之一（石甜，2021）。

紅皮云杉形成于第三紀中晚期（4.12～7.64 Ma, 百萬年）(Shaoet al.，2019)，說明最先引起紅皮云杉群體遺傳分化的因素是第三紀的地質和氣候事件。在第三紀中晚期，青藏高原出現快速抬升（8～10 Ma），加速了亞洲大陸內部的干燥度和亞洲季風的產生（Anet al.，2001）。第三紀末發生地殼運動，使得大興安嶺沿東側斷層翹起；第四紀初期的構造運動又引起火山噴發，在小興安嶺地區形成大片玄武巖熔巖臺地，同時有輕微翹起，沿北西方向斷裂隆起（樊祺誠等，2015；劉建強等，2017）；第四紀中晚期至全新世，長白山發生多次火山噴發，長白山主體于全新世形成后，地殼相對穩定（王團華，2006）。紅皮云杉的分化時間段發生于第三紀末至第四紀期間，東亞季風氣候增強和地殼變動的時間基本一致。與其分化時間接近的還有第三紀植物楓香（Liquidambar formosana）（10.30 Ma）（孫榮喜，2017）、水青樹（Tetracentron sinense）（9.6 Ma）（Sunet al.，2014）、青岡棟屬（Cyclobalanopsis）（9.07 Ma）（Denget al.，2018）等。

植物在冰期會退縮到避難所，耐寒植物在冰期往往會擴張且在間冰期發生種群收縮，故生物避難所還應包括間冰期避難所（Stewartet al.，2010）。紅皮云杉屬于耐寒植物，通過分子證據、適生分布區模擬結果以及孢粉化石證據可推測紅皮云杉的間冰期避難所。分子證據顯示，在紅皮云杉不同地理種群中，長白山腹地的臨江種群（π= 0.002 35）和天橋嶺種群（π =0.002 32）均具有較高的核苷酸多樣性（nucleotide diversity, π）(Wanget al.，2023)；此外，由于多個歷史時期的物種分布重疊區往往是植物潛在避難所（Carnavalet al.，2008），結合對末次間冰期和末次盛冰期的適生分布區模擬結果推測，長白山腹地可能是紅皮云杉第四紀間冰期避難所，且在該種群進化過程中發揮了重要作用。從生境特征上看，長白山地區屬于紅皮云杉分布區內海拔最高的區域（海拔100 0～274 9 m），較大的海拔落差為紅皮云杉等山地植物適應氣候變化提供了充足的緩沖環境，且該區域末次間冰期時的氣候濕冷（孫廣友等，2012），為紅皮云杉提供了生長所需的充沛水分。此外，小興安嶺腹地也可能是紅皮云杉的避難所之一。分子證據顯示小興安嶺腹地的種群（沾河）核苷酸多樣性較高（π= 0.002 31）（Wanget al.，2023）；此地又是多個時期適生分布區的重疊區域；從孢粉數據來看，松遼盆地發現大量中更新世晚期的云杉屬孢粉化石（石寧，1996），與末次間冰期的時間吻合；小興安嶺處于華力西褶皺帶，第三紀時山體北部與東北方的結雅河低地、西南方的松嫩平原聯成一片，在上新世末更新世初沿北東向和北西向斷裂抬升為山嶺，與結雅河低地及松嫩平原分離（王歡等，2018），其復雜多變的山地峽谷地貌為紅皮云杉提供了豐富的避難場所。因此，結合分子證據、適生區模擬、孢粉化石和地形地貌，推測紅皮云杉在間冰期的生物避難所在小興安嶺腹地和長白山腹地。

由于氣候數據是根據數學方法估計所得，與古代實際氣候存在差異，且未考慮非生物因素的影響，因而MaxEnt 模型在適生分布區預測方面存在一定局限性。然而，通過與其他物種分布模式的比較，MaxEnt具有較高的可靠性和穩定性（Phillipset al.，2006；Elithet al.，2010；朱耿平等，2016）。在模擬物種當代適生分布區時，其結果和物種實際分布范圍具有高度一致性（高文強等，2016；趙光華等，2021；王東升等，2022）。該方法可以呈現物種在不同時期的地理分布格局，判斷環境因素對物種分化和形成的影響，推測物種的生物避難所，預測物種在未來氣候下的分布范圍變化趨勢并為物種保護提供依據（李國慶等，2013；朱耿平等，2013）。

4 結論

MaxEnt 模型預測紅皮云杉適生分布區結果可靠，訓練集和測試集的AUC 均大于0.99。與溫度因子相比，降水對紅皮云杉的地理分布影響更重要，在其推廣種植或者異地保護時，可選擇降水充沛地區。溫度因子可能限制紅皮云杉向南擴張，均一穩定的生長環境對促進其種群擴張有積極作用。由于第三紀末至第四紀期間的地質和氣候事件，紅皮云杉的適生區經歷多次轉換，包括末次盛冰期種群擴張、全新世中期及之后的種群持續收縮。紅皮云杉的適生分布區質心轉移與氣候變暖條件下物種向溫度更低的區域遷移的規律相符。全球氣候變暖不利于紅皮云杉的種群擴張，未來情景下紅皮云杉的適生區呈縮減趨勢。應考慮優先保護大興安嶺北部到中部、長白山西部和張廣才嶺北部等高適生區的邊緣種群。根據分子證據、適生區模擬，結合孢粉化石和地形地貌，推測小興安嶺腹地和長白山腹地是紅皮云杉的間冰期避難所。