基于GARP模型對(duì)松材線蟲在重慶適生區(qū)的預(yù)測(cè)

李宏群， 蒲永蘭， 劉曉莉， 劉曉梅， 楊清鈺， 丁世敏

(1.長(zhǎng)江師范學(xué)院生命科學(xué)與技術(shù)學(xué)院，重慶 408100； 2.重慶市林業(yè)科學(xué)研究院森林資源保護(hù)與林業(yè)信息技術(shù)研究所，重慶 400036；3.長(zhǎng)江師范學(xué)院圖書館，重慶 408100； 4.長(zhǎng)江師范學(xué)院武陵山片區(qū)綠色發(fā)展協(xié)同創(chuàng)新中心，重慶 408100)

松材線蟲病是由外來(lái)入侵種松材線蟲(Bursaphelenchusxylophilus)引起的，目前被列為我國(guó)第一大林業(yè)外來(lái)有害生物，被稱作松樹“癌癥”“松樹艾滋病”“松樹禽流感”，已給發(fā)生地區(qū)造成巨大的經(jīng)濟(jì)、生態(tài)和社會(huì)損失。我國(guó)自1982年在江蘇省南京市中山陵的黑松上首次發(fā)現(xiàn)該病以來(lái)，短短30多年，疫情已經(jīng)迅速擴(kuò)展到我國(guó)14個(gè)省(市、區(qū))，192個(gè)縣(市、區(qū))，累計(jì)致死松樹5億多株，毀滅松林33萬(wàn)hm2，嚴(yán)重影響了疫區(qū)經(jīng)濟(jì)和社會(huì)的可持續(xù)發(fā)展[1-2]。目前松材線蟲病對(duì)廬山、黃山和三峽庫(kù)區(qū)等的生態(tài)安全構(gòu)成了嚴(yán)重威脅，且到目前為止沒有十分有效的防治方法，形勢(shì)十分嚴(yán)峻。

物種潛在分布區(qū)預(yù)測(cè)模型在瀕危物種保護(hù)、生物多樣性保護(hù)區(qū)設(shè)計(jì)的優(yōu)先性評(píng)估、外來(lái)入侵物種擴(kuò)散潛能以及全球氣候變化對(duì)物種分布區(qū)的影響等研究中得到越來(lái)越廣泛的應(yīng)用[3]。松材線蟲的入侵很難被控制，消耗大量的人力、物力和財(cái)力，以致預(yù)防是比防治更經(jīng)濟(jì)的行為。在過去，國(guó)內(nèi)外對(duì)松材線蟲治病機(jī)制、發(fā)生現(xiàn)狀以及防治措施的研究較多[4-6]，而對(duì)于其潛在生境預(yù)測(cè)的研究相對(duì)較少[1，7-8]。三峽庫(kù)區(qū)重慶段還沒有人研究，且目前重慶的萬(wàn)州區(qū)、涪陵區(qū)、長(zhǎng)壽區(qū)、巴南區(qū)、云陽(yáng)縣和忠縣已經(jīng)被確定是松材線蟲病的高發(fā)疫區(qū)，已表現(xiàn)出擴(kuò)大蔓延的趨勢(shì)，引起了我國(guó)學(xué)者的重視。用來(lái)預(yù)測(cè)物種潛在分布區(qū)的預(yù)設(shè)預(yù)測(cè)規(guī)則的遺傳算法(genetic algorithm for rule-set production，簡(jiǎn)稱GARP)生態(tài)位模型最初由Stockwell創(chuàng)建，是預(yù)設(shè)規(guī)則的遺傳算法，利用物種的已知分布數(shù)據(jù)和環(huán)境數(shù)據(jù)模擬物種的基本生態(tài)位需求，探索物種已知分布區(qū)的環(huán)境特征與研究區(qū)域的非隨機(jī)關(guān)系[9-11]。GARP模型應(yīng)用時(shí)間較早，是目前最準(zhǔn)確的物種潛在分布區(qū)預(yù)測(cè)模型之一，應(yīng)用較為廣泛[9-11]。因此，本研究以松材線蟲在重慶地區(qū)的分布數(shù)據(jù)為基礎(chǔ)，運(yùn)用GARP生態(tài)模型軟件，預(yù)測(cè)其潛在入侵區(qū)域，為制定合理的防治措施提供科學(xué)依據(jù)。

1 材料與方法

1.1 材料

1.1.1 松材線蟲病分布點(diǎn)數(shù)據(jù)來(lái)源 依據(jù)重慶市涪陵區(qū)林業(yè)局病蟲害防治站提供的松材線蟲發(fā)生地分布數(shù)據(jù)，采用現(xiàn)場(chǎng)觀測(cè)，手持GPS定位的方式，提取發(fā)生地的經(jīng)緯度坐標(biāo)，核對(duì)位置后得到63個(gè)有效分布點(diǎn)，根據(jù)GARP軟件要求，將分布點(diǎn)經(jīng)度和緯度儲(chǔ)存成.csv格式的文件。

1.1.2 軟件工具 GARP軟件：從http：//www.nhm.ku.edu/desktopgarp/免費(fèi)下載，版本為1.1.6；GIS軟件由中國(guó)科學(xué)院重慶綠色智能技術(shù)研究院生態(tài)過程與重建研究中心提供，版本為ArcGIS 9.3。

1.1.3 基礎(chǔ)地理數(shù)據(jù) 從國(guó)家基礎(chǔ)地理信息系統(tǒng)(http：Mnfgis.nsdi.gov.cn/)下載獲得1 ∶400萬(wàn)的我國(guó)國(guó)界和省界以及縣界行政區(qū)劃圖，作為分析的底圖。

1.1.4 環(huán)境數(shù)據(jù) 本研究選取影響松材線蟲病分布的19個(gè)生物氣候變量作為主要環(huán)境數(shù)據(jù)(表1)。數(shù)據(jù)來(lái)自Berkeley大學(xué)Worldclim下載中心國(guó)際網(wǎng)站(http.www.worldclim.org)免費(fèi)下載空間，分辨率為30″當(dāng)前狀態(tài)(current conditions 1950—2000)19個(gè)生物氣候變量，所有層數(shù)據(jù)的投影格式均轉(zhuǎn)換為GCS-WGS-1984，然后通過重慶市的地圖對(duì)所獲得數(shù)據(jù)進(jìn)行掩膜獲得三峽庫(kù)區(qū)重慶段數(shù)據(jù)，此過程在ArcGIS 9.3軟件中進(jìn)行，并把所有環(huán)境數(shù)據(jù)轉(zhuǎn)換為ASCⅡ格式。

1.2 模型模擬

在使用該GARP模型時(shí)，每次運(yùn)算隨機(jī)選取50%的數(shù)據(jù)用于建立模型，剩下50%的數(shù)據(jù)用于驗(yàn)證模型誤差。創(chuàng)建模型規(guī)則時(shí)反復(fù)運(yùn)算1 000次或其精度達(dá)到0.01而終止。在選取最優(yōu)的環(huán)境因子時(shí)，首先將19個(gè)環(huán)境因子一起進(jìn)行 1 000 次模型運(yùn)算，將所獲得的溢出誤差95%置信區(qū)間的算術(shù)平均值作為基礎(chǔ)溢出誤差。然后使用刀切法(jackknife)依次省略每一個(gè)環(huán)境因子，創(chuàng)建19組不同的環(huán)境因子組。分析每一組環(huán)境因子與溢出誤差(omission error)之間的關(guān)系，如果一個(gè)環(huán)境因子的缺失導(dǎo)致溢出誤差顯著降低(低于基礎(chǔ)溢出誤差×95%)，則該環(huán)境因子將不包含在最終的分析中。反之，如果一個(gè)環(huán)境因子的缺失導(dǎo)致溢出誤差高于基礎(chǔ)溢出誤差，則認(rèn)為這一因子對(duì)預(yù)測(cè)結(jié)果影響顯著[12]。

在最終分析中，重復(fù)創(chuàng)建1 000次模型，創(chuàng)建模型規(guī)則時(shí)重復(fù)運(yùn)算次數(shù)最大值(max iteration)設(shè)為5 000，收斂極限(convergence limit)為0.01，規(guī)則類型(rule types)使用默認(rèn)設(shè)置，并生成arc格式的圖層。運(yùn)算完畢后，選取所有內(nèi)部遺漏誤差和外部遺漏誤差為0的預(yù)測(cè)中溢出誤差最小的50個(gè)模型，即為最優(yōu)模型。利用ArcGIS的空間擴(kuò)展模塊(spatial analyst tools)，將最優(yōu)模型的預(yù)測(cè)結(jié)果進(jìn)行等權(quán)重疊加。疊加的結(jié)果為每一個(gè)柵格的數(shù)值為0～50，0表示50個(gè)預(yù)測(cè)模型都預(yù)測(cè)該柵格不適合松材線蟲繁殖，而50表示50個(gè)預(yù)測(cè)模型都預(yù)測(cè)該柵格適合松材線蟲繁殖。同時(shí)，把每一個(gè)柵格的值除以50來(lái)表示該區(qū)域的適宜性，稱之為重疊系數(shù)(overlap index，簡(jiǎn)稱OI)。重疊系數(shù)值越大，表示該柵格的環(huán)境參數(shù)越適合物種的生存。因此，根據(jù)重疊系數(shù)的大小可以將預(yù)測(cè)的區(qū)域劃分為：最佳適生區(qū)(0.9

2 結(jié)果與分析

2.1 環(huán)境因子選擇及模型選擇

19個(gè)環(huán)境因子共同進(jìn)行模型運(yùn)算所獲得的基礎(chǔ)溢出誤差為28.530 7。選取最優(yōu)模型時(shí)，依次去除19個(gè)環(huán)境因子所獲得溢出誤差的平均值如表1所示。其中溫度季節(jié)性變化、最熱月最高溫度、最濕潤(rùn)季度平均溫度、最冷季度平均溫度、最干旱月降水量、降水量季節(jié)變化和最熱季度降水量被去除時(shí)，溢出誤差明顯下降(<28.530 7×0.95=27.104 2)，而去除年平均溫度、月均值晝夜溫差、等溫性、最冷月最低溫度、氣溫年變化范圍、最干旱季度平均溫度、最熱季度平均溫度、年降水量、最濕潤(rùn)月降水量、最濕潤(rùn)季度降水量、最干旱季度降水量和最冷季度降水量則會(huì)使溢出誤差明顯升高。因此，溫度季節(jié)性變化、最熱月最高溫度、最濕潤(rùn)季度平均溫度、最冷季度平均溫度、最干旱月降水量、降水量季節(jié)變化和最熱季度降水量在最終的分析中被剔除，使用其他12個(gè)環(huán)境因子建立 1 000個(gè)模型，并選取所有內(nèi)部遺漏誤差和外部遺漏誤差為0的預(yù)測(cè)中72個(gè)模型，即為最優(yōu)模型。以獨(dú)立外部檢驗(yàn)數(shù)據(jù)集為基礎(chǔ)的卡方檢驗(yàn)顯示這72個(gè)模型在統(tǒng)計(jì)上都高度顯著(P<0.01)，即模型的預(yù)測(cè)能力遠(yuǎn)遠(yuǎn)好于隨機(jī)模型。

表1 生態(tài)位模型參數(shù)

注：“+”、“-”分別表示該環(huán)境因子包含、不包含在最終的分析中。

2.2 潛在分布區(qū)的預(yù)測(cè)

把這72個(gè)模型作為最優(yōu)的模型，利用ArcGIS的空間擴(kuò)展模塊，將從最優(yōu)模型中隨機(jī)選取的50個(gè)預(yù)測(cè)結(jié)果進(jìn)行等權(quán)重疊加。模型預(yù)測(cè)圖與我國(guó)重慶行政區(qū)圖疊加的結(jié)果如圖1所示，紅、橙、綠、白4種顏色依次表示潛在分布的可能性由高到低。由圖1可看出，松材線蟲病適生區(qū)主要集中在重慶中部，即三峽庫(kù)區(qū)腹地，然后向重慶東北、東南以及南方向傳播。具體分布為：最佳適生區(qū)有涪陵區(qū)、長(zhǎng)壽區(qū)、豐都縣、墊江縣、忠縣、梁平區(qū)和萬(wàn)州區(qū)；高度適生區(qū)有南川區(qū)、武隆區(qū)、萬(wàn)盛區(qū)、綦江區(qū)、南岸區(qū)、巴南區(qū)、江北區(qū)、渝北區(qū)、彭水縣、石柱縣、奉節(jié)縣、開縣和云陽(yáng)縣；中度適生區(qū)有巫山縣、巫溪縣、酉陽(yáng)縣、黔江區(qū)、江津區(qū)、永川區(qū)、永昌縣；非適生區(qū)有沙坪壩區(qū)、北碚區(qū)、合川區(qū)、九龍坡區(qū)、大渡口區(qū)、大足縣、壁山縣、銅梁縣、潼南縣、秀山縣等。

3 討論

3.1 環(huán)境因子分析

在篩選用于入侵生物的最優(yōu)環(huán)境因子集合時(shí)，通常會(huì)采用一系列的刀切法分析，通過多次重復(fù)計(jì)算，每次省略1個(gè)環(huán)境參數(shù)，然后分析某一環(huán)境參數(shù)與遺漏誤差之間的相關(guān)性，確定最優(yōu)的環(huán)境參數(shù)，即如果1個(gè)環(huán)境參數(shù)的存在與遺漏誤差呈正相關(guān)關(guān)系，表明該環(huán)境參數(shù)會(huì)降低模型預(yù)測(cè)的準(zhǔn)確性，該環(huán)境參數(shù)會(huì)在最終的分析中被去掉。此方法有普遍適用性[10-12]。對(duì)于GARP模型來(lái)說，雖然可以產(chǎn)生成功率較高且具有顯著統(tǒng)計(jì)意義的結(jié)果，但預(yù)測(cè)產(chǎn)生了破碎化的過預(yù)測(cè)生境[13]。有研究表明，GARP可以成功地預(yù)測(cè)到物種大部分的分布點(diǎn)，同時(shí)也會(huì)預(yù)測(cè)到遠(yuǎn)遠(yuǎn)大于物種已知分布區(qū)的分布范圍[14-15]。為使預(yù)測(cè)范圍更加精確，以便實(shí)施重點(diǎn)地段的監(jiān)測(cè)和預(yù)警，本研究在對(duì)環(huán)境因子進(jìn)行篩選時(shí)，并未使用遺漏誤差，而是使用溢出誤差[12]。在Desktop Garp模型運(yùn)算的結(jié)果列表中，溢出誤差代表的是預(yù)測(cè)結(jié)果超出實(shí)際分布區(qū)所占的百分比。本研究在篩選最優(yōu)環(huán)境因子集合時(shí)，根據(jù)刀切法分析，溫度季節(jié)變化、最熱月最高溫度、最濕潤(rùn)季度平均溫度、最冷季度平均溫度、最干旱月降水量、降水量季節(jié)變化和最熱季度降水量被去除時(shí)溢出誤差明顯下降，表明這7個(gè)環(huán)境參數(shù)與溢出誤差呈正相關(guān)，說明這7個(gè)環(huán)境參數(shù)會(huì)降低模型預(yù)測(cè)的準(zhǔn)確性。因此，這7個(gè)環(huán)境參數(shù)在最終的分析中被去掉。其他的12個(gè)環(huán)境參數(shù)被用來(lái)建立最佳模型，也說明這12個(gè)環(huán)境參數(shù)對(duì)其分布影響顯著，分別為年平均溫度、月均值晝夜溫差、等溫性、最冷月最低溫度、氣溫年變化范圍、最干旱季度平均溫度、最熱季度平均溫度、年降水量、最濕潤(rùn)月降水量、最濕潤(rùn)季度降水量、最干旱季度降水量和最冷季度降水量。Graham在研究大量化石生物的基礎(chǔ)上認(rèn)為，物種的生境更多的與溫度和降水的極值有關(guān)，而不是與平均值有關(guān)[16]，以及外來(lái)物種入侵種往往具有較強(qiáng)的抗逆性[10-11]。本研究結(jié)果正好驗(yàn)證了上述研究結(jié)果，也符合三峽大壩修建后極端氣候在重慶反復(fù)出現(xiàn)的事實(shí)。另外，最熱季度平均溫度、降水季節(jié)性變化以及最溫暖月最高溫度是影響松材線蟲傳媒昆蟲——松墨天牛(Monochamusalternatus)的主要環(huán)境因子[1]，因此，在研究松材線蟲的潛在生境時(shí)，可以綜合考慮寄主植被和傳媒昆蟲的生態(tài)環(huán)境因子。

3.2 潛在分布區(qū)分析

用GARP生態(tài)位模型模擬外來(lái)入侵物種在入侵地區(qū)的潛在分布是在展開較深入的生物學(xué)研究之前探知入侵物種可能分布的有效手段，它不需要物種的生物學(xué)資料，因此對(duì)具有潛在入侵威脅或初侵入的檢疫性物種的風(fēng)險(xiǎn)分析和預(yù)防控制有重要的參考意義[11]。松材線蟲生態(tài)位模型的預(yù)測(cè)結(jié)果主要集中在重慶中部，即三峽庫(kù)區(qū)腹地，然后向重慶東北、東南以及南方向傳播，與其在重慶已入侵地區(qū)的地理生態(tài)式樣基本相符，所有已經(jīng)入侵的地區(qū)幾乎均被生態(tài)位模型預(yù)測(cè)為最佳適生區(qū)域和高度適生區(qū)域，其中涪陵區(qū)、長(zhǎng)壽區(qū)、豐都縣、墊江縣、忠縣、梁平區(qū)和萬(wàn)州區(qū)等均處在最佳適生區(qū)內(nèi)，即三峽庫(kù)區(qū)腹地，并與高度適生區(qū)緊密相連，即南川區(qū)、武隆區(qū)、萬(wàn)盛區(qū)、綦江區(qū)、南岸區(qū)、巴南區(qū)、江北區(qū)、渝北區(qū)、彭水縣、石柱縣、奉節(jié)縣、開縣和云陽(yáng)縣，在這些高度適生區(qū)中也發(fā)現(xiàn)零星松材線蟲病適生區(qū)，說明這些高度適生區(qū)就是未來(lái)松材線蟲病可能爆發(fā)的區(qū)域。因此在最佳適生區(qū)和高度適生區(qū)域應(yīng)該進(jìn)行嚴(yán)格的檢疫和防控措施。另外，控制松材線蟲病的傳媒昆蟲——松墨天牛，因?yàn)榭刂扑赡炫Ｊ乔袛嗨刹木€蟲自然傳播途徑的有效方法。

[1]巨云為，李明陽(yáng)，吳文浩. 江蘇省松材線蟲發(fā)生的預(yù)測(cè)方法[J]. 林業(yè)科學(xué)，2010，46(12)：91-96.

[2]Robinet C，Roques A，Pan H，et al. Role of human-mediated dispersal in the spread of the pinewood nematode in China[J]. PLoS One，2009，4(2)：e4646.

[3]Guisan A，Thuiller W. Predicting species distribution：offering more than simple habitat models[J]. Ecology Letters，2005，8(9)：993-1009.

[4]徐福元，鄭華英，劉云鵬，等. 馬尾松種源對(duì)松褐天牛成蟲取食、松材線蟲病的抗性分析[J]. 林業(yè)科學(xué)，2011，47(1)：101-106.

[5]Kanetani S，Kikuchi T，Akiba M，et al. Detection ofBursaphelenchusxylophilusfrom old discs of deadPinusarmandiivar.amamianatrees using a new detection kit[J]. Forest Pathology，2011，41(5)：387-391.

[6]何龍喜，吳小芹，吉 靜，等. 不同松樹感染松材線蟲后NO和核酸酶變化與其抗病性的關(guān)系[J]. 林業(yè)科學(xué)，2012，48(1)：109-114.

[7]吳文浩，李明陽(yáng). 基于生態(tài)位模型的松材線蟲潛在生境預(yù)測(cè)方法研究[J]. 林業(yè)調(diào)查規(guī)劃，2009，34(5)：33-38.

[8]Robinet C，Opstal N V，Baker R，et al. Applying a spread model to identify the entry points from which the pine wood nematode，the vector of pine wilt disease，would spread most rapidly across Europe[J]. Biological Invasions，2011，13(12)：2981-2995.

[9]Stockwell D，Peters D. The GARP modelling system：problems and solutions to automated spatial prediction[J]. International Journal of Geographic Information Systems，1999，13(2)：143-158.

[10]余 巖，陳立立，何興金. 基于GARP的加拿大一枝黃花在中國(guó)的分布區(qū)預(yù)測(cè)[J]. 云南植物研究，2009，31(1)：57-62.

[11]周先容，余 巖，周頌東，等. 巴山榧樹地理分布格局及潛在分布區(qū)[J]. 林業(yè)科學(xué)，2012，48(2)：1-8.

[12]王 瑞. 我國(guó)嚴(yán)重威脅性外來(lái)入侵植物入侵與擴(kuò)散歷史過程重建及其潛在分布區(qū)的預(yù)測(cè)[D]. 北京：中國(guó)科學(xué)院研究生院，2006.

[13]馬松梅，張明理，張宏祥，等. 利用最大熵模型和規(guī)則集遺傳算法模型預(yù)測(cè)孑遺植物裸果木的潛在地理分布及格局[J]. 植物生態(tài)學(xué)報(bào)，2010，34(11)：1327-1335.

[14]Hernandez P A，Graham C H，Master L L，et al. The effect of sample size and species characteristics on performance of different species distribution modeling methods[J]. Ecography，2006，29(5)：773-785.

[15]Pearson R G，Raxworthy C J，Nakamura M，et al. Predicting species distributions from small numbers ofoccurrence records：a test case using cryptic geckos in Madagascar[J]. Journal of Biogeography，2007，34(1)：102-117.

[16]Graham A. The current status of the legume fossil record in the Caribbean region[M]. Lodon：the Royal Botanic Gardens Press，1992：161-167.