基于COⅠ基因的我國懸鈴木方翅網蝽(半翅目： 網蝽科) 種群遺傳多樣性和遺傳結構分析*

張元臣 盧紹輝 龔東風 馬興莉

(1.安陽工學院生物與食品工程學院 安陽 455000； 2.河南省太行山林業有害生物野外科學觀測研究站 林州 456550；3.河南省林業科學研究院 鄭州 450002； 4.北卡羅來納州立大學農業與生命科學學院有害生物綜合治理中心 羅利27695)

物種遺傳變異與其適應環境的能力和進化程度密切相關，也是物種在新環境下能夠定殖和遷移的關鍵因素，而生物入侵能夠引起種群遺傳變異，這為研究入侵生物的生態和進化問題提供了理想的模式(Amourouxetal., 2013; Lau, 2015)。分子標記是研究物種遺傳變異非常有效的工具，能夠為物種的起源和分化提供直接證據，可揭示入侵物種的擴散路徑和結果，并推測未來的種群發生情況(李菁等, 2010; 賓淑英等, 2014; 楊順義等, 2018)。目前已經開發出了多種分子標記技術，而線粒體基因具有進化迅速、垂直母系遺傳、重組困難等特點，使其廣泛應用于昆蟲種群的遺傳變異和遺傳結構的研究中(王克勤等， 2018; Cameron, 2014; Menetal., 2017)。很多研究已經利用昆蟲線粒體細胞色素氧化酶亞基Ⅰ(cytochrome C oxidase subunit Ⅰ，COⅠ)基因對種群遺傳分化、入侵路徑和入侵來源等進行了深入研究(王蒙等, 2014; 王戎疆等， 2018; Seraphimetal., 2016; Ariasetal., 2018; Chinetal., 2018; Aslametal., 2019)。

懸鈴木方翅網蝽(Corythuchaciliata)(半翅目Hemiptera: 網蝽科Tingidae)(以下簡稱網蝽)具有快速適應新棲息地的特點，在新的生境里經常大面積暴發成災(盧紹輝等, 2018; Maceljskietal., 1974)； 目前已經成為懸鈴木屬(Platanus)樹木上的首要害蟲(Bures, 1997)，防治非常困難(Cravedi, 2000)。2002年在長沙首次發現該網蝽(Streito, 2006)， 2006年開始在武漢市區大面積暴發(李傳仁等, 2007)，此后，該網蝽在長江中下游的多個區域相續被發現，嚴重為害懸鈴木屬植物，給城市綠化帶來了非常大的危害(Juetal., 2009)。自該網蝽入侵我國以來，昆蟲研究者對其生態生物學特性、防治方法、重要基因的功能等進行了深入的研究(王福蓮等, 2008; 郝德君等, 2012; 付寧寧等, 2017; 李峰奇等, 2018; Juetal., 2011)。Yang 等(2014)發掘了9個該網蝽微衛星標記位點，并對這9個微衛星標記位點在不同地理種群中的有效性進行了驗證。進一步利用這些微衛星標記位點和3個線粒體基因的部分核苷酸序列，對收集于2009年到2014年間我國的21個網蝽地理種群和國外的2個網蝽地理種群進行了遺傳結構研究，結果發現這些種群間具有高度的遺傳分化，并推測該網蝽從東海岸入侵我國，之后向全國范圍進行擴散(Yangetal., 2017)。然而，隨著時間的推移，懸鈴木方翅網蝽在我國入侵范圍急劇擴大，因此有必要對該網蝽遺傳結構進行分析和擴散路徑的重建。基于此，本研究采用線粒體COⅠ標記基因對采集于2015年我國12個不同地理種群網蝽的遺傳多樣性與遺傳結構進行分析，旨在探討網蝽在中國不同地理種群間的遺傳變異情況、分化程度及基因流水平，據此分析遺傳分化產生的原因及在我國擴散的中心區域，并進行了入侵趨勢的分析，從而為有效防治網蝽提供基礎理論依據。

1 材料與方法

1.1 試驗材料的采集

選擇網蝽為害的懸鈴木，隨機取下一片葉片，用毛筆將葉片上的網蝽成蟲挑到含有100%酒精的2 mL離心管中，帶回室內用顯微鏡進行確認，之后將確定的網蝽成蟲樣品放到-20 ℃冰箱中保存。每個點隨機選取間距在200 m以上的懸鈴木20棵，每一棵樹收集1頭成蟲，每個收集點共20頭網蝽成蟲被收集。網蝽成蟲收集時間在2015年7月，收集地分別在長沙、武漢、上海、北京、石家莊、濟南、太原、杭州、成都、重慶、貴陽和鄭州共12個城市。

1.2 基因組DNA的提取

將100%酒精保存的網蝽成蟲從離心管中挑出來，放到10 cm培養皿中，用ddH2O沖洗2次。將網蝽放到吸水紙上將水吸干，按照DNA提取試劑盒(DP304-03)(天根生化科技公司)提取DNA。用1.2%～1.5%的瓊脂糖電泳和NanoDrop(賽默飛世爾科技公司)檢測基因組DNA的質量和濃度后，將其稀釋同樣的濃度(30 μg·μL-1)，放置于-20 ℃條件下保存備用。

1.3 COⅠ基因擴增、測序及拼接

根據Yang等(2017)的方法對COⅠ基因進行擴增，將PCR產物送往上海生工進行雙向測序。將測序獲得的核苷酸序列去除掉引物，剩余核苷酸序列放到NCBI數據庫中進行BLASTn比對，根據比對結果驗證網蝽COⅠ基因核苷酸序列的正確性，隨后以網蝽COⅠ核苷酸序列為基礎進行后續的數據分析。

1.4 數據分析

參考呼曉慶等(2019)的方法對不同地理種群網蝽線粒體COⅠ序列進行核苷酸組分、保守位點、多態性位點、單倍型組成、單倍型數量、單倍型分布(Hap)、單倍型多樣性(haplotype diversity,Hd)、核苷酸多樣性(nucleotide diversity,π)、核苷酸平均差異數(average number of nucleotide differences,k)、不同地理種群間的固定系數FST，并依據(1-FST)/2FST公式計算種群間基因流動值(Nm)。根據Qin等(2016)的方法構建了9個網蝽單倍型網狀進化圖，并進行網蝽不同種群間的遺傳距離與地理距離的相關性分析。為了方便對網蝽在主要分布區的遺傳多樣性進行分析，本文根據不同地理種群間地理位置將主要分布區分為華中組群(長沙、武漢)、華東組群(上海、杭州)、華北組群(北京、石家莊、太原、鄭州和濟南)、西南組群(成都、重慶和貴陽)。對4個組群進行分子方差分析(analysis of molecular variance,AMOVA)，并計算不同組群網蝽間的FST值和Nm值，同時對不同組群間進行種群歷史動態分析得到擴張前后網蝽有效種群的大小(region mutation parameter,θ)、網蝽種群的擴張時間(T)、Tajima’sD值、Fu’sFS值及期望錯配分布值與觀察錯配分布值之間的平方和(SSD)。通過MIGRATE V.3.2.1軟件(Beerli, 2006)計算網蝽組群有效遷移率(migration rate,M)； 利用MEGA v7.0.26軟件根據本研究中12個網蝽種群和NCBI上其他國內外網蝽種群的COⅠ基因核苷酸序列，基于Kimura 2-Parameter模型計算距離矩陣并根據鄰接法建立系統發育樹。

2 結果與分析

2.1 網蝽COⅠ基因核苷酸的堿基構成及序列

本研究共獲得了240條761 bp的網蝽COⅠ基因序列，A、T、G和C核苷酸含量分別為35.0%、32.6%、17.5%和14.9%。A+T核苷酸含量(67.6%)明顯高于G+C核苷酸含量(32.4%)，核苷酸組成類型符合昆蟲線粒體核苷酸的特征。240個COⅠ序列含有13個多態性位點(polymorphic sites)，為總堿基數的1.71%。發現的多態性位點全部是核苷酸間的相互替換引起的，含有4個自裔位點和9個簡約信息位點。

2.2 網蝽地理種群的單倍型和核苷酸多樣性

在240條網蝽COⅠ基因核苷酸中發現了9個單倍型，分別有4個共享單倍型和5個獨享單倍型(表1)，表明網蝽地理種群間已經出現了遺傳分化。所有單倍型中分布最多的是Hap1和Hap3單倍型，都被9個種群共享，而Hap4和Hap7單倍型分別被3個和2個種群共享，其他單倍型為獨享單倍型。12個地理種群單倍型在1～6個之間，其中GY種群擁有最高的單倍型數量(6個)，而SH和CS種群擁有最低的單倍型數量(1個)(表1)。單倍型網絡圖表明，Hap1單倍型居于網絡圖的中心，與其相關聯的單倍型數量最多(圖1)。Hap3單倍型共享個體數量也較多，然而與其相關聯的單倍型數量比較少，而且Hap3單倍型是從Hap7單倍型分化出來的，Hap7單倍型直接從Hap1單倍型分化而來。基于此，本研究認為Hap1為祖先單倍型。

表1 網蝽COⅠ基因單倍型的分布①Tab.1 Distribution information of haplotype from COⅠ gene in C. ciliata

圖1 網蝽單倍型網絡進化Fig. 1 Haplotype Median-Joining network using COⅠ gene from C. ciliata

網蝽總種群單倍型多樣度Hd值為0.584，核苷酸多樣度π值為0.003，核苷酸平均差異數k值為2.276(表2)。不同地理種群間Hd值在0.000～0.758之間；π值在0.000～0.005之間；k值在0.000～3.853之間(表2)。在分析的12個地理種群中，GY種群擁有最高的單倍型遺傳多樣性(Hd值為0.758)，而CS和SH種群擁有最低的單倍型多樣性(Hd值均為0)(表2)。總體上，12個網蝽地理種群中遺傳多樣性差異比較明顯，CS、SH種群無多態性，12個種群遺傳多樣性從中國南部——貴陽到華中地區——鄭州等地區逐漸減少，而從華中到華北——北京等地區遺傳多樣性逐漸增大。

表2 12個網蝽地理種群COⅠ 基因單倍型多樣度及核苷酸多樣性①Tab.2 Haplotype diversity and nucleotide diversity of COⅠ gene in twelve geographic populations of C. ciliata

2.3 系統發育

種群系統發育樹結果表明，法國、斯洛文尼亞與德國種群聚集在一起，表明這3個區域網蝽種群的親緣關系非常近(圖2)。在國內種群中，CS、SH、WH種群與這3個歐洲種群親緣關系也非常近(圖2)，如果我國網蝽來源于歐洲種群，CS、WH和SH都有可能是我國網蝽的首次入侵地。

圖2 根據COⅠ基因核苷酸序列建立的網蝽種群系統發育樹Fig. 2 Phylogenetic tree of the nucleotide sequence of COⅠ gene from C. ciliata populations

其他國內外網蝽種群并未呈現出地理區域的相關性，甚至在同一區域內聚集有多個遺傳關系比較遠的多個種群，例如斯洛文尼亞1種群與韓國種群聚在一支，而4種群與法國種群聚在一支，5種群卻與我國保定1種群聚在一支。作為網蝽發源地的美國種群與國內保定1種群和揚州種群聚集在一支，說明它們的親緣關系比較近(圖2)。

2.4 分子變異

組群間AMOVA分析發現，其在組群間、組群內種群間和種群內個體間的固定系數，即FCT、FSC和FST的P值都小于0.05，說明上述分組群進行分析是合理的。組群間、組群內種群間和種群內個體間變異分別占總變異的26.03%、21.78%和52.19%(表3)。這樣的結果說明，本研究中網蝽種群間變異主要是由種群內個體之間的變異引起的。

表3 12個網蝽地理種群的分子方差分析①Tab.3 AMOVA analysis of twelve geographic populations of C. ciliata

2.5 網蝽種群間的遺傳關系

種群間遺傳分化程度及基因流動分別采用FST及Nm來衡量，FST和Nm值代表的意義參考盧紹輝(2020)。4個組群間，華中與華東組群間網蝽具有較低的遺傳分化，而華中與西南、華北組群間網蝽具有非常高的遺傳分化，而華東與西南、華北組群間網蝽也具有較高程度的遺傳分化(表4)。在66個種群對間FST值在-0.05～0.89中，其中有43個種群對間有顯著差異(表5)。有10個種群對間的FST值小于0.05(表5)，說明這些種群對間遺傳分化程度很低，甚至沒有發生遺傳分化，而其他33個種群對間的FST值大于0.05(表5)，說明這些種群對間遺傳分化程度比較大。不同地理種群對間Nm值范圍在-85.25～47.58間，有10個種群對間的Nm值大于4(表5)，暗示這些種群對間的基因交流經常發生； 17個種群對間Nm值在1～4間(表5)，表明這些種群對間有中等程度的基因交流； 其他種群對間的Nm值均小于1(表5)，表明種群對間幾乎不發生基因交流。

表4 網蝽組群對的FST值(對角線下)和Nm值(對角線上)①Tab.4 Pairwise FST(below the diagonal) and Nm(above the diagonal) among group populations of C. ciliata

表5 12個網蝽地理種群對間的FST值(對角線下)和Nm值(對角線上)Tab.5 Pairwise FST(below the diagonal) and Nm(above the diagonal) among twelve geographical populations of C. ciliata

2.6 網蝽遺傳距離與地理距離的相關性

對不同網蝽種群的遺傳距離與地理距離間進行相關性分析，發現二者之間相關性系數為0.15，但不存在顯著差異(圖3)。說明網蝽不同地理種群間的遺傳分化不是由地理隔離引起的，沒有距離隔離的現象。網蝽種群位置的經度和緯度與單倍型多樣性的相關系數分別為-0.127 和0.103，都沒有表現出顯著的相關性。

圖3 網蝽種群遺傳距離和地理距離間的相關性Fig. 3 Relevance between genetic and geographic distance from twelve populations of C. ciliata

2.7 不同區域種群歷史動態和有效遷移率

對4個區域種群的歷史動態進行了分析。Tajima’sD和Fu’sFS分別為負值且P<0.05時，顯著偏離了中性突變，認為在歷史上發生過種群擴張事件，否則認為種群在歷史上沒有擴張事件；PSSD>0.05，表明不能拒絕群體擴張的假說，即符合原來群體發生了擴張，而PSSD<0.05說明種群偏離突然擴張假設模型，即原來群體不能發生擴張(Tajima, 1989)。華東和華中組群內只有2種單倍型，因此這2個組群不能繪制錯配圖。在主要分布區，θ0和θ1值沒有顯著的改變，Tajima’sD及Fu’sFS值均不顯著(表6)，說明中國區域網蝽種群沒有發生擴張事件。同時，錯配分布圖有非常顯著的雙峰(圖4C)以及PSSD值小于0.05(表6)，也證實中國區域網蝽種群沒有發生擴張事件。

表6 網蝽種群歷史分析①Tab.6 Demographic history parameters of C. ciliata

圖4 中國主要分布區網蝽的錯配分布Fig. 4 Mismatch distributions of the main distribution range of C. ciliata in ChinaA： 西南地區 Southwest China； B： 華北地區 North China； C： 中國主要分布區 Main distribution areas in China.

華北、西南和華東組群網蝽擴張前有效種群θ0和擴張后有效種群θ1間數值變化不大(表6)，表明這3個組群網蝽數量增加不明顯，且Tajima’sD和Fu’sFS值均沒有顯著性(表6)，呈現雙峰或者是下降型的錯配分布圖(圖4A和B)，表明這3個地區網蝽未經歷過明顯的種群擴張(盧紹輝， 2020)。華中區域網蝽θ1值比較大，表明該地區網蝽種群增加較大(表 6)。同時，Tajima’sD為負值且Fu’sFS為顯著的負值，表明該區域網蝽種群經歷過擴張事件(表 6)。同時，PSSD值大于0.05(表 6)，也證實了該組群網蝽經歷過種群擴張事件。

4個網蝽區域組群間的有效遷移率(M)都比較高，在421.9～648.0間(表7)。華中組群向華東、華北、西南的M值分別為580.9、610.4和598.4，而華東、華北、西南組群向華中組群的M值分別為449.8、421.9和455.7，表明華中組群遷出要高于遷入(表7)，這與網蝽種群歷史分析中得到的“華中地區種群經歷過種群擴張的事件”結論一致。

表7 中國網蝽組群間的遷移率①Tab.7 Immigration rate among group populations of C. ciliata in China

3 討論

本研究對收集的12個地理種群的240個懸鈴木方翅網蝽個體進行COⅠ基因測序，獲得了240條761 bp的核苷酸序列，根據這些序列進一步分析了不同地理種群間的核苷酸序列差異及遺傳分化。本研究測定的COⅠ基因核苷酸序列間未發現有堿基缺失和插入的情況。腺嘌呤(A)和胸腺嘧啶(T)的含量(67.6%)明顯高于鳥嘌呤(G)和胞嘧啶(C)的含量(32.4%)，表現出A+T堿基偏倚性，符合昆蟲線粒體構成的典型特征(Jermiinetal., 1994)。

多種因素影響生物向新棲息地擴散轉移，通常這些因素也都可以導致生物種群發生遺傳分化(陳圣賓等， 2010)。種群分化程度與其棲息地的距離存在著正相關，距離越遠分化程度越大，這種現象稱為距離隔離(Kimuraetal., 1964； Slatkin, 1993)。然而，本文發現中國分布區內的網蝽不同地理種群間并不存在距離隔離，說明中國區域網蝽發生的遺傳分化不是由距離隔離產生的。對其他昆蟲研究也發現了類似的結果，如草地螟(Loxostegesticticalis)(呼曉慶等， 2019)。研究發現，一些山脈和高原，比如云貴高原等，能明顯影響物種的遺傳結構(Yeetal., 2014; Xieetal., 2017)。GY、CQ和CD種群位于我國云貴高原區域，該區具有復雜的地形和獨特的氣候特點，被認為是生物高度富集的重要區域(Barthlottetal., 2007)。西南區域的GY、CQ和CD網蝽種群，具有很高的單倍體多樣性，可能是因為更適應這里的環境所造成的，因為生物在適宜的棲息地往往有高的遺傳多樣性(Petitetal., 2003； 盧紹輝, 2020)。華中區域的WH、CS種群和華東區域的SH、HZ種群具有非常少的單倍型個數，并且2種群間未發現遺傳分化。網蝽可以通過人類活動傳播(鞠瑞亭等， 2011)，長江中下游地區是我國經濟高度發達的區域之一，人員和貨物流動頻繁(姚瑞華等， 2014)，因此華中和華東組群間沒有遺傳分化很有可能是因為人類活動導致經常發生基因交流造成的。另外，華中和華東區域處于我國長江中下游，湖泊河流密集。每年的夏季為網蝽的發生高峰期(李峰奇等， 2018； 盧紹輝, 2020)，也是長江中下游地區洪水頻繁發生的時期(許有鵬等， 2005; 劉蕾等， 2018)，網蝽有可能以懸鈴木葉片為載體順著河流進行擴散，加快了網蝽種群間基因交流。因此，地理障礙、人為傳播和氣候條件等都能影響種群之間的基因交流，從而導致種群間出現遺傳分化。

網蝽地理種群之間的單倍型多樣數量隨著種群棲息地經度的減小而增加，且位于華東和華中組群具有非常低的單倍型多樣性。研究發現，瓶頸和自然選擇作用通常能影響入侵生物的遺傳多樣性，往往初始種群有更高的多態性，相反有比較低的多態性(Lanzavecchiaetal.， 2008; 盧紹輝, 2020)。筆者的研究結果表明網蝽在西南地區具有最高的遺傳多態性，按此推理，西南種群應該為原始種群。然而，對4個組群的歷史動態和錯配圖進行分析發現，西南、華北和華東組群并沒有發生過種群向外擴張事件，僅僅華中組群發生了種群向外擴張事件。華中組群遷出率遠高于遷入率，也暗示華中組群經歷過種群遷出事件。結合網蝽2002年在華中地區首次報道的事實(Streito, 2006)，筆者推測長沙和武漢可能是蟲源地之一。同時系統發育分析表明，WH、SH和CS種群聚集在一支，具有非常高的親緣關系，因此上海也可能是蟲源地。Yang等(2017)也發現網蝽首先入侵中國的東海岸，之后向全國擴散，結合這項研究結果，本研究中WH、SH和CS種群僅有SH種群居于東海岸，那么更有可能上海為網蝽的蟲源地，并向四周擴散。根據系統發育分析結果，SH種群與歐洲種群的親緣關系最近，推測我國網蝽首先是從歐洲入侵的。上海是我國對外開放的最大窗口，與歐洲貿易頻繁發生，因此很有可能通過人類活動進行長距離傳播。之后以上海為蟲源地，通過人類活動的長距離傳播和飛行等短距離傳播在全國范圍內進行擴散(李峰奇等, 2018)。然而，網蝽在我國并沒有最先發現于上海(肖娛玉等， 2006; Streito, 2006)。入侵昆蟲在新發地初期常因種群數量有限而往往不容易被發現, 如白蠟窄吉丁(Agrilusplanipennis)(Siegertetal., 2014)、喜菘蝽(Bagradahilaris)(Palumboetal., 2016)和斑翅果蠅(Drosophilasuzukii)(Bieńkowskietal., 2020)。網蝽也很有可能在上海首先發生入侵后，由于種群數量有限而未被發現，導致了這種現象的出現。

SH、WH和CS種群對間沒有發生種群的分化(FST=0)，說明它們之間的親緣關系非常近。依據地理距離推測網蝽遷移路線是上海→武漢→長沙。研究認為網蝽能夠借助自然風向較遠距離進行傳播(李峰奇等, 2018; 盧紹輝, 2020)。夏季在副熱帶高壓的作用下，我國沿海地區風向以東南風為主(李慧等， 2013)。因此網蝽可能依靠東南風向西北方向遷移，擴散至武漢、長沙等地區。另外，網蝽廣泛分布于城市懸鈴木上，交通工具從懸鈴木樹下通過時也可能將網蝽運走，從而產生了遠距離的擴散。CS與GY、CQ、CD 種群FST值分別為0.36、0.42和0.79，說明西南地區網蝽種群可能是CS種群以貴陽為中轉站遷移進入的，形成了長沙→貴陽→重慶→成都的遷移路線。然而云貴高原地理障礙能力非常強，外部種群很難進入(Wangetal., 2020)。再加上長沙與貴陽距離遙遠，網蝽不太可能通過飛行進入貴陽，更可能是通過人為傳播進入貴陽地區。研究發現，人類活動，比如船、高速公路、苗木運輸等，都能遠距離傳播網蝽(Halbertetal., 1983; 虞國躍等, 2014; 鞠瑞亭等, 2011)。這些因素都可能導致網蝽擴散到我國西南地區。與其他種群相比，WH與ZZ種群間分化非常低(FST=0.10)，具有非常高的親緣性，說明ZZ種群很有可能從武漢傳入的。ZZ與SJZ、JN、TY種群間的FST值分別為0.69、0.57、0.70，說明ZZ種群與這3個地區的種群產生了顯著的遺傳分化。假設SJZ、JN和TY種群是從鄭州地區遷入形成的，之后擴散到北京，那么BJ與ZZ種群應該也有顯著的遺傳分化，然而研究結果表明，ZZ與BJ種群的FST值為0.15，說明2個種群間遺傳分化程度非常低。近年來國內外頻繁的貿易往來，很可能在這些地區發生了網蝽多次入侵事件，從而導致這種現象的出現。多次入侵事件的發生，導致很難推測網蝽從鄭州向北擴散的路線。因此，網蝽向西的主要遷移路線為上海→武漢→長沙→貴陽→重慶→成都； 向北的主要遷移路線為上海→武漢→鄭州，而鄭州向北遷移途徑需要進一步研究。

4 結論

本文對12個懸鈴木方翅網蝽種群遺傳結構進行了研究，在不同組群間和不同地理種群間具有顯著分化的遺傳結構，而網蝽種群個體間的變異是種群發生遺傳分化的主要因素，同時地理障礙、人類活動和氣候條件均有可能影響網蝽種群間的遺傳分化程度。我國主要分布區內網蝽可能以上海為中心向西、向北往全國進行擴散。