高橋仁蚧線粒體全基因組測(cè)序與分析

鄧 鋆, 王 剛, 盧聰聰, 張江濤, 黃曉磊,*

(1.福建農(nóng)林大學(xué)植物保護(hù)學(xué)院, 閩臺(tái)作物有害生物生態(tài)防控國(guó)家重點(diǎn)實(shí)驗(yàn)室, 福州 350002;2.江西農(nóng)業(yè)大學(xué)林學(xué)院, 鄱陽(yáng)湖流域森林生態(tài)系統(tǒng)保護(hù)與修復(fù)國(guó)家林業(yè)和草原局重點(diǎn)實(shí)驗(yàn)室, 南昌 330045)

線粒體是真核生物細(xì)胞中具有獨(dú)立遺傳物質(zhì)的細(xì)胞器，是細(xì)胞制造能量的場(chǎng)所，參與細(xì)胞分化、細(xì)胞信息傳遞和細(xì)胞凋亡等過(guò)程(潘寶平和卜文俊, 2005; Cameron, 2014)。截至2020年，在GenBank數(shù)據(jù)庫(kù)中已發(fā)布489種半翅目(Hemiptera)昆蟲完整線粒體基因組。作為半翅目一員，蚧總科(Coccoidea)現(xiàn)已知35科8 310種蚧蟲(García Moralesetal., 2016)，但在GenBank數(shù)據(jù)庫(kù)中僅有3個(gè)蚧蟲完整的線粒體基因組序列，均為蚧科(Coccidae)物種(Dengetal., 2019; Luetal., 2020; Xuetal., 2021)。相比現(xiàn)有的蚧蟲物種數(shù)量，蚧蟲線粒體基因組的測(cè)序工作明顯滯后，GenBank數(shù)據(jù)庫(kù)中嚴(yán)重缺少除cox1基因外的其他蚧蟲線粒體標(biāo)記基因數(shù)據(jù)，這阻礙了分子進(jìn)化、系統(tǒng)發(fā)育和生物地理學(xué)等相關(guān)工作在蚧蟲中的開(kāi)展，也阻礙了線粒體基因組在半翅目昆蟲研究中進(jìn)一步的應(yīng)用。Deng等(2019)公布首個(gè)蚧蟲即日本龜蠟蚧Ceroplastesjaponicus完整的線粒體基因組序列，并于2020年完整注釋了蚧科咖啡黑盔蚧Saissetiacoffeae線粒體基因組，彌補(bǔ)了GenBank數(shù)據(jù)庫(kù)因長(zhǎng)期缺失蚧蟲線粒體基因組數(shù)據(jù)而造成的空白(Luetal., 2020)。研究發(fā)現(xiàn)蚧蟲線粒體基因組中存在高度的基因重排現(xiàn)象、高AT含量以及非典型的tRNAs二級(jí)結(jié)構(gòu)(例如縮減tRNAs，主要表現(xiàn)在缺失DHU臂或TΨC臂)(Luetal., 2020)。事實(shí)上，GenBank數(shù)據(jù)庫(kù)缺失的不僅僅是完整的蚧蟲線粒體基因組序列，數(shù)據(jù)庫(kù)中只有兩條超過(guò)5 000 bp的部分線粒體基因序列。這些現(xiàn)象反映想要獲得完整的蚧蟲線粒體基因組可能具有較高的難度。

高橋仁蚧Aclerdatakahashii隸屬于半翅目蚧總科仁蚧科(Aclerdidae)，分布在巴西、中國(guó)、埃及、印度、印度尼西亞、馬來(lái)西亞、美國(guó)等14個(gè)國(guó)家和地區(qū)，涵蓋古北界、新北界、東洋界、新熱帶界(García Moralesetal., 2016)。高橋仁蚧各齡期均生活在植物的葉鞘下，寄主主要為禾本科芒屬M(fèi)iscanthus、甘蔗屬Saccharum、粽葉蘆屬Thysanolaena植物(García Moralesetal., 2016)。本研究對(duì)高橋仁蚧粒體基因組進(jìn)行二代高通量測(cè)序，獲得其完整的線粒體基因組序列，對(duì)其堿基組成、密碼子使用情況、tRNA結(jié)構(gòu)進(jìn)行了預(yù)測(cè)分析，并基于線粒體基因組13個(gè)蛋白編碼基因構(gòu)建半翅目系統(tǒng)發(fā)育樹(shù)。研究結(jié)果將為蚧總科分類及仁蚧科線粒體基因組研究提供科學(xué)依據(jù)。

1 材料與方法

1.1 昆蟲樣本采集與DNA的提取

本研究所用的高橋仁蚧成蟲樣本于2019年8月5日采集自福建省建甌市云際山公園內(nèi)(27°1′19″N, 118°18′36″E)，利用《中國(guó)動(dòng)物志》(王子清, 2001)進(jìn)行形態(tài)鑒定。測(cè)序樣本置于95%的酒精中-80℃保存。使用DNeasy DNA Extraction Kit試劑盒(Qiagen, Hilden, 德國(guó))提取10頭成蟲混合樣本的總DNA，經(jīng)NanoDrop 2000分光光度計(jì)和1%的瓊脂糖凝膠電泳檢測(cè)總DNA的質(zhì)量和濃度。

1.2 高通量測(cè)序及線粒體基因組的拼接

將1.1節(jié)質(zhì)量合格濃度大于2 μg的DNA樣本送至北京諾禾致源科技股份有限公司進(jìn)行高通量測(cè)序，構(gòu)建350 bp的短片段文庫(kù)，利用Illumina HiSeq測(cè)序平臺(tái)進(jìn)行雙末端(paired-end, PE)測(cè)序。 將下機(jī)的原始數(shù)據(jù)去接頭和質(zhì)控。利用NOVOPlasty 3.7.1(Dierckxsensetal., 2017)軟件以該物種線粒體cox1基因片段為種子序列進(jìn)行線粒體基因組的拼接，在K-mer=23, 27, 33下均獲一條長(zhǎng)為16 599 bp的成環(huán)序列；為保證拼接的準(zhǔn)確性，使用MEGAHIT 1.0(Lietal., 2016)對(duì)clean data進(jìn)行從頭組裝，建立本地Blast數(shù)據(jù)庫(kù)，提取一條長(zhǎng)為16 735 bp的線粒體基因組序列；隨后將兩種軟件所拼接的序列進(jìn)行比對(duì)驗(yàn)證，未發(fā)現(xiàn)有堿基不同。最終使用Pilon 1.3.2(Walkeretal., 2014)對(duì)16 599 bp拼接序列進(jìn)行校正，獲得高橋仁蚧線粒體基因組全長(zhǎng)序列。

1.3 線粒體基因組的注釋與序列分析

利用MITOS2在線服務(wù)器(Berntetal., 2013)對(duì)獲得的線粒體基因組序列進(jìn)行初步注釋。剩余未發(fā)現(xiàn)的tRNA由ARWEN 1.2注釋(Laslett and Canb?ck, 2008)。使用TRF(Tandem Repeats Finder)在線服務(wù)器(https:∥tandem.bu.edu/trf/trf.html)查找串聯(lián)重復(fù)序列(Benson, 1999)。運(yùn)用Phylosuite 1.2.1(Zhangetal., 2020)分析核苷酸密碼子使用頻率和核酸組成。利用公式AT-skew=(A-T)/(A+T)和GC-skew=(G-C)/(G+C)分別計(jì)算AT偏斜和GC偏斜，使用OGDRAM在線可視化軟件繪制高橋仁蚧線粒體全基因組(Greineretal., 2019)。

1.4 系統(tǒng)發(fā)育分析

從GenBank數(shù)據(jù)庫(kù)中下載半翅目15科31種已報(bào)道的線粒體全基因組序列，基于13個(gè)蛋白編碼基因，在Phylosuite軟件中使用MACSE 2.03(Ranwezetal., 2018)進(jìn)行序列多重比對(duì)，并進(jìn)行串聯(lián)，獲得蛋白編碼基因核苷酸序列數(shù)據(jù)集。以蜚蠊目維州散白蟻Reticulitermesvirginicus線粒體全基因組為外群，應(yīng)用IQ-TREE1.6.8(Nguyenetal., 2015)和MrBayes 3.2.6(Ronquistetal., 2012)分別構(gòu)建最大似然法(maximum likelihood, ML)和貝葉斯推斷法(Bayesian inference, BI)系統(tǒng)發(fā)育樹(shù)。IQ-TREE系統(tǒng)樹(shù)分支節(jié)點(diǎn)的置信水平用ultrafast bootstrap approximation approach (UFBoot)(重復(fù)抽樣10 000次)估計(jì)。利用ModelFinder(Kalyaanamoorthyetal., 2017)篩選出最佳核苷酸替換模型為GTR+I+G+F，MrBayes在GTR+I+G+F模型進(jìn)行貝葉斯推斷，運(yùn)算2 000 000代，每1 000代抽樣一次，burn-in參數(shù)設(shè)置為25%。

2 結(jié)果

2.1 高橋仁蚧線粒體基因組結(jié)構(gòu)和堿基組成

高橋仁蚧線粒體基因組序列全長(zhǎng)16 599 bp(GenBank登錄號(hào): MW839575)，由一段控制區(qū)和37個(gè)基因(13個(gè)蛋白編碼基因、22個(gè)tRNA基因和2個(gè)rRNA基因)組成，呈閉合雙鏈環(huán)狀的DNA分子結(jié)構(gòu)(圖1)。通常將多數(shù)編碼基因所在的鏈定義為主要編碼鏈J鏈(majority strand)，則另外一條鏈定義為次要編碼鏈N鏈(minority strand)。線粒體基因組中9個(gè)蛋白編碼基因(cox1,cox2,cox3,cob,atp8,atp6,nad2,nad3和nad6)和14個(gè)tRNA基因(trnM,trnW,trnL2,trnK,trnD,trnG,trnR,trnS1,trnE,trnN,trnA,trnI,trnT和trnS2)在J鏈上，剩余14個(gè)編碼基因在N鏈上。線粒體基因組中A, T, C和G堿基含量分別為47.54%, 36.97%, 11.05%和4.44%?？侫+T含量為84.51%，AT偏斜為0.125，呈現(xiàn)明顯的AT偏向性(表1)。線粒體基因組中共有9處基因重疊(表2)，共76 bp，其中最長(zhǎng)的重疊區(qū)域位于trnE-trnF之間，為37 bp；其次位于trnI-nad2之間，長(zhǎng)為15 bp。存在基因間隔區(qū)24個(gè)，共1 250 bp；最長(zhǎng)間隔區(qū)在trnT-cob處，長(zhǎng)度為184 bp。最短間隔區(qū)位于trnN-trnQ處長(zhǎng)度為2 bp。無(wú)間隔無(wú)重復(fù)區(qū)有4處(表2)。

圖1 高橋仁蚧線粒體基因組結(jié)構(gòu)Fig.1 Structure of the mitochondrial genomeof Aclerda takahashii

表1 高橋仁蚧線粒體基因組核苷酸組成Table 1 Nucleotide composition of the mitochondrial genome of Aclerda takahashi

續(xù)表1 Table 1 continued

表2 高橋仁蚧線粒體基因組結(jié)構(gòu)Table 2 Organization of the mitochondrial genome of Aclerda takahashii

2.2 高橋仁蚧線粒體蛋白質(zhì)編碼基因和相對(duì)同義密碼子使用頻率(RSCU)

高橋仁蚧線粒體基因組中13個(gè)蛋白編碼基因序列總長(zhǎng)為10 632 bp，約占線粒體基因組全序列的64.05%。nad5的序列最長(zhǎng)，為1 620 bp。atp8的序列最短，長(zhǎng)141 bp；均采用標(biāo)準(zhǔn)的起始密碼子(ATN)和終止密碼子(TAA),其中，3個(gè)蛋白編碼基因(cox1,atp8和atp6)的起始密碼子為ATA， 5個(gè)蛋白編碼基因(cox2,cox3,cob,nad1和nad4)的起始密碼子為ATG，其余5個(gè)蛋白編碼基因的起始密碼子為ATT(表2)。共有62個(gè)不同的密碼子，其中有5個(gè)密碼子使用次數(shù)最多：AUA(Met)439次，AUU(Ile)408次，UUU(Phe)390次，UUA(Leu2)366次和AAU(Asn)260次。氨基酸占比由高到低的排序?yàn)镸et(13.62%)，Ile(12.89%)，Phe(12.21%)，Leu2(11.95%)和Asn(9.23%)(圖2)。

圖2 高橋仁蚧線粒體基因組的相對(duì)同義密碼子使用頻率(RSCU)Fig.2 Relative synonymous codon usage (RSCU) in the mitochondrial genome of Aclerda takahashii

2.3 高橋仁蚧線粒體基因組tRNA基因和rRNA基因

22個(gè)tRNA基因中，MITOS2成功注釋了15個(gè)，剩余7個(gè)由ARWEN注釋；長(zhǎng)度變化范圍在48(trnM)～79 bp(trnE)之間；只有10個(gè)tRNA基因(trnR,trnN,trnC,trnE,trnI,trnL1,trnL2,trnK,trnF和trnM)的二級(jí)結(jié)構(gòu)是典型的三葉草結(jié)構(gòu)，其余12個(gè)皆為非典型的結(jié)構(gòu)：trnY和trnW缺失DHU臂，trnS1和trnS2的DHU臂和TΨC臂缺失，剩下8個(gè)tRNA基因都缺失TΨC臂。此外在一些tRNA基因的氨基酸接受臂上出現(xiàn)堿基錯(cuò)配現(xiàn)象，例如trnR(U-U),trnE(U-U),trnN(A-A)以及trnF(U-C)。除了trnN的反密碼子臂上出現(xiàn)1處堿基錯(cuò)配現(xiàn)象(C-A)，其余都能形成完全配對(duì)的臂(圖3)。

圖3 高橋仁蚧線粒體基因組tRNA基因二級(jí)結(jié)構(gòu)Fig.3 Secondary structure of tRNA genes in the mitochondrial genome of Aclerda takahashii

2.4 系統(tǒng)發(fā)育

運(yùn)用最大似然法和貝葉斯推斷法構(gòu)建的系統(tǒng)發(fā)育樹(shù)均展示了相同的拓?fù)浣Y(jié)構(gòu)(圖4)，并且每個(gè)分支節(jié)點(diǎn)上都有高的支持率(自舉值>65，后驗(yàn)概率>0.90)。二者結(jié)果充分展示了這些總科的如下系統(tǒng)發(fā)育關(guān)系：(((跳蝽總科(Saldoidea)+細(xì)蝽總科(Leptopodoidea))+(臭蟲總科(Cimicoidea)+獵蝽總科(Reduvioidea)))+(角蟬總科(Membracoidea)+(蠟蟬總科(Fulgoroidea)+(木虱總科(Psylloidea)+(粉虱總科(Aleyrodoidea)+(蚜總科(Aphidoidea)+蚧總科(Coccoidea)))))))。兩個(gè)系統(tǒng)發(fā)育樹(shù)展示蚧科的日本龜蠟蚧C.japonicus、咖啡黑盔蚧S.coffeae、朝鮮球堅(jiān)蚧Didesmococcuskoreanus均分布在一個(gè)分支上，高橋仁蚧作為仁蚧科代表，單獨(dú)形成一個(gè)分支。

圖4 最大似然法構(gòu)建的基于線粒體基因組的13個(gè)蛋白質(zhì)編碼基因核苷酸序列的半翅目系統(tǒng)發(fā)育樹(shù)Fig.4 Phylogenetic tree of Hemiptera based on the sequences of 13 protein-coding genesof mitochondrial genome using maximum likelihood method節(jié)點(diǎn)旁左側(cè)數(shù)字代表Bootstrap支持率，右側(cè)數(shù)字代表貝葉斯后驗(yàn)概率。Numbers on nodes refer to bootstrap support values (left) and Bayesian posterior probabilities (right).

3 討論

高橋仁蚧tRNA二級(jí)結(jié)構(gòu)具有特殊性，一部分tRNA無(wú)法通過(guò)MITOS2得到注釋(表2)。主要是由于這些tRNA是非典型的三葉草結(jié)構(gòu)，缺失DHU臂或TΨC臂，其中tRNAser(S1)和tRNAser(S2)更是同時(shí)缺失DHU臂和TΨC臂(圖3)。同時(shí)在部分tRNA的氨基酸臂上，有個(gè)別堿基不配對(duì)，這些原因?qū)е买幌x的tRNA無(wú)法通過(guò)線粒體注釋軟件MITOS2準(zhǔn)確注釋。一般情況下，動(dòng)物的trnS1(AGN)通常會(huì)缺失DHU臂，而不具有經(jīng)典的三葉草結(jié)構(gòu)(Boore, 1999)。這種trnS1缺失DHU臂的現(xiàn)象在半翅目中均有發(fā)現(xiàn)(郭仲龍和袁明龍, 2016)，白背飛虱Sogatellafurcifera的線粒體基因組的2個(gè)trnS基因和1個(gè)trnG基因也有缺臂的現(xiàn)象存在(Zhangetal., 2014)。在高橋仁蚧基因組中存在嚴(yán)重的tRNA缺臂現(xiàn)象，22個(gè)tRNA中的12個(gè)均存在缺臂現(xiàn)象(圖3)，這一現(xiàn)象在蚧科的線粒體基因組中也有發(fā)現(xiàn)(Luetal., 2020)。然而在其他類群的昆蟲中，大量tRNA的缺臂現(xiàn)象鮮有發(fā)生(Luetal., 2020)。類似的報(bào)道主要來(lái)自于線蟲(Jühlingetal., 2012; Lorenzetal., 2017)、螨類(Xueetal., 2016; 趙亞男和李朝品, 2020)和蛛形綱動(dòng)物(Masta and Boore, 2008)。在后生動(dòng)物中，線粒體tRNA通過(guò)編譯可以修復(fù)這些缺失的臂，形成典型的三葉草結(jié)構(gòu)的tRNA(Lavrovetal., 2000; Segoviaetal., 2011)。這種轉(zhuǎn)錄后的編譯可能廣泛存在于仁蚧科的類群中。同時(shí)，我們也發(fā)現(xiàn)高橋仁蚧線粒體蛋白編碼基因具有很強(qiáng)的密碼子偏好性(圖2)，這與半翅目其他昆蟲的密碼子使用情況類似(Luetal., 2020)。

半翅目昆蟲的AT含量一般在70%～80%之間，最低AT含量的半翅目線粒基因組為粉虱科Bemisiaafer(65.67%)，最高的是粉虱科Aleurodicusdugesii(86.33%)(Wangetal., 2015)。其中胸喙亞目(Sternorrhyncha)是半翅目4個(gè)亞目中AT含量最高的亞目，平均AT含量達(dá)到80.97%(郭仲龍和袁明龍, 2016)。目前GenBank中最高AT含量(88.02%)的昆蟲線粒體基因組為膜翅目錘角細(xì)蜂科Trichopriadrosophilae線粒體基因組(GenBank登錄號(hào): NC_048491)。作為胸喙亞目的一員，3個(gè)已知的蚧蟲線粒體基因組AT含量均超過(guò)了80.00%(Dengetal., 2019; Luetal., 2020)，高橋仁蚧的AT含量達(dá)到84.51%(表1)。超高的AT含量造成測(cè)序上的困難，AT重復(fù)片段以及poly (A)的特殊結(jié)構(gòu)經(jīng)常導(dǎo)致一代測(cè)序的質(zhì)量不佳，出現(xiàn)套峰或測(cè)序信號(hào)的衰減中斷。在高橋仁蚧的線粒體基因組中，一個(gè)非常極端的例子是在nad5基因的位置上，連續(xù)出現(xiàn)了31個(gè)A堿基，這或可以解釋為何沒(méi)有一個(gè)蚧蟲線粒體基因組是通過(guò)一代測(cè)序獲得的原因。二代測(cè)序技術(shù)，也就是高通量測(cè)序技術(shù)，本身由于對(duì)于GC有偏好性，導(dǎo)致在AT含量高的地方往往測(cè)序覆蓋度不一致，甚至缺失覆蓋(Chenetal., 2013; Browneetal., 2020)。這導(dǎo)致利用二代測(cè)序技術(shù)拼接蚧蟲線粒體基因組時(shí)，獲得的線粒體基因組常碎片化。再加上由于蚧蟲本身個(gè)體小，采集困難，易被寄生蜂寄生(Dengetal., 2012)，這些原因都導(dǎo)致很難獲得完整蚧蟲基因組，也成為GenBank數(shù)據(jù)庫(kù)中嚴(yán)重缺少蚧蟲線粒體基因組數(shù)據(jù)的主要原因。

在構(gòu)建的半翅目代表物種系統(tǒng)發(fā)育樹(shù)上(圖4)，仁蚧科與蚧科物種親緣關(guān)系最近，這與之前基于18S rRNA, 28S rRNA和EF-1α的3個(gè)核基因片段的蚧總科系統(tǒng)發(fā)育研究的結(jié)果(Vea and Grimaldi, 2016)一致。由于蚧蟲線粒體基因組數(shù)據(jù)有限，在蚧蟲中基于線粒體多基因的系統(tǒng)發(fā)育研究還有待開(kāi)展。已知的蚧蟲線粒體基因組構(gòu)建的系統(tǒng)發(fā)育樹(shù)也無(wú)法真實(shí)反映蚧蟲的進(jìn)化的格局和系統(tǒng)發(fā)育關(guān)系。半翅目蚧總科中線粒體基因組數(shù)據(jù)的缺乏，嚴(yán)重阻礙了對(duì)蚧蟲本身以及整個(gè)半翅目昆蟲利用線粒體基因組開(kāi)展進(jìn)化、系統(tǒng)發(fā)育和生物地理等問(wèn)題的研究，獲得更多不同科的蚧蟲線粒體基因組序列將成為解決這一問(wèn)題的關(guān)鍵。

本研究獲得了高橋仁蚧線粒體基因組全序列，是仁蚧科首個(gè)線粒體基因組序列，提供了仁蚧科線粒體基因結(jié)構(gòu)和組成的基礎(chǔ)數(shù)據(jù)，為后續(xù)獲得仁蚧科更多物種的線粒體基因組提供了參考。