全蝽線粒體全基因組序列及其進化分析

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

（山西農業大學 植物保護學院，山西 太谷 030801）

昆蟲線粒體基因組由一個長度為15～20 kb的閉合環狀雙鏈DNA分子組成[1]，通常編碼37個基因（13個蛋白質編碼基因（PCGs）、22個轉運RNA基因（tRNAs）、2個核糖體RNA基因（rRNAs））和1個控制區域[2]。與核基因相比，線粒體基因組小且多態，進化速度快，很少進行基因重組[3]。因此，許多研究者用它鑒定物種并進行系統發育關系重建[4-6]。在過去的10 a中，隨著測序技術的高速發展，線粒體基因組測序量顯著增加，被廣泛應用于系統地理學、分子進化和系統發育領域[7-8]。

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

本研究首次對全蝽的完整線粒體基因組進行了測序、組裝及注釋分析，同時對線粒體基因組和其他蝽總科物種進行了核酸組成、密碼子使用、基因順序以及tRNA二級結構的比較分析，研究結果可進一步豐富蝽科線粒體基因組學數據，通過全蝽線粒體基因組序列并聯合其他蝽總科物種的線粒體基因組數據，確定全蝽在蝽科系統發育關系中的位置，并討論分析蝽科在蝽總科中的分類地位。

1 材料和方法

1.1 標本采集與測序

本試驗所用的全蝽標本于2019年8月采自安徽省安慶市岳西縣坪鷂落村，并在提取DNA之前儲存在100%乙醇中；證據標本（編號SXAU2-6）保存在山西農業大學昆蟲學研究所。

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

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

利用Geneious R9軟件進行序列的組裝，并對組裝后的完整線粒體基因組序列進行注釋。13個蛋白質編碼基因由開放閱讀框查找工具ORF依照無脊椎動物線粒體遺傳密碼子進行查找，并通過與其他昆蟲同源基因的比對來確定每個蛋白質基因的起始和終止密碼子的位置。利用MITOS Web（http://mitos.bioinf.uni-leipzig.de/）預測22個tRNAs的位置和二級結構。2個rRNAs基因的邊界通過與其他已知的蝽科rRNA基因的比較來確定。

蛋白編碼基因的密碼子使用情況和核苷酸組成通 過Mega v6.0進行了 統計[21]，AT和GC偏斜 通過公式（1）和（2）進行計算。同時，利用DnaSP v5.0計算每個非同義位點的非同義替換數（Ka）、每個同義位點的同義替換數（Ks）。

1.3 系統發育樹的構建

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

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

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

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

2 結果與分析

2.1 線粒體基因組結構

全蝽的線粒體基因組序列是一個15 479 bp的環狀DNA分子（圖1）。這個線粒體基因組有37個基 因，包 括2個rRNA、22個tRNA和13個PCG（Protein Coding Gene），其基因順序與亞庫巴果蠅Drosophila yakuba一致（表3）。37個基因中有23個基 因 由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個 基 因 由N鏈 編 碼（trnQ、trnC、trnY、trnF、nad5、trnH、nad4、nad4L、trnP、nad1、trnL1、rrnL、trnV、rrnS）。全蝽的線粒體基因組結構緊湊，基因重疊長度范圍為1～8 bp，有7處重疊。其中atp8/atp6和nad4L/nad4這2對基因的重疊均為ATGATAA 7個核苷酸。rrnS和trnI之間的非編碼區長度為750 bp，在trnA和trnR之間也有一個較大的非編碼區。

圖1 全蝽的線粒體基因組環狀結構Fig.1 Circular structure of the mitochondrial gemome of Homalogonia obtusa

表3 全蝽線粒體基因組結構Tab.3 Structure of mitochondrial genome of Homalogonia obtusa

續表3全蝽線粒體基因組結構Tab.3（Continued）Structure of mitochondrial genome of Homalogonia obtusa

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

從表4可以看出，全蝽線粒體基因組的AT含量為76.18%，AT-skew為0.108，具有明顯的AT偏好性，位于編碼N鏈蛋白質編碼基因中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個蛋白編碼基因中除cox1基因AT含量為69.08%外，其余12個蛋白編碼基因中的AT含量均高于70%，且偏向性最大的均出現在密碼子第2位。全蝽的13個蛋白編碼基因中AT含量最高的3個基因依次是atp8、nad6和nad4L；AT含量比較低的蛋白質基因是cox1、cox2和cox3，故推測高AT含量的基因進化速率快于低AT含量的基因，即atp8的進化速率最高，而cox1的進化速率最低。

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

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

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

對于單個蛋白質編碼基因，本研究計算了它們的同義位點的非同義替代率（Ka）和同義替代率（Ks）（圖3）。在核苷酸水平上coxl、cytb、cox3和cox2這4個基因最為保守，進化最慢。非同義與同義置換率的比率（ω=Ka/Ks）在13個基因間變異顯 著，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）。這與通常認為atp8基因是線粒體基因組中進化最快的觀點一致，表明這些基因承受不同的選擇性功能約束。然而所有蛋白基因的ω值均小于1，表明這些基因的進化受到純化選擇作用。

圖3 蝽總科的13個蛋白質編碼的進化速率Fig.3 Evolutionary rates of 13 protein-coding genes in Pentatomoidea

2.3 tRNAs基因 與rRNAs基因

全蝽線粒體基因組中共有22個tRNAs基因，可攜帶20種氨基酸，其中trnL存在2個不同的tRNA受 體TAA和TAG，trnS存 在2個 不 同 的tRNA受 體GCT和TGA。22個tRNAs基因大多數可以折疊成三葉草狀的二級結構（圖4）。與大多數報道的蝽科昆蟲線粒體基因組一樣，trnS1（GCU）的DHU臂被一個簡單的環取代，特殊的是，在全蝽中trnA沒有TΨC環和TΨC臂。同時，在繪制全蝽的22個tRNAs的二級結構時發現，tRNA臂結構在維持tRNA二級結構穩定中起著關鍵性的作用，因為它們的核苷酸序列更為保守。tRNA臂連接處以及額外環最為多變，二者都為進化速度比較快的區域。

圖4 全蝽的22個tRNA的二級結構Fig.4 Secondary structure of 22 tRNA of Homalogonia obtusa

rrnL和rrnS的長度分別為1 278、816 bp；其中，rrnL位于trnV和trnL（TAG）間隙，rrnS位于trnV和控制區間隙。在rrnL的二級結構中有6個結構域，其中第3結構域缺失，rrnS的二級結構中有3個結構域。

2.4 系統發育分析

本研究以2條長蝽科的物種序列作為外群、蝽總科26個的物種作為內群，基于貝葉斯推斷BI（圖5、6）和最大似然法ML（圖7、8），進行蝽總科內的系統發育分析。選擇13個PCGs和13個PCGs的第1位和第2位密碼子（PCG12）串聯序列數據集構建系統發育樹。

基于13個PCGs在貝葉斯分析中，蝽總科內的關系是（異蝽科+（蝽科+（同蝽科+（（龜蝽科+盾蝽科）+（土蝽科+（兜蝽科+荔蝽科）））））），這與基于PCG12在貝葉斯分析得到的拓撲結構一致；在最大似然法分析中，基于13PCGs的蝽總科內的關系是（異蝽科+（同蝽科+（蝽科+（（龜蝽科+盾蝽科）+（土蝽科+（兜蝽科+荔蝽科））））））；基于PCG12的蝽總科內的關系卻是（蝽科+（（異蝽科+同蝽科）+（（龜蝽科+盾蝽科）+（土蝽科+（兜蝽科+荔蝽科））））））。

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

圖5 基于13個蛋白編碼基因利用貝葉斯分析構建的蝽總科系統發育樹Fig.5 Phylogenetic trees of Pentatomoidea constructed based on 13 protein-coding genes using Bayesian inference

圖7 基于13個蛋白編碼基因利用最大似然法構建的蝽總科系統發育樹Fig.7 Phylogenetic trees of Pentatomoidea constructed based on 13 protein-coding genes using RaxML inference

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

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

3 結論與討論

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

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

本研究中，BI和ML系統發育樹的結果均顯示，全蝽位于蝽科支系，在分子證據上說明全蝽隸屬于蝽科。同時，結果還很好地證明了蝽科的單系性，揭示了蝽科和（龜蝽科+盾蝽科）是姐妹群，這與以前的研究結果一致[18]。文中構建的4棵系統發育樹結果顯示，蝽科內部的短喙蝽亞科和益蝽亞科是單系性，但是舌盾蝽亞科和蝽亞科是并系，而且蝽亞科內部的物種關系更為混亂。1833年，YANG[25]基于盾片形狀將Graphosoma歸于蝽亞科Graphosomini，但是其他的研究中都是把它置于舌盾蝽亞科[26]。本研究中，也將Graphosoma歸于舌盾蝽亞科，但是它并不與舌盾蝽亞科的Scotinophara聚為一支。造成這樣的研究結果可能是因為蝽科物種的形態復雜，變異較大，特征不穩定，對其進行系統的劃分比較困難。學者劃分蝽科各亞科的標準也僅僅是根據此類昆蟲的外部形態的某一個或某幾個特征，這并不足以確定系統進化關系；也可能是因為本研究中選取的線粒體基因組基因含量不如核基因多，所以，不能充分表達各亞科之間的進化遺傳信息，因而劃分的結果會有分歧；當然本研究中也受到了樣本不足的限制。所以，要想解決蝽亞科和舌盾蝽亞科的分類關系存在的問題，還需要更多的形態和分子數據構建蝽科內各個亞科的系統發育位置。

本研究通過對全蝽線粒體的全基因組測序和生物信息學分析，發現全椿線粒體基因組全長15 479 bp，AT含量為76.18%。其中，13個編碼基因的起始密碼子多為ATN和TTG，終止密碼子為TAN或單個T堿基，atp8進化最快，cox1進 化 最慢；22個tRNA中trnA沒有TΨC環和TΨC臂，trnS1（GCU）沒有DHU臂，其他20個tRNA為典型的三葉草結構。系統發育樹結果顯示，全蝽位于蝽科支系且蝽科為單系。本研究結果豐富了蝽科線粒體基因組學和NCBI數據庫的基本數據，為進一步深化蝽科和蝽總科的系統發育關系討論提供了參考。