基于廣義可加模型分析溫濕度因子對斑翅果蠅田間種群動態(tài)的影響

陳 曉，趙遠鵬，李 欣，肖 春，李正躍，張立敏**

(1.云南農(nóng)業(yè)大學植物保護學院，昆明 650201；2. 省部共建云南生物資源保護與利用國家重點實驗室，昆明 650201；3. 昆明英茂花卉產(chǎn)業(yè)有限公司，昆明 650000；4. 云南農(nóng)業(yè)大學大數(shù)據(jù)(信息工程)學院，昆明 650201)

斑翅果蠅DrosophilasuzukiiMatsumura，又稱櫻桃果蠅，隸屬雙翅目Diptera環(huán)裂亞目Cylorrhapha果蠅科Drsophilidae果蠅屬Drosophila水果果蠅亞屬Sophophora(孫鵬等，2011)。斑翅果蠅寄主范圍非常廣泛，已知寄主涉及100多種，主要包括藍莓、櫻桃、草莓、葡萄等(Kenisetal.，2016)，且該果蠅的主要寄主植物為高密度單一種植模式，并且該果蠅能以其它野生寄主和觀賞植物作為庇護場所和再侵染源(Klicketal.，2012)，所以容易導致該果蠅種群數(shù)量在種植區(qū)域內(nèi)快速增長甚至爆發(fā)。近年來，隨著各國之間的貿(mào)易往來加深，在世界范圍內(nèi)，斑翅果蠅除了在適宜的自然或管理環(huán)境中擴散，也隨著人為的商品運輸及交通工具的攜帶進行傳播，還有以具潛在風險的有害生物為載體進行協(xié)同擴散，利用廣泛分布于邊境地帶的貿(mào)易商品進行擴散，更有甚者附帶于農(nóng)副產(chǎn)品和廢棄物中進行快速地分布擴散(Hauseretal.，2009)?，F(xiàn)已有多個國家相關檢疫部門高度重視斑翅果蠅日漸增長的危害，發(fā)布了相關預警或將其列入檢疫性有害生物名單(孫鵬等，2011)。

光照、溫度、濕度、降雨等氣象因子可綜合作用于昆蟲的生長發(fā)育、生殖繁衍、行為活動、生存分布等，其中溫度是影響上述生物學指標的關鍵氣象因子，自然環(huán)境中溫度的變化直接影響昆蟲新陳代謝、生長以及繁殖等生命活動的速率，從而對昆蟲的發(fā)生數(shù)量、體形大小、行為變異以及活動期長短產(chǎn)生十分重要的直接作用(石保坤，2012)；環(huán)境濕度(包括空氣相對濕度和土壤含水量)的變化可直接導致昆蟲體內(nèi)含水量的變化，影響蟲體的水分平衡、生長發(fā)育及繁殖速率、抗逆性、行為誘導信號等，進而對其個體發(fā)育及群體發(fā)生產(chǎn)生影響(常曉娜等，2008)。因此，外界環(huán)境溫度、濕度的變化是昆蟲種群隨著季節(jié)消長變化的主要因素，且溫度與濕度因子聯(lián)合作用下可增強上述現(xiàn)象的表達，在高溫和低相對濕度的雙重影響下，不育或可育斑翅果蠅的生存和繁殖均受到負面影響(Krügeretal.，2021)。

在解釋生物學原理和自然種群規(guī)律方面，非線性模型均優(yōu)于線性模型，非線性模型已經(jīng)廣泛應用于昆蟲生長發(fā)育和種群動態(tài)研究。早期研究人員對昆蟲種群動態(tài)研究以描述性為主，并進行種群聚集與擴散遷移的研究，現(xiàn)在則應用模型模擬來研究昆蟲與環(huán)境之間的相互作用關系，目的是為了加強理解昆蟲種群的生命系統(tǒng)，改進對昆蟲種群的監(jiān)測調(diào)查和防治的具體措施。而關于斑翅果蠅種群動態(tài)模型研究的文章較少，例如，整合與溫度相關的生命表數(shù)據(jù)并通過矩陣投影模型來估算斑翅果蠅的種群動態(tài)(Wimanetal.，2014)；或應對入侵北美地區(qū)的斑翅果蠅開發(fā)了一個連續(xù)的時間階段種群結構模型(Langilleetal.，2016)。Hastie和Tibshirani(1986)提出并發(fā)展了廣義可加模型(generalized additive models, GAM)，該模型可擬合響應變量與解釋變量之間的非線性關系，但模型需要考慮曲線的擬合優(yōu)度、光滑度，需要克服解釋變量之間共曲線性關系導致的系數(shù)標準誤偏差，并利用兩個解釋變量之間的Pearson相關系數(shù)R來判別兩者之間的相關性程度。若要求得到的擬合曲線平穩(wěn)變化，不出現(xiàn)非劇烈快速波動的情況，則需要通過響應變量與解釋變量之間的關系曲線粗糙度來度量擬合光滑度(歐陽芳和戈峰，2013)。

本研究在昆明市云南農(nóng)業(yè)大學后山設置誘捕點，調(diào)查統(tǒng)計2016年至2017年連續(xù)兩年斑翅果蠅成蟲的種群動態(tài)變化，明確調(diào)查期內(nèi)斑翅果蠅成蟲種群數(shù)量動態(tài)，分析斑翅果蠅成蟲種群數(shù)量分布，通過廣義可加模型(GAM)分析斑翅果蠅成蟲種群數(shù)量與溫濕度因子間的關系，并判別是何種關鍵因子顯著影響斑翅果蠅成蟲種群數(shù)量動態(tài)變化。不僅從生態(tài)學角度探尋了溫濕度因子對斑翅果蠅成蟲種群動態(tài)消長規(guī)律起到重要作用，也為后續(xù)研究斑翅果蠅產(chǎn)生何種適應性應對環(huán)境因子變化提供了重要的依據(jù)。

1 材料與方法

1.1 試驗地點

2016至2017年，在云南農(nóng)業(yè)大學校園后山自然果林(N25°8′16″，E102°45′15″，海拔1 970 m)進行連續(xù)兩年的斑翅果蠅成蟲種群動態(tài)調(diào)查。所選果林樹勢良好，主要以梨樹為主，夾雜少量桃樹，果林無化學藥劑噴施，地面生長有龍葵、蛇莓等野生漿果。

1.2 田間誘捕及溫濕度記錄

誘劑配制比例參照優(yōu)化配比，即糖(g)∶酒(mL)∶醋(mL)∶水(mL)=1∶3∶1∶6(伍蘇然等，2007)。采用五點取樣法，將自制誘捕器放置在離地面1.2～1.5 m處，不同取樣點之間距離至少間隔10 m。溫濕度記錄儀(GGL-20)由深圳宇問加壹傳感系統(tǒng)有限公司生產(chǎn)提供，懸掛位于果林較為隱蔽處。

春、冬兩季因自然環(huán)境氣溫較低，誘劑不易揮發(fā)，故收集頻率為1周一次，夏、秋兩季因自然環(huán)境氣溫較高，誘劑易揮發(fā)，故收集頻率為 2或 3 d一次。將誘捕到的所有昆蟲樣本帶回實驗室，在解剖鏡下進行分類鑒定，并統(tǒng)計斑翅果蠅雌雄成蟲數(shù)量。

1.3 廣義可加模型

廣義可加模型分析可分為以下3個步驟，依次為：變量預分析、模型構建、模型預測。廣義可加模型是根據(jù)具體情況將廣義線性模型和可加模型結合起來，其中，廣義線性模型是改變響應變量條件期望的函數(shù)形式，將它記為g(muY)，其中muY=E(Y/X1,…,Xn)，可表示為：g(muY)=α0+α1X1…+αnXn+ε，還可用非參數(shù)的形式來描述響應變量條件期望與解釋變量的對應關系，用f(x)表示，可表示為：E(Y/X)=f(x)，將其推廣到多個解釋變量時，則記為可加模型E(Y/X1，…，Xn)=f(X1)…+f(Xn)。將兩種方式結合起來，則廣義可加模型可表達為：g(muY)=α0+f1(X1)…+fn(Xn)+ε。其中，在g(muY)項中，muY是Y的期望值，g(…)是連接函數(shù)，α0是截矩，在fn(Xn)項中，fn(…)是解釋變量Xn的單變量函數(shù)，ε為隨機變量。模型評價指標選用廣義交叉驗證值(Generalized Cross-Validation, GCV)，GCV值最小的模型最優(yōu)(歐陽芳和戈峰，2013)。

1.4 數(shù)據(jù)分析

使用Excel對數(shù)據(jù)進行整理和預處理，利用R(version 4.0.5)對斑翅果蠅種群數(shù)量數(shù)據(jù)進行對數(shù)預處理，mgcv數(shù)據(jù)包擬合廣義可加模型并進行分析。

2 結果與分析

2.1 成蟲田間種群發(fā)生動態(tài)

2016年至2017年，兩年中全年均有斑翅果蠅發(fā)生，斑翅果蠅成蟲種群數(shù)量均從6月開始明顯增加至10月種群數(shù)量回落，高峰期分別出現(xiàn)在7月、8月、9月，在11月至次年5月斑翅果蠅成蟲種群數(shù)量都維持在較低水平，尤其在2-5月種群數(shù)量呈全年最低水平。雌、雄成蟲的種群動態(tài)變化也呈現(xiàn)先升高后降低的趨勢；其中，雄成蟲種群動態(tài)變化與總種群動態(tài)變化趨勢一致。2016年種群變化趨勢呈單峰型，而2017年呈雙峰型；其中，雌成蟲種群變化趨勢在兩年中均呈單峰型，最高峰分別出現(xiàn)在2016年7月、2017年8月。2017年斑翅果蠅成蟲數(shù)量，尤其是高峰期數(shù)量相較于 2016年明顯降低。2016年、2017年最高峰期斑翅果蠅成蟲總量分別為587頭、236頭，雄成蟲數(shù)量分別為356頭、130頭；雌成蟲數(shù)量分別為312頭、106頭(圖1)。

圖1 2016-2017年斑翅果蠅成蟲種群動態(tài)Fig.1 Population dynamics of Drosophila suzukii adults from 2016 to 2017

2.2 廣義可加模型分析溫濕度因子與斑翅果蠅成蟲種群數(shù)量之間的關系

2.2.1變量預分析

對2016-2017年斑翅果蠅雌成蟲、雄成蟲和總種群數(shù)量取對數(shù)并繪制QQ圖和頻次分布圖(圖2、圖3、圖4)，并初步判定種群數(shù)量(響應變量)的分布特征為正態(tài)分布類型，Identity link可作為連接函數(shù)，進而將種群數(shù)量作為相應變量Y，溫度(Temperature)、濕度(Humidity)因子作為解釋變量構建模型。

圖2 斑翅果蠅雌成蟲種群數(shù)量正態(tài)QQ圖和頻次分布圖Fig.2 Normal QQ plot and frequency distribution diagram of Drosophila suzukii female adults

圖3 斑翅果蠅雄成蟲種群數(shù)量正態(tài)QQ圖和頻次分布圖Fig.3 Normal QQ plot and frequency distribution diagram of Drosophila suzukii male adults

圖4 斑翅果蠅總種群數(shù)量正態(tài)QQ圖和頻次分布圖Fig.4 Normal QQ plot and frequency distribution diagram of Drosophila suzukii total adults

2.2.2模型構建

分析調(diào)查周期內(nèi)溫度與濕度因子間Pearson相關系數(shù)，得知兩個解釋變量間Pearson相關系數(shù)小于0.5，故本研究將斑翅果蠅雌成蟲、雄成蟲、總種群數(shù)量作為響應變量，溫度(Temperature)、濕度(Humidity)作為解釋變量構建廣義可加模型。對于多變量數(shù)據(jù)集，模型構建時可構造一個或多個可行模型，以不同解釋變量的組合構成多個分析模型。因此，初步構建響應變量與解釋變量間的分析模型是包括多個可能模型的模型組，再通過模型評估篩選最優(yōu)擬合結果。

g(Female)=αF+f1F(Temperature)+f2F(Humidity)+εF

g(Male)=αM+f1M(Temperature)+f2M(Humidity)+εM

g(Total)=αT+f1T(Temperature)+f2T(Humidity)+εT

其中，g(Female)、g(Male)、g(Total)分別是指斑翅果蠅雌成蟲、雄成蟲以及總種群動態(tài)的連接函數(shù)；α是函數(shù)截矩；在fn(Xn)項中，fn(…)是解釋變量Xn的單變量函數(shù)；ε是隨機變量。3個下角標F、M、T分別是指雌成蟲、雄成蟲以及總種群。

模型參數(shù)擬合結果表明，溫濕度因子顯著影響了2016-2017年斑翅果蠅田間成蟲種群動態(tài)(P<0.05)，雌成蟲、雄成蟲、總種群數(shù)量的GAM模型參數(shù)α分別為2.21346、2.55606和3.08316；雌成蟲、雄成蟲和總種群數(shù)量與溫濕度因子關聯(lián)的GCV值分別為0.57299、0.74501和0.64611，其中雌成蟲的GCV值最小，因此說明斑翅果蠅雌成蟲種群數(shù)量與溫濕度因子進行擬合的模型最優(yōu)(表1)。

表1 基于溫濕度因子的斑翅果蠅種群動態(tài)廣義可加模型參數(shù)擬合及檢驗

2.2.3模型預測

斑翅果蠅雌成蟲種群動態(tài)與溫濕度因子之間呈非線性相關，其中雌成蟲種群發(fā)生在溫度23℃以下呈局部負相關，在23℃以上呈局部正相關且上漲趨勢較為平緩；雌成蟲種群發(fā)生與濕度呈正相關，但在濕度變幅大的情況下種群發(fā)生增量小(圖5)。斑翅果蠅雄成蟲的種群動態(tài)與溫濕度因子之間呈非線性相關，雄成蟲種群發(fā)生在溫度23℃以下呈局部負相關，在23℃以上呈局部正相關，且上漲趨勢較為明顯；雄成蟲種群發(fā)生與濕度呈正相關，但在濕度變幅大的情況下種群發(fā)生增量小，斑翅果蠅總種群發(fā)生與雄成蟲種群發(fā)生表現(xiàn)基本一致(圖6，圖7)。

圖5 斑翅果蠅雌成蟲種群數(shù)量與溫濕度因子之間的相關性Fig.5 Population sizes of Drosophila suzukii female adults correlated with temperature-humidity factor注：Y軸標題括號中的數(shù)值為模型擬合的光滑度參數(shù)，該數(shù)值接近1，預測變量與目標變量間的相關關系越接近水平直線，反之曲線呈現(xiàn)不同相關關系。圖6-圖7同。Notes : The value in the brackets of Y-axis heading is the smoothness parameters of the prediction model. If the smoothness parameter is close to 1, the relationship between the predictor variables and dependent variables would be closer to a horizontal straight line, otherwise, the relationship would be different correlations. Same to Fig.6-Fig.7.

圖6 斑翅果蠅雄成蟲種群數(shù)量與溫濕度因子之間的相關性Fig.6 Population sizes of Drosophila suzukii male adults correlated with temperature-humidity factor

圖7 斑翅果蠅總種群數(shù)量與溫濕度因子之間的相關性Fig.7 Population sizes of Drosophila suzukii total adults correlated with temperature-humidity factor

3 結論與討論

3.1 斑翅果蠅成蟲全年種群數(shù)量動態(tài)規(guī)律及分析

種群動態(tài)調(diào)查是昆蟲種群生態(tài)學研究的核心問題，尤其是為研究昆蟲田間種群的適應性在環(huán)境因子變化下發(fā)生何種改變提供了基礎和前提。國外已有對斑翅果蠅的種群動態(tài)變化及寄主的調(diào)查研究，例如，在比利時對斑翅果蠅種群發(fā)生進行調(diào)查，結果表明2014年的早春(非常溫和的冬季之后)，已經(jīng)發(fā)現(xiàn)了斑翅果蠅的存在(Belienetal.，2014)；研究發(fā)現(xiàn)從1月開始，當休眠狀態(tài)下的斑翅果蠅雌成蟲有效積溫達到50～800℃·d時，則可在季節(jié)早期進行繁殖，而且可能會受到早期寄主植物可用性的限制(Wiman etal.，2016)。而斑翅果蠅在國內(nèi)的研究主要集中在果園內(nèi)進行種群動態(tài)的監(jiān)測調(diào)查以及為害特性觀察的研究，例如，調(diào)查研究發(fā)現(xiàn)陜西灞橋櫻桃果園里斑翅果蠅為優(yōu)勢種(郭麗娜等，2014)；陜西楊凌櫻桃園中斑翅果蠅種群數(shù)量在10月中旬達到高峰期(?；廴旱?，2017)；在斑翅果蠅發(fā)生的高峰期，云南石屏果園中的楊梅果實被害率可達100%(趙超等，2017)。在本研究中，2016-2017年調(diào)查發(fā)現(xiàn)斑翅果蠅成蟲的種群數(shù)量通常在野生果林的果實進入成熟期后大幅增長，高峰期集中在 7月至9月，與不同引誘劑對貴州麻江田間藍莓果蠅誘捕結果基本一致(徐芳玲等，2012)，即8月藍莓果實進入盛果期，成熟的藍莓對果蠅的誘集數(shù)量呈現(xiàn)高峰值；在2月至5月，斑翅果蠅成蟲種群數(shù)量動態(tài)均維持在較低水平，與陜西楊凌櫻桃園中斑翅果蠅發(fā)生特點基本一致(楊靜，2017)；斑翅果蠅雄成蟲以及總種群數(shù)量動態(tài)變化趨勢基本同步，其中總種群發(fā)生動態(tài)與王華弟等研究浙江地區(qū)楊梅果蠅種群發(fā)生動態(tài)時監(jiān)測得出的結果基本一致(王華弟等，2017)。

3.2 溫濕度因子與斑翅果蠅成蟲種群動態(tài)之間的相關性分析

溫濕度變化對昆蟲的生長發(fā)育速率、生殖能力及存活率等種群數(shù)量構成的關鍵因子產(chǎn)生重大影響，從而導致昆蟲的種群數(shù)量發(fā)生變化(陳瑜和馬春森，2010)。此前有研究報道，結合2016-2018年氣溫和降水量資料，分析大豆蚜蟲Aphisglycines田間動態(tài)與氣溫和降水量的總體關系，結果得出氣溫以20℃為節(jié)點起伏，田間溫度越高大豆蚜蟲發(fā)生量越大，而降水量對大豆蚜蟲田間種群動態(tài)影響不明顯(邵珊珊等，2019)。滇中菜區(qū)小菜蛾Plutellaxylostella成蟲種群消長動態(tài)與溫度顯著相關r=0.442，與降雨量成正相關r=0.283，但相關程度弱(趙雪晴等，2016)。假眼小綠葉蟬Empoascavitis的成蟲蟲口數(shù)量受田間氣候條件影響較大，其中溫度影響最大(李慧玲和林乃銓，2012)。本研究中以溫度23℃為節(jié)點，斑翅果蠅成蟲種群發(fā)生數(shù)量與23℃以下呈局部負相關，與23℃以上呈局部正相關，與濕度呈正相關，但在濕度變幅大的情況下種群發(fā)生增量小，得出濕度對斑翅果蠅的種群影響較溫度的影響低，溫度是影響斑翅果蠅成蟲種群數(shù)量動態(tài)的關鍵生態(tài)因子，這一結果與上述溫濕度因子對其他昆蟲種群發(fā)生的影響研究結果趨向一致。

本研究中斑翅果蠅種群數(shù)量在高峰期的時段，正處于氣溫較高的夏季，同時果園也具備了數(shù)量、品種較多的寄主植物，而在種群數(shù)量維持在較低水平的其它季節(jié)，寄主植物的數(shù)量及品種也減少，后續(xù)深入探索可參考寄主等因子對斑翅果蠅種群發(fā)生的影響研究，例如，墨西哥果實結實季節(jié)(8月至次年5月)是影響當?shù)靥镩g斑翅果蠅季節(jié)性活動的最主要因素之一(Dominguez-Mongeetal.，2021)；相同溫度條件下，豇豆蚜蟲Aphiscraccivora對于大豆(首選寄主)和棉花(次選寄主)表現(xiàn)不同，即在大豆中的內(nèi)在增長率、繁殖率、每成年子代數(shù)和壽命均顯著高于棉花，且溫度波動到極端水平時，導致以棉花為食的蚜蟲的上述各項參數(shù)迅速下降，而以大豆為食的蚜蟲則沒有下降，在一定程度上解釋了夏初新疆棉區(qū)氣溫升高時豇豆蚜蟲從棉花向大豆進行生態(tài)位轉變的現(xiàn)象(Luetal.，2017)。因此，在全球變暖的背景下，除了對害蟲進行種群動態(tài)建模和風險分析之外，同時也應探索寄主等因子影響害蟲如何調(diào)節(jié)自身應對極端環(huán)境，或發(fā)生了何種生態(tài)位轉變的研究。