全蝽線粒體全基因組序列及其進(jìn)化分析

高曉云，魏久鋒，丁曉飛，劉迎香，趙清

（山西農(nóng)業(yè)大學(xué) 植物保護(hù)學(xué)院，山西 太谷 030801）

昆蟲(chóng)線粒體基因組由一個(gè)長(zhǎng)度為15～20 kb的閉合環(huán)狀雙鏈DNA分子組成[1]，通常編碼37個(gè)基因（13個(gè)蛋白質(zhì)編碼基因（PCGs）、22個(gè)轉(zhuǎn)運(yùn)RNA基因（tRNAs）、2個(gè)核糖體RNA基因（rRNAs））和1個(gè)控制區(qū)域[2]。與核基因相比，線粒體基因組小且多態(tài)，進(jìn)化速度快，很少進(jìn)行基因重組[3]。因此，許多研究者用它鑒定物種并進(jìn)行系統(tǒng)發(fā)育關(guān)系重建[4-6]。在過(guò)去的10 a中，隨著測(cè)序技術(shù)的高速發(fā)展，線粒體基因組測(cè)序量顯著增加，被廣泛應(yīng)用于系統(tǒng)地理學(xué)、分子進(jìn)化和系統(tǒng)發(fā)育領(lǐng)域[7-8]。

蝽總科是異翅亞目中最大的總科之一，廣泛分布于世界各地，大約包含8 000多個(gè)物種[9]。近年來(lái)，許多學(xué)者研究了蝽總科的系統(tǒng)發(fā)育關(guān)系，試圖為蝽總科的分類提供更多的見(jiàn)解。例如，Gapud首次通過(guò)系統(tǒng)發(fā)育方法研究了蝽總科，并利用形態(tài)學(xué)數(shù)據(jù)推斷了各科間的系統(tǒng)發(fā)育關(guān)系[9]。在不同的分析方法下，蝽總科被認(rèn)為是單系的[7,10-12]。GRAZIA等學(xué)者根據(jù)分子數(shù)據(jù)集[11-14]或形態(tài)學(xué)[9,14]認(rèn)為，異蝽科是蝽總科中其他科的姐妹群[14]。在隨后的研究中，其他人使用部分?jǐn)?shù)據(jù)來(lái)推斷兜蝽科[15]和荔蝽科之間的關(guān)系[15-16]。后來(lái)，越來(lái)越多學(xué)者經(jīng)常發(fā)現(xiàn)兜蝽科和荔蝽科之間的姐妹關(guān)系[11-13,17-19]。研究發(fā)現(xiàn)，蝽總科內(nèi)除了兜蝽科和荔蝽科之外，其他科級(jí)關(guān)系混亂復(fù)雜，尤其是蝽科的位置[19-20]。比如，1991年Gapud基于形態(tài)構(gòu)建的蝽總科系統(tǒng)發(fā)育樹(shù)中蝽科與土蝽科是姐妹關(guān)系；2008年GRZAIA等[14]通過(guò)形態(tài)特征和COI、12S、16S、18S、28S相結(jié)合證實(shí)了部分Phloeidae和蝽科是姐妹群，隨后僅基于形態(tài)學(xué)卻發(fā)現(xiàn)Lestoniidae和蝽科的姐妹關(guān)系；2010年KMENT等[16]指出，蝽科和龜蝽科是姐妹群；2017年LIS等[17]基于18S和28S，證實(shí)了蝽科和（兜蝽科+荔蝽科）是姐妹關(guān)系；2019年有學(xué)者指出，蝽科和盾蝽科是姐妹組[7]；2021年XU等[18]基于線粒體基因組構(gòu)建的系統(tǒng)發(fā)育關(guān)系發(fā)現(xiàn)，蝽科和（異蝽科+同蝽科）是姐妹關(guān)系。

本研究首次對(duì)全蝽的完整線粒體基因組進(jìn)行了測(cè)序、組裝及注釋分析，同時(shí)對(duì)線粒體基因組和其他蝽總科物種進(jìn)行了核酸組成、密碼子使用、基因順序以及tRNA二級(jí)結(jié)構(gòu)的比較分析，研究結(jié)果可進(jìn)一步豐富蝽科線粒體基因組學(xué)數(shù)據(jù)，通過(guò)全蝽線粒體基因組序列并聯(lián)合其他蝽總科物種的線粒體基因組數(shù)據(jù)，確定全蝽在蝽科系統(tǒng)發(fā)育關(guān)系中的位置，并討論分析蝽科在蝽總科中的分類地位。

1 材料和方法

1.1 標(biāo)本采集與測(cè)序

本試驗(yàn)所用的全蝽標(biāo)本于2019年8月采自安徽省安慶市岳西縣坪鷂落村，并在提取DNA之前儲(chǔ)存在100%乙醇中；證據(jù)標(biāo)本（編號(hào)SXAU2-6）保存在山西農(nóng)業(yè)大學(xué)昆蟲(chóng)學(xué)研究所。

取全蝽的胸部和腿部肌肉組織，將其研磨成粉末狀后嚴(yán)格按照基因組DNA提取試劑盒中使用手冊(cè)（德國(guó)Hilden Qiagen）進(jìn)行基因組的提取，使用分光光度計(jì)檢驗(yàn)合格后，送至百邁克生物技術(shù)公司（中國(guó)上海）測(cè)序。全蝽的線粒體全基因組上傳至GenBank（登錄號(hào)：OL884612）。

1.2 序列的拼接、注釋與分析

利用Geneious R9軟件進(jìn)行序列的組裝，并對(duì)組裝后的完整線粒體基因組序列進(jìn)行注釋。13個(gè)蛋白質(zhì)編碼基因由開(kāi)放閱讀框查找工具ORF依照無(wú)脊椎動(dòng)物線粒體遺傳密碼子進(jìn)行查找，并通過(guò)與其他昆蟲(chóng)同源基因的比對(duì)來(lái)確定每個(gè)蛋白質(zhì)基因的起始和終止密碼子的位置。利用MITOS Web（http://mitos.bioinf.uni-leipzig.de/）預(yù)測(cè)22個(gè)tRNAs的位置和二級(jí)結(jié)構(gòu)。2個(gè)rRNAs基因的邊界通過(guò)與其他已知的蝽科rRNA基因的比較來(lái)確定。

蛋白編碼基因的密碼子使用情況和核苷酸組成通 過(guò)Mega v6.0進(jìn)行了 統(tǒng)計(jì)[21]，AT和GC偏斜 通過(guò)公式（1）和（2）進(jìn)行計(jì)算。同時(shí)，利用DnaSP v5.0計(jì)算每個(gè)非同義位點(diǎn)的非同義替換數(shù)（Ka）、每個(gè)同義位點(diǎn)的同義替換數(shù)（Ks）。

1.3 系統(tǒng)發(fā)育樹(shù)的構(gòu)建

基于13個(gè)蛋白編碼基因PartitionFinder得到的分區(qū)和模型如表1所示。

表1 基于13個(gè)蛋白編碼基因PartitionFinder得到的分區(qū)和模型Tab.1 Partitions and models based on Partition Finder of 13 protein-coding genes

從數(shù)據(jù)庫(kù)下載9科27種的線粒體全基因組序列，分別基于13個(gè)蛋白編碼基因和13個(gè)基因密碼子前2位，進(jìn)行多重比對(duì)，利用PartitionFinder v1.1.1確定劃分?jǐn)?shù)據(jù)的最佳分區(qū)及替換模型。系統(tǒng)發(fā)育樹(shù)法的構(gòu)建選取最大似然法（ML）和貝葉斯法（BI）。ML樹(shù)在RAxML-7.0.3中生成，通過(guò)dos環(huán)境輸入命令，命令中需指定外群名稱和分區(qū)文件，1 000次快速自舉檢驗(yàn)分析。BI樹(shù)在MrBayes 3.2.2中運(yùn)行計(jì)算，將PartitionFinder中的分區(qū)信息（表1和表2）加到建樹(shù)命令末尾，以隨機(jī)起始樹(shù)開(kāi)始運(yùn)算。參數(shù)設(shè)置為：ngenerations=2 000 000，samplefreq=1 000，printfreq=1 000。當(dāng)收斂值低于0.01時(shí)，程序?qū)⑼Ｖ惯\(yùn)行。

表2 基于13個(gè)蛋白編碼基因的第一位和第二位密碼子PartitionFinder得到的分區(qū)和模型Tab.2 Partitions and models based on Partition Finder of first codon and second codon of 13 protein-coding genes

2 結(jié)果與分析

2.1 線粒體基因組結(jié)構(gòu)

全蝽的線粒體基因組序列是一個(gè)15 479 bp的環(huán)狀DNA分子（圖1）。這個(gè)線粒體基因組有37個(gè)基 因，包 括2個(gè)rRNA、22個(gè)tRNA和13個(gè)PCG（Protein Coding Gene），其基因順序與亞庫(kù)巴果蠅Drosophila yakuba一致（表3）。37個(gè)基因中有23個(gè)基 因 由J鏈 編 碼（trnI、trnM、nad2、trnW、cox1、trnL2、cox2、trnK、trnD、atp8、atp6、cox3、trnG、nad3、trnA、trnR、trnN、trnS1、trnE、trnT、nad6、cob、trnS2），14個(gè) 基 因 由N鏈 編 碼（trnQ、trnC、trnY、trnF、nad5、trnH、nad4、nad4L、trnP、nad1、trnL1、rrnL、trnV、rrnS）。全蝽的線粒體基因組結(jié)構(gòu)緊湊，基因重疊長(zhǎng)度范圍為1～8 bp，有7處重疊。其中atp8/atp6和nad4L/nad4這2對(duì)基因的重疊均為ATGATAA 7個(gè)核苷酸。rrnS和trnI之間的非編碼區(qū)長(zhǎng)度為750 bp，在trnA和trnR之間也有一個(gè)較大的非編碼區(qū)。

圖1 全蝽的線粒體基因組環(huán)狀結(jié)構(gòu)Fig.1 Circular structure of the mitochondrial gemome of Homalogonia obtusa

表3 全蝽線粒體基因組結(jié)構(gòu)Tab.3 Structure of mitochondrial genome of Homalogonia obtusa

續(xù)表3全蝽線粒體基因組結(jié)構(gòu)Tab.3（Continued）Structure of mitochondrial genome of Homalogonia obtusa

2.2 堿基組成與蛋白編碼基因

從表4可以看出，全蝽線粒體基因組的AT含量為76.18%，AT-skew為0.108，具有明顯的AT偏好性，位于編碼N鏈蛋白質(zhì)編碼基因中T偏斜值較大,而J鏈的T偏斜值較小。PCGs、tRNA、rRNA的AT含量均值分別為75.40%、77.70%、78.70%。

表4 全蝽線粒體基因核酸組成情況Tab.4 Nucleic acid composition of mitochondrial gene of Homalogonia obtusa

13個(gè)蛋白編碼基因中除cox1基因AT含量為69.08%外，其余12個(gè)蛋白編碼基因中的AT含量均高于70%，且偏向性最大的均出現(xiàn)在密碼子第2位。全蝽的13個(gè)蛋白編碼基因中AT含量最高的3個(gè)基因依次是atp8、nad6和nad4L；AT含量比較低的蛋白質(zhì)基因是cox1、cox2和cox3，故推測(cè)高AT含量的基因進(jìn)化速率快于低AT含量的基因，即atp8的進(jìn)化速率最高，而cox1的進(jìn)化速率最低。

昆蟲(chóng)mtDNA的高A+T含量在很大程度上影響了蛋白密碼子的使用，本研究中對(duì)同義密碼子使用情況進(jìn)行了分析。全蝽中主要是以T或A為第3位的同義密碼子的出現(xiàn)頻率遠(yuǎn)高于其他密碼子。全蝽的13個(gè)PCGs中，A+T含量高于C+G含量，密碼子的使用也反映了核酸組成的偏斜。由同義密碼子使用度發(fā)現(xiàn)，最常用的密碼子是NNA，最不常用的密碼子是NNG（圖2）。在全蝽中最常用的4個(gè)氨基酸密碼子是UUA（亮氨酸），其次是AUU（異亮氨酸），再次是UUU（苯丙氨酸），最后是AUA（甲硫氨酸）。

圖2 全蝽的相對(duì)密碼子使用率Fig.2 The relative synonymous codon usage（RSCU）in Homalogonia obtuse

在全蝽的起始密碼子和終止密碼子中除cox1和nad1以TTG啟動(dòng)外，其他基因起始密碼子均為ATA和ATG（ATN）。nad2、atp8、atp6、cox3、nad4、nad5、nad6、nad4l、nad1、cob等10個(gè)蛋白質(zhì)基因的終止密碼子是TAA，另外3個(gè)cox1、cox2、nad3以單個(gè)T終止。起始密碼子和終止密碼子也主要由A和T組成，這與蝽科其他昆蟲(chóng)的情況一致[11]。

對(duì)于單個(gè)蛋白質(zhì)編碼基因，本研究計(jì)算了它們的同義位點(diǎn)的非同義替代率（Ka）和同義替代率（Ks）（圖3）。在核苷酸水平上coxl、cytb、cox3和cox2這4個(gè)基因最為保守，進(jìn)化最慢。非同義與同義置換率的比率（ω=Ka/Ks）在13個(gè)基因間變異顯 著，coxl、cytb、cox3和cox2偏 低（0.13～0.26），nad1、atp6和nad3中等（0.34～0.35），nad4L、nad5、nad4、nad6、nad2和atp8最高（0.51～0.83）。這與通常認(rèn)為atp8基因是線粒體基因組中進(jìn)化最快的觀點(diǎn)一致，表明這些基因承受不同的選擇性功能約束。然而所有蛋白基因的ω值均小于1，表明這些基因的進(jìn)化受到純化選擇作用。

圖3 蝽總科的13個(gè)蛋白質(zhì)編碼的進(jìn)化速率Fig.3 Evolutionary rates of 13 protein-coding genes in Pentatomoidea

2.3 tRNAs基因 與rRNAs基因

全蝽線粒體基因組中共有22個(gè)tRNAs基因，可攜帶20種氨基酸，其中trnL存在2個(gè)不同的tRNA受 體TAA和TAG，trnS存 在2個(gè) 不 同 的tRNA受 體GCT和TGA。22個(gè)tRNAs基因大多數(shù)可以折疊成三葉草狀的二級(jí)結(jié)構(gòu)（圖4）。與大多數(shù)報(bào)道的蝽科昆蟲(chóng)線粒體基因組一樣，trnS1（GCU）的DHU臂被一個(gè)簡(jiǎn)單的環(huán)取代，特殊的是，在全蝽中trnA沒(méi)有TΨC環(huán)和TΨC臂。同時(shí)，在繪制全蝽的22個(gè)tRNAs的二級(jí)結(jié)構(gòu)時(shí)發(fā)現(xiàn)，tRNA臂結(jié)構(gòu)在維持tRNA二級(jí)結(jié)構(gòu)穩(wěn)定中起著關(guān)鍵性的作用，因?yàn)樗鼈兊暮塑账嵝蛄懈鼮楸Ｊ?。tRNA臂連接處以及額外環(huán)最為多變，二者都為進(jìn)化速度比較快的區(qū)域。

圖4 全蝽的22個(gè)tRNA的二級(jí)結(jié)構(gòu)Fig.4 Secondary structure of 22 tRNA of Homalogonia obtusa

rrnL和rrnS的長(zhǎng)度分別為1 278、816 bp；其中，rrnL位于trnV和trnL（TAG）間隙，rrnS位于trnV和控制區(qū)間隙。在rrnL的二級(jí)結(jié)構(gòu)中有6個(gè)結(jié)構(gòu)域，其中第3結(jié)構(gòu)域缺失，rrnS的二級(jí)結(jié)構(gòu)中有3個(gè)結(jié)構(gòu)域。

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

本研究以2條長(zhǎng)蝽科的物種序列作為外群、蝽總科26個(gè)的物種作為內(nèi)群，基于貝葉斯推斷BI（圖5、6）和最大似然法ML（圖7、8），進(jìn)行蝽總科內(nèi)的系統(tǒng)發(fā)育分析。選擇13個(gè)PCGs和13個(gè)PCGs的第1位和第2位密碼子（PCG12）串聯(lián)序列數(shù)據(jù)集構(gòu)建系統(tǒng)發(fā)育樹(shù)。

基于13個(gè)PCGs在貝葉斯分析中，蝽總科內(nèi)的關(guān)系是（異蝽科+（蝽科+（同蝽科+（（龜蝽科+盾蝽科）+（土蝽科+（兜蝽科+荔蝽科）））））），這與基于PCG12在貝葉斯分析得到的拓?fù)浣Y(jié)構(gòu)一致；在最大似然法分析中，基于13PCGs的蝽總科內(nèi)的關(guān)系是（異蝽科+（同蝽科+（蝽科+（（龜蝽科+盾蝽科）+（土蝽科+（兜蝽科+荔蝽科））））））；基于PCG12的蝽總科內(nèi)的關(guān)系卻是（蝽科+（（異蝽科+同蝽科）+（（龜蝽科+盾蝽科）+（土蝽科+（兜蝽科+荔蝽科））））））。

4種拓?fù)浣Y(jié)構(gòu)的結(jié)果均支持全蝽隸屬于蝽科，且對(duì)于盾蝽科和龜蝽科之間的姐妹群關(guān)系、兜蝽科和荔蝽科的姐妹群關(guān)系均具有很高的支持率。蝽科的單系性也得到很好的支持。4種拓?fù)浣Y(jié)構(gòu)的不同之處主要是蝽科與同蝽科和異蝽科之間的系統(tǒng)發(fā)育關(guān)系。貝葉斯推斷結(jié)果顯示（圖5、6），蝽科進(jìn)化比同蝽科早，但是最大似然法分析得到的結(jié)果卻是同蝽科是繼異蝽科后進(jìn)化最早的科（圖7），基于PCG12的分析中（圖8），蝽科比異蝽科進(jìn)化還早。

圖5 基于13個(gè)蛋白編碼基因利用貝葉斯分析構(gòu)建的蝽總科系統(tǒng)發(fā)育樹(shù)Fig.5 Phylogenetic trees of Pentatomoidea constructed based on 13 protein-coding genes using Bayesian inference

圖7 基于13個(gè)蛋白編碼基因利用最大似然法構(gòu)建的蝽總科系統(tǒng)發(fā)育樹(shù)Fig.7 Phylogenetic trees of Pentatomoidea constructed based on 13 protein-coding genes using RaxML inference

圖8 基于13個(gè)蛋白編碼基因的第1位和第2位密碼子利用最大似然法構(gòu)建的蝽總科系統(tǒng)發(fā)育樹(shù)Fig.8 Phylogenetic trees of Pentatomoidea constructed based on first codon and second codon of 13 protein-coding genes using RaxML inference

圖6 基于13個(gè)蛋白編碼基因的第1位和第2位密碼子利用貝葉斯分析構(gòu)建的蝽總科系統(tǒng)發(fā)育樹(shù)Fig.6 Phylogenetic trees of Pentatomoidea constructed based on first codon and second codon of 13 protein-coding genes using Bayesian inference

3 結(jié)論與討論

本研究首次對(duì)全蝽的線粒體基因組進(jìn)行了完整測(cè)序和報(bào)道，線粒體基因組序列的結(jié)構(gòu)和順序緊湊，與其他蝽科的物種的結(jié)構(gòu)和順序相同[12,19,22]。通過(guò)對(duì)全蝽線粒體基因組13個(gè)蛋白編碼基因的同義密碼子使用度、起始密碼子和終止密碼子及A+T偏斜性等方面分析可以看出，堿基的A+T含量偏斜與進(jìn)化速率相關(guān)。在一些進(jìn)化較慢的基因中A+T含量比 較 低，比如cox1、cox2和cox3在大多數(shù)物種中的A+T含量普遍偏低，其中cox1尤為明顯；而在一些進(jìn)化較快的基因中A+T含量都比較高，如atp8、nad4L和nad6在所選物種中的A+T含量基本偏高。在堿基偏斜方面，N鏈蛋白基因都是T偏斜度很高，而J鏈蛋白基因基本為T偏斜但是偏斜度較小。

ATN是13個(gè)蛋白編碼基因中最常見(jiàn)的起始密碼子，此外，TTG也存在于大多數(shù)昆蟲(chóng)的線粒體基因組中。在全蝽中，除了cox1和nad1之外，所有蛋白編碼基因的起始密碼子都是ATN，以TTG開(kāi)始。全蝽的10個(gè)蛋白編碼基因的終止密碼子是TAA，其余3個(gè)蛋白編碼基因（cox1、cox2和nad3）的終止密碼子是單個(gè)T。大多數(shù)蛋白編碼基因的終止密碼子是TAN，然而，一些基因的終止密碼子也是單個(gè)T或TA[23]。根據(jù)研究發(fā)現(xiàn)，不完全終止密碼子T或TA在轉(zhuǎn)錄過(guò)程中進(jìn)行多聚腺苷酸化，最終成為TAA[24]。

本研究中，BI和ML系統(tǒng)發(fā)育樹(shù)的結(jié)果均顯示，全蝽位于蝽科支系，在分子證據(jù)上說(shuō)明全蝽隸屬于蝽科。同時(shí)，結(jié)果還很好地證明了蝽科的單系性，揭示了蝽科和（龜蝽科+盾蝽科）是姐妹群，這與以前的研究結(jié)果一致[18]。文中構(gòu)建的4棵系統(tǒng)發(fā)育樹(shù)結(jié)果顯示，蝽科內(nèi)部的短喙蝽亞科和益蝽亞科是單系性，但是舌盾蝽亞科和蝽亞科是并系，而且蝽亞科內(nèi)部的物種關(guān)系更為混亂。1833年，YANG[25]基于盾片形狀將Graphosoma歸于蝽亞科Graphosomini，但是其他的研究中都是把它置于舌盾蝽亞科[26]。本研究中，也將Graphosoma歸于舌盾蝽亞科，但是它并不與舌盾蝽亞科的Scotinophara聚為一支。造成這樣的研究結(jié)果可能是因?yàn)轵砜莆锓N的形態(tài)復(fù)雜，變異較大，特征不穩(wěn)定，對(duì)其進(jìn)行系統(tǒng)的劃分比較困難。學(xué)者劃分蝽科各亞科的標(biāo)準(zhǔn)也僅僅是根據(jù)此類昆蟲(chóng)的外部形態(tài)的某一個(gè)或某幾個(gè)特征，這并不足以確定系統(tǒng)進(jìn)化關(guān)系；也可能是因?yàn)楸狙芯恐羞x取的線粒體基因組基因含量不如核基因多，所以，不能充分表達(dá)各亞科之間的進(jìn)化遺傳信息，因而劃分的結(jié)果會(huì)有分歧；當(dāng)然本研究中也受到了樣本不足的限制。所以，要想解決蝽亞科和舌盾蝽亞科的分類關(guān)系存在的問(wèn)題，還需要更多的形態(tài)和分子數(shù)據(jù)構(gòu)建蝽科內(nèi)各個(gè)亞科的系統(tǒng)發(fā)育位置。

本研究通過(guò)對(duì)全蝽線粒體的全基因組測(cè)序和生物信息學(xué)分析，發(fā)現(xiàn)全椿線粒體基因組全長(zhǎng)15 479 bp，AT含量為76.18%。其中，13個(gè)編碼基因的起始密碼子多為ATN和TTG，終止密碼子為TAN或單個(gè)T堿基，atp8進(jìn)化最快，cox1進(jìn) 化 最慢；22個(gè)tRNA中trnA沒(méi)有TΨC環(huán)和TΨC臂，trnS1（GCU）沒(méi)有DHU臂，其他20個(gè)tRNA為典型的三葉草結(jié)構(gòu)。系統(tǒng)發(fā)育樹(shù)結(jié)果顯示，全蝽位于蝽科支系且蝽科為單系。本研究結(jié)果豐富了蝽科線粒體基因組學(xué)和NCBI數(shù)據(jù)庫(kù)的基本數(shù)據(jù)，為進(jìn)一步深化蝽科和蝽總科的系統(tǒng)發(fā)育關(guān)系討論提供了參考。