王國(guó)崢,耿其芳,肖孟陽(yáng),張孟源,張?jiān)蒲?王中生

南京大學(xué)生命科學(xué)學(xué)院,南京 210023

金錢松(Pseudolarixamabilis(Nelson) Rehd.)是我國(guó)特有單屬種松科植物,也是著名的第三紀(jì)孑遺植物,現(xiàn)主要分布于江蘇南部、安徽南部、福建北部、浙江、江西、湖南、湖北利川至四川萬(wàn)縣交界地區(qū)[1]。最早的化石記錄出現(xiàn)于西伯利亞?wèn)|北部的晚白堊紀(jì)地層中,始新世時(shí)出現(xiàn)于挪威斯匹次卑爾根島西部,美國(guó)、歐洲、亞洲中部以及日本在第三紀(jì)的不同時(shí)期都曾發(fā)現(xiàn)有該屬化石分布[2]。據(jù)化石資料推測(cè),金錢松屬可能起源于白堊紀(jì)中世紀(jì)晚期,第三紀(jì)時(shí)曾發(fā)展為北半球中高緯度廣布屬,其中在日本和歐洲一直維持至更新世第一次冰期降臨[3]。在新生代末期冰期,它從歐洲發(fā)生地向我國(guó)南方遷移,最后定居于長(zhǎng)江流域,形成現(xiàn)今的地理分布格局[4]。

由于金錢松特殊的進(jìn)化及遷移歷史,其在松科植物系統(tǒng)發(fā)育研究中備受關(guān)注,同樣對(duì)古生態(tài)以及古氣候的研究具有重要意義[5]。從其化石記錄可發(fā)現(xiàn)金錢松的分布與地質(zhì)變化過(guò)程密切相關(guān),地質(zhì)變化導(dǎo)致氣候劇烈變化,金錢松分布范圍隨之發(fā)生改變。一些學(xué)者已從金錢松植物學(xué)特征[5]、種群生態(tài)學(xué)[6]、以及遺傳多樣性[7- 8]等方面開(kāi)展了相關(guān)研究,其中高燕會(huì)[7]和劉俊[8]等利用分子標(biāo)記發(fā)現(xiàn)金錢松天然種群具有較高的遺傳多樣性,認(rèn)為遺傳因素不是金錢松的致瀕因子,其呈狹域分布的主要原因可能來(lái)自外界自然因素。

近年來(lái)生態(tài)位模型在預(yù)測(cè)、解析物種分布范圍及其變化等方面應(yīng)用廣泛,利用物種已知的分布數(shù)據(jù)和環(huán)境變量,根據(jù)一定的算法運(yùn)算來(lái)構(gòu)建模型,并將運(yùn)算結(jié)果投射到不同的時(shí)間和空間中以預(yù)測(cè)物種的潛在適生區(qū)[9]。目前較為常用的生態(tài)位模型較多,每種模型都可以獨(dú)立預(yù)測(cè)出潛在適生區(qū),但各模型都存在一定的偏好性[10],本文參考Miguel B. Araújo[10]的方法,采用了4種基于不同算法的模型(GARP、Bioclim、Domain和Maxent)預(yù)測(cè)金錢松的全國(guó)潛在適生區(qū),利用集合預(yù)測(cè)系統(tǒng)思想,綜合多種模型(3種以上)預(yù)測(cè)結(jié)果,盡可能降低由經(jīng)驗(yàn)性選擇一種模型導(dǎo)致的假陰性或假陽(yáng)性影響,同時(shí)一種模型的缺陷可能被另一種模型所彌補(bǔ),從而提高預(yù)測(cè)的科學(xué)性。GARP是基于遺傳算法的規(guī)則組合進(jìn)行局域環(huán)境空間建模[11]；Bioclim基于矩形框架模型,每一種環(huán)境因子被抽象成一個(gè)矩形框,絕大多數(shù)分布點(diǎn)(95%以上點(diǎn))位于這些框內(nèi),多個(gè)矩形框構(gòu)成限制范圍,若某點(diǎn)位于這個(gè)范圍內(nèi),便可以認(rèn)為其適宜該物種生存[12]；Domain模型基于Gower算法,通過(guò)計(jì)算所有已知分布點(diǎn)之間的最大相似度,來(lái)評(píng)估預(yù)測(cè)點(diǎn)對(duì)于環(huán)境變量的適宜程度[13]；Maxent模型通過(guò)物種的分布數(shù)據(jù)和環(huán)境數(shù)據(jù),找出物種分布規(guī)律的最大熵,從而對(duì)物種的分布進(jìn)行預(yù)測(cè)[14- 20]。

1 數(shù)據(jù)獲取和研究方法

1.1 金錢松樣本分布數(shù)據(jù)

金錢松的地理分布位點(diǎn)主要通過(guò)中國(guó)數(shù)字植物標(biāo)本館(http://www.cvh.org.cn/)以及中國(guó)科學(xué)院北京植物所標(biāo)本館(http://pe.ibcas.ac.cn/)查閱,共獲取44個(gè)分布記錄數(shù)據(jù),準(zhǔn)確到縣。用Google Earth 轉(zhuǎn)換出分布點(diǎn)的經(jīng)緯度坐標(biāo),輸入Excel保存成.csv格式。

1.2 環(huán)境變量數(shù)據(jù)及預(yù)處理

環(huán)境變量共22個(gè),包括19個(gè)氣候因子和3個(gè)地形因子,19個(gè)氣候因子來(lái)源于Worldclim (http://www.worldclim.org/),包括末次盛冰期(Last Glacial Maximum,LGM)、當(dāng)前(Current)和2070年的數(shù)據(jù),分辨率2.5 arc-minutes,將環(huán)境因子數(shù)據(jù)用ArcGIS 10.2統(tǒng)一轉(zhuǎn)化為ASCII格式以便使用。

由于環(huán)境變量之間具有多重共線性,會(huì)導(dǎo)致預(yù)測(cè)分布過(guò)度擬合[21],因此本文參考張琴[15]對(duì)環(huán)境共線性診斷的方法,對(duì)環(huán)境因子進(jìn)行Spearman秩相關(guān)分析,篩選出相關(guān)性較低的環(huán)境變量(Spearman系數(shù)<0.75),并在Spearman系數(shù)大于0.75的變量中選出具有生態(tài)學(xué)重要意義的環(huán)境變量,共得到10個(gè)變量(表1)。

表1 環(huán)境變量

1.3 相關(guān)軟件

4種生態(tài)位模型所用軟件：基于GARP模型預(yù)測(cè)的Desktopgarp(Version1.1.6),基于Bioclim和Domain模型預(yù)測(cè)DIVA-GIS(Version 7.5)；基于最大熵模型預(yù)測(cè)軟件Maxent(Version 3.3.3)。

數(shù)據(jù)處理軟件有：Excel用于分布數(shù)據(jù)記錄以及數(shù)據(jù)格式轉(zhuǎn)換,DIVA-GIS(Version 7.5)用于圖層格式轉(zhuǎn)換,ArcGIS(Version 10.2)用于圖層的處理和數(shù)據(jù)的格式轉(zhuǎn)換。

1.4 模型預(yù)測(cè)與檢驗(yàn)方法

1.4.1GARP模型預(yù)測(cè)

將物種已知分布點(diǎn)數(shù)據(jù)轉(zhuǎn)換成Desktopgarp模型軟件支持的格式,通過(guò)模型軟件中 Upload Data Points功能加載,選擇70%的分布數(shù)據(jù)作為訓(xùn)練數(shù)據(jù),其余數(shù)據(jù)作為檢驗(yàn)數(shù)據(jù)[14]。環(huán)境數(shù)據(jù).asc文件通過(guò)Desktopgarp中的Dataset Manager處理轉(zhuǎn)換成Desktopgarp可識(shí)別的格式.raw,以數(shù)據(jù)集的形式加載到軟件中。

1.4.2Maxent模型預(yù)測(cè)

將已知分布點(diǎn)數(shù)據(jù)文件格式在Excel中轉(zhuǎn)為.csv格式文件導(dǎo)出,通過(guò)Browse加載到Maxent軟件中。再將10個(gè)環(huán)境變量ASCⅡ文件的環(huán)境數(shù)據(jù)通過(guò)Browse加載到Maxent軟件中。并進(jìn)一步利用Maxent軟件對(duì)金錢松在末次盛冰期(LGM)氣候和未來(lái)氣候(2070年)條件下進(jìn)行模擬分布[22],再用ArcGIS 10.2 SDMs Tool模擬金錢松從末次盛冰期到當(dāng)前的分布遷移變化[23]。

1.4.3Bioclim和Domain模型分析

在DIVA-GIS中首先添加.shp格式的訓(xùn)練數(shù)據(jù)集,再將10個(gè)環(huán)境變量文件ASCⅡ轉(zhuǎn)換成.grd格式,生成stack數(shù)據(jù)集。在Modeling-Bioclim/Domain模塊中添加stack格式的環(huán)境數(shù)據(jù)集,依次進(jìn)行Bioclim和Domain模型的預(yù)測(cè)[15]。

1.5 模型預(yù)測(cè)結(jié)果檢驗(yàn)

受試者工作特征曲線(Receiver Operating Characteristic,ROC)和Kappa統(tǒng)計(jì)量廣泛用于模型精度的檢驗(yàn)。本文參考張琴[15]的方法,利用DIVA-GIS軟件,將10組不同的訓(xùn)練和測(cè)試數(shù)據(jù)用4種模型預(yù)測(cè)得到的40個(gè)預(yù)測(cè)分布圖生成stack數(shù)據(jù)集,再導(dǎo)入軟件創(chuàng)建的驗(yàn)證點(diǎn).shp文件,最后得出各模型的AUC值和Kappa值。

受試者工作特征曲線下的面積即為AUC值(Area under recriver operating characteristic curve,AUC),AUC的數(shù)值范圍為0.5—1.0,“0.5”說(shuō)明預(yù)測(cè)結(jié)果為隨機(jī)分布,預(yù)測(cè)值越大表示該模型預(yù)測(cè)結(jié)果越精確。一致性檢驗(yàn)統(tǒng)計(jì)量(Kappa)是對(duì)兩種方法測(cè)定結(jié)果的一致部分進(jìn)行檢驗(yàn),取值范圍是[-1,1],“-1”說(shuō)明兩次判斷的結(jié)果完全一致,“1”說(shuō)明兩次判斷的結(jié)果完全不一致,值越大一致性越好,通用的經(jīng)驗(yàn)法則是Kappa大于0.75表示一致性好,小于0.4表示一致性差[24]。

1.6 影響適生分布的環(huán)境因子分析

利用Maxent軟件中的“刀切法”(Do jackknife to measure variable importance)可以顯示出環(huán)各個(gè)境變量對(duì)預(yù)測(cè)金錢松分布的貢獻(xiàn)。其他3種模型(GARP、Bioclim、Domain)預(yù)測(cè)軟件不具備分析環(huán)境因子對(duì)預(yù)測(cè)分布貢獻(xiàn)的功能,因此將GARP、Bioclim、Domain等3種模型基于10個(gè)環(huán)境因子預(yù)測(cè)的適生區(qū)作為對(duì)照組,分別以“刀切法”分析影響較大的三個(gè)環(huán)境因子為單因子,或作為共同因子來(lái)預(yù)測(cè)金錢松的適生區(qū),進(jìn)行預(yù)測(cè)結(jié)果的比較。

2 結(jié)果與分析

2.1 4種模型預(yù)測(cè)結(jié)果

從GARP模型預(yù)測(cè)圖中(圖1)可以看出金錢松在我國(guó)的適生區(qū)集中在江蘇、浙江、安徽、江西、河南、湖北、湖南、貴州、重慶、四川東部、陜西南部以及云南北部地區(qū),集中分布在24.43°—33.35°N和106.41°—123.42°E之間,而搜集的樣本數(shù)據(jù)在27.27°—34.05°N和108.41°—121.80°E之間。

Maxent預(yù)測(cè)分布圖(圖1)顯示金錢松高適區(qū)主要集中在江蘇南部、安徽南部、湖北南部、江西北部以及浙江西北部。綠色代表中度適生區(qū),依次擴(kuò)散開(kāi)來(lái),最大范圍向北達(dá)江蘇省潥陽(yáng)市,向南至湖南衡山磨鏡臺(tái),最西可達(dá)重慶市萬(wàn)州區(qū),最東可達(dá)浙江東天目山。

通過(guò)DIVA-GIS的Modeling功能運(yùn)算得出的Bioclim模型預(yù)測(cè)圖(圖1)和Domain模型預(yù)測(cè)圖(圖1),可觀察到代表最適區(qū)的紅色區(qū)域在這兩種模型中占較小比例,Bioclim預(yù)測(cè)分布范圍較小,顯示金錢松可在浙江、湖南、安徽南部、湖北南部以及江西西北部高度適生。Domain模型預(yù)測(cè)結(jié)果分布范圍較大,總體與GARP模型預(yù)測(cè)分布范圍接近,高適與最適分布范圍與Maxent相似,但在河南新鄉(xiāng)、南陽(yáng)以及云南昆明出現(xiàn)最適分布區(qū)。

2.2 不同模型預(yù)測(cè)準(zhǔn)確性評(píng)估

本文采用10組訓(xùn)練數(shù)據(jù)和測(cè)試數(shù)據(jù)對(duì)4個(gè)模型進(jìn)行ROC曲線分析和Kappa一致性檢驗(yàn)。GARP模型AUC平均值為0.922,Maxent模型AUC平均值為0.980,Bioclim模型AUC平均值為0.957,Domain模型的AUC平均值為0.940(見(jiàn)表2),均遠(yuǎn)遠(yuǎn)超過(guò)隨機(jī)模型(AUC=0.5)。并且4種模型的Kappa均值都大于0.75,表明4種模型預(yù)測(cè)的一致性均較顯著,預(yù)測(cè)精度較高。4種模型AUC值的標(biāo)準(zhǔn)偏差排序?yàn)锽ioclim > Domain > Maxent >GARP。

表2 4種模型的AUC值和Kappa值

2.3 影響適生區(qū)的環(huán)境因子分析

各個(gè)境變量對(duì)預(yù)測(cè)金錢松分布的貢獻(xiàn)如圖2,在10個(gè)環(huán)境因子中,年均溫(Bio1)對(duì)于金錢松分布的貢獻(xiàn)最大,年降水量(Bio12)、最濕月降水量(Bio13)和最熱月極高溫度(Bio5)對(duì)金錢松分布的影響也比較大,溫度日較差(Bio2)和海拔對(duì)金錢松分布的影響比較小。

圖1 4種生態(tài)位模型的金錢松適生區(qū)預(yù)測(cè)圖Fig.1 Potential distribution of P.amabilis based on four ecological niche models

圖2 刀切法檢測(cè)生物氣候變量對(duì)分布的重要性Fig.2 Import of environmental variables for prediction based on jackknife test

從單因子預(yù)測(cè)結(jié)果來(lái)看(圖3),3種模型基于年均溫(Bio1)的模擬結(jié)果與各自預(yù)測(cè)適生區(qū)最為接近,但范圍較適生區(qū)均偏大。橫向?qū)Ρ?種模型分別基于Bio1、Bio12、Bio13等3種單因子的模擬結(jié)果,基于年均溫(Bio1)的預(yù)測(cè)結(jié)果相對(duì)較為接近,說(shuō)明年均溫(Bio1)對(duì)這3種模型的預(yù)測(cè)結(jié)果也都有顯著的影響?；谀杲邓?Bio12)與最濕月降水量(Bio13)兩種因子的預(yù)測(cè)分布圖都與適生區(qū)差異較大,且三種模型之間對(duì)比預(yù)測(cè)分布圖差異也較大,說(shuō)明這兩種因子不是最主要的預(yù)測(cè)分布影響因子。將Bio1、Bio12、Bio13等作為共同影響因子,3種模型的預(yù)測(cè)分布圖與其對(duì)應(yīng)的適生區(qū)都十分接近,說(shuō)明這3種因子是預(yù)測(cè)分布結(jié)果的確存在顯著影響,尤其是Bioclim模型,二者基本完全重疊。此外,基于年均溫(Bio1)的預(yù)測(cè)分布與3因子聯(lián)合預(yù)測(cè)分布最為接近,進(jìn)一步證明年均溫(Bio1)對(duì)金錢松預(yù)測(cè)分布起“框架”作用,而其余環(huán)境因子對(duì)于預(yù)測(cè)分布起到“修飾”作用。

圖3 三種環(huán)境因子模擬結(jié)果Fig.3 The results of the other three models圖a1—a5分別為GARP模型基于年均溫(Bio1)、年降水量(Bio12)、最濕月降水量(Bio13)、三因子(Bio1&12&13)以及對(duì)照組結(jié)果,圖b1—b5分別為Domain模型基于年均溫(Bio1)、年降水量(Bio12)、最濕月降水量(Bio13)、三因子(Bio1&12&13)以及對(duì)照組結(jié)果,圖c1—c5分別為Bioclim模型基于年均溫(Bio1)、年降水量(Bio12)、最濕月降水量(Bio13)、三因子(Bio1&12&13)以及對(duì)照組結(jié)果

2.4 金錢松不同時(shí)期分布格局變化模擬

金錢松從末次盛冰期到當(dāng)前,分布范圍主要向南遷移,分布擴(kuò)張面積達(dá)到178510.08 km2,北部分布范圍有部分收縮,達(dá)到122865.39 km2,凈分布范圍面積擴(kuò)大了55644.69 km2。從當(dāng)前到2070年分布范圍呈現(xiàn)收縮趨勢(shì),主要是南部適生區(qū)大面積收縮,達(dá)到246396.07 km2,占現(xiàn)生分布范圍的14.45%,而向北擴(kuò)張68311.77 km2,僅為收縮面積的27.72%,分布范圍收縮了178084.30 km2(圖4)。

圖4 基于Maxent和ArcGIS金錢松不同時(shí)期的分布格局變化Fig.4 Comparison of the distribution changes of P.amabilis based on ArcGIS and Maxent

3 討論

3.1 金錢松適生范圍預(yù)測(cè)

近年來(lái)物種分布模型發(fā)展迅速,多個(gè)可用的方法模型相繼被開(kāi)發(fā)[21]。本文通過(guò)模型檢驗(yàn)得出四個(gè)模型都可以較為準(zhǔn)確地預(yù)測(cè)出金錢松的潛在適生區(qū)(AUC>0.9),但GARP預(yù)測(cè)結(jié)果的AUC均值最小,精度較差,這可能是金錢松樣本數(shù)據(jù)較少導(dǎo)致。楊會(huì)楓[19]用不同樣本容量檢驗(yàn)GARP預(yù)測(cè)結(jié)果精度,發(fā)現(xiàn)受到樣本容量干擾,只有當(dāng)樣本容量超過(guò)一定閾值,AUC均值才會(huì)趨于穩(wěn)定(大于0.9)。Domain和Bioclim預(yù)測(cè)結(jié)果相比其他模型適生區(qū)范圍較小,這可能與模型受樣點(diǎn)信息有關(guān),金錢松分布數(shù)據(jù)主要來(lái)自于標(biāo)本館,樣本數(shù)據(jù)存在一定偏好性,科研工作者一般根據(jù)自己研究所需進(jìn)行采樣或者單純?yōu)闃?biāo)本館收集標(biāo)本,樣本信息較為零散,缺乏系統(tǒng)性和代表性；其次金錢松樣本數(shù)據(jù)較少可能導(dǎo)致生態(tài)位空間縮小的現(xiàn)象[26],并且由于某些地區(qū)金錢松種群遺傳多樣性較為豐富[5- 6],其物種在不同區(qū)域的生態(tài)位也可能發(fā)生漂移[26]。

Maxent模型預(yù)測(cè)結(jié)果比其他模型更為精確(AUC平均值為0.980),相關(guān)文獻(xiàn)也證明Maxent在大、小樣本數(shù)據(jù)條件下均能很好地預(yù)測(cè)物種適生區(qū)[14- 20]。但不同模型預(yù)測(cè)結(jié)果可以相互補(bǔ)充[21,27- 29],如在Maxent預(yù)測(cè)分布圖上,河南鄭州、信陽(yáng)大別山以及云南昆明均屬于金錢松的低適生區(qū),而Domain模型中這一區(qū)域則表現(xiàn)為高適區(qū),經(jīng)查閱文獻(xiàn)發(fā)現(xiàn)河南鄭州、信陽(yáng)大別山及云南昆明等地區(qū)均有人工成功引種金錢松的記錄[28- 29],可確認(rèn)是金錢松適生區(qū)。

3.2 金錢松分布范圍變化

金錢松為高大落葉喬木,為深根系樹(shù)種[5],其一年生幼苗生長(zhǎng)旺盛期為4—5月和8—9月,這兩個(gè)階段的月均溫和月降水量接近,為金錢松幼苗生長(zhǎng)的最適條件[29]。本文發(fā)現(xiàn)年均溫(Bio1)與年降水量(Bio12)是金錢松潛在適生區(qū)預(yù)測(cè)的主要影響因子,也決定了金錢松分布格局的形成以及遷移方向。末次冰期氣候條件惡劣[23,30],而在金錢松分布區(qū)的浙江天目山、大明山、安徽黃山以及湖南、江西、湖北三省交界處的羅霄山脈是其良好的避難所[31],這些地區(qū)山峰地形多變,小生境水熱條件優(yōu)越,有效的降低了冰期氣候?qū)χ参锓植嫉挠绊憽哪┐伍g冰期(LGM)到當(dāng)前氣候條件下金錢松預(yù)測(cè)分布范圍變化可知,因冰期之后氣溫回暖,金錢松的分布格局逐漸向溫?zé)岬哪戏綌U(kuò)張。從當(dāng)前到2070年金錢松的預(yù)測(cè)分布范圍變化可知,未來(lái)其分布會(huì)一定程度向北擴(kuò)張,但在南方卻大面積收縮,總縮小面積達(dá)到178084.30km2。隨著全球氣候變化加劇[32],溫室效應(yīng)將大大影響植物分布格局以及其物候期[33-35],而植物分布變化是個(gè)較為緩慢的過(guò)程,全球氣溫上升速度遠(yuǎn)遠(yuǎn)超過(guò)金錢松自身遷移的能力,因此必須適時(shí)考慮人為輔助金錢松的北向遷移[33]。

3.3 金錢松潛在適生區(qū)與金錢松保護(hù)

金錢松作為瀕危物種,自然分布范圍有限,且生境受到各種威脅,建立金錢松自然保護(hù)區(qū)與種子園可以有效保護(hù)金錢松[5],自然保護(hù)區(qū)可以最大程度防止人為破壞,種子園可以有效保障種質(zhì)資源和遺傳水平,而二者的建立與物種適宜生境密切相關(guān)[36- 38]。本文結(jié)合多種生態(tài)位模型對(duì)金錢松的潛在適生區(qū)進(jìn)行預(yù)測(cè)及其分布格局變化分析,對(duì)金錢松保護(hù)提出以下建議：除現(xiàn)有浙江西天目山金錢松自然保護(hù)區(qū)外,還可在江西銅鼓縣、湖南張家界和衡陽(yáng)這幾處當(dāng)前金錢松高適分布區(qū)域內(nèi)建立自然保護(hù)區(qū)與種子園,保護(hù)與育苗、培種相結(jié)合,最大程度保護(hù)現(xiàn)有種群；在目前金錢松中等適宜區(qū)以及未來(lái)氣候條件下高適分布區(qū),如安徽北部、河南南部、湖北東南部等地區(qū)可以考慮人工引種,進(jìn)行原土栽植協(xié)助金錢松北移實(shí)驗(yàn),未來(lái)可將其作為園林樹(shù)種,發(fā)展為園林綠植,而4種模型都預(yù)測(cè)為低適區(qū)或不適區(qū)的地方則不適于選擇種植金錢松。