茸毛小長蝽的線粒體基因組測序和分析

林興雨,宋 南

(河南農業(yè)大學 植物保護學院,河南 鄭州 450046)

長蝽總科(Lygaeoidea)是半翅目(Hemiptera)蝽次目(Pentatomomorpha)中物種數(shù)量較為豐富的類群之一。目前,全世界有17個科超過4 290個物種已經(jīng)被描述[1]。長蝽總科的大多數(shù)物種屬于植食性昆蟲,也有少部分是天敵昆蟲。茸毛小長蝽(Nysiusgraminicola)隸屬于長蝽總科長蝽科(Lygaeidae),是一些蔬菜和水果等作物上的重要害蟲,尤其是在高粱、葡萄、番茄和桃等作物上危害較為嚴重[2]。

隨著廉價且快速的二代測序技術的飛速發(fā)展,越來越多的昆蟲高質量的線粒體基因組通過測序得到。昆蟲線粒體基因組具有母系遺傳以及進化速率較快等獨特的優(yōu)勢,作為一種分子遺傳標記被廣泛應用于昆蟲的系統(tǒng)發(fā)育研究中[3-5]。截至目前,NCBI數(shù)據(jù)庫中已經(jīng)有超過560條蝽次目的部分或完整的線粒體基因組序列被釋放。近年來,利用分子學手段對蝽次目的高階元的系統(tǒng)發(fā)育研究較多[6-7],而針對長蝽總科的內部的系統(tǒng)發(fā)育關系卻研究較少。對長蝽總科系統(tǒng)發(fā)育研究的匱乏,不利于人們從分子學角度深入理解長蝽總科的系統(tǒng)發(fā)育關系。為了使用線粒體基因組數(shù)據(jù)更好地去了解長蝽總科的系統(tǒng)發(fā)育關系,本研究利用二代測序技術首次獲得了茸毛小長蝽的線粒體全基因組序列,對其線粒體基因組序列進行了詳細注釋,分別分析了它的蛋白質編碼基因、轉運RNA基因、RNA基因,以及整個線粒體基因組的核苷酸組成、AT偏斜和GC偏斜、密碼子使用頻率情況、種間遺傳距離,預測了茸毛小長蝽的轉運RNA的二級結構。

此外,以GenBank數(shù)據(jù)庫中已經(jīng)釋放的長蝽總科的17個科(象蝽科Idiostolidae、梭長蝽科Pachygronthidae、地長蝽科Rhyparochromidae、室翅長蝽科Heterogastridae、保長蝽科Artheneidae、尖長蝽科Oxycarenidae、皮蝽科Piesmatidae、Cryptorhamphidae、束蝽科Colobathristidae、Meschiidae、蹺蝽科Berytidae、莎長蝽科Cymidae、Ninidae、桿長蝽科Blissidae、束長蝽科Malcidae、大眼長蝽科Geocoridae和長蝽科)52個物種的線粒體基因組序列作為內群,緣蝽總科Coreoidea的2個物種(Melanacanthusmarginatus和Notobitusmontanus)線粒體基因組序列作為外群,本研究基于線粒體基因組中的13個蛋白質編碼基因的核苷酸序列矩陣和13個蛋白質編碼基因的氨基酸序列矩陣,分別使用最大似然法(maximum likelihood,ML)和貝葉斯法(Bayesian inference,BI)重建長蝽總科內部的系統(tǒng)發(fā)育關系,明確茸毛小長蝽的系統(tǒng)發(fā)育位置,以期為更好地理解長蝽總科和茸毛小長蝽的系統(tǒng)發(fā)育關系提供線粒體基因組數(shù)據(jù)。

1 材料與方法

1.1 標本采集和DNA提取

本研究所使用的茸毛小長蝽的標本于2018年8月采集于河南省鄭州市。將采集到的茸毛小長蝽標本浸泡在95%乙醇中,存放于河南農業(yè)大學植物保護學院標本館并置于-20 ℃低溫冰箱保存直到用于DNA提取。在實驗室內使用天根公司的動物基因組DNA提取試劑盒(TIANamp Genomic DNA Kit,中國)按照說明書步驟提取茸毛小長蝽胸部肌肉的總DNA。使用NanoDrop 2000分光光度計對總DNA濃度檢測以及配置1.5%的瓊脂糖凝膠電泳進行總DNA的質量檢測。

1.2 二代測序及線粒體基因組的組裝

將檢測質量合格且濃度大于50 ng/μL的高質量的DNA送往上海派森諾生物科技有限公司進行350 bp的小片段文庫構建并進行二代測序,基于Illumina HiSeq 2500測序平臺進行雙末端(Paired-end,PE)測序。最后,將下機的原始數(shù)據(jù)去除接頭,刪除低質量序列,得到所需序列。利用GetOrganelle v1.7.5.2軟件[8]對測序獲得的clean reads進行線粒體基因組的組裝和拼接,從而得到高質量的線粒體全基因組數(shù)據(jù)。

1.3 線粒體基因組、注釋和分析

將組裝得到的茸毛小長蝽的線粒體基因組數(shù)據(jù)在MITOS[9]網(wǎng)站上進行初步的基因注釋,具體參數(shù)設置為:Genetic Code:05-inverterbrate;Reference:RefSeq 63 Metazoa,其余均按照MITOS的默認參數(shù)進行設置。通過MITOS預測獲得茸毛小長蝽的22個轉運RNA的二級結構。然后,通過近緣種線粒體基因組的序列比對和驗證,手動調整蛋白質編碼基因和RNA編碼基因的基因邊界。將注釋完整的線粒體基因組序列上傳至GenBank,得到茸毛小長蝽的GenBank登錄號:OQ553818。利用系統(tǒng)發(fā)育分析學軟件MEGA 7.0[10]分析茸毛小長蝽線粒體基因組的蛋白質編碼基因、轉運RNA基因和核糖體RNA基因的核苷酸組成以及蛋白質編碼基因的密碼子的使用頻率情況。基于A、T含量之差與A、T含量之和的比計算AT偏斜,基于G、C含量之差與G、C含量之和的比計算GC偏斜。

1.4 多序列比對

采用MAFFT軟件中的E-INS-I策略進行蛋白質編碼基因的多序列比對[11]。使用TimAl v1.4程序對多序列比對中質量較差的序列進行剔除[12]。最后,運用FASconCAT-G_v1.04程序進行串聯(lián)獲得系統(tǒng)發(fā)育分析所用的數(shù)據(jù)矩陣[13]。用于系統(tǒng)發(fā)育分析的數(shù)據(jù)矩陣包括:(1) 13 個蛋白質編碼基因的核苷酸序列構建的數(shù)據(jù)矩陣分別用于最大似然法和貝葉斯法系統(tǒng)發(fā)育分析; (2) 13 個蛋白質編碼基因的氨基酸序列構建的數(shù)據(jù)矩陣分別用于最大似然法和貝葉斯法系統(tǒng)發(fā)育分析;

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

本研究使用二代測序得到的茸毛小長蝽的線粒體基因組數(shù)據(jù),并結合GenBank中獲得的54種昆蟲(52種長蝽總科昆蟲作為內群,以及緣蝽總科2個物種作為外群)的線粒體基因組數(shù)據(jù)的蛋白質編碼基因構建的2個數(shù)據(jù)矩陣。通過最大似然法(Maximum Likelihood, ML)和貝葉斯法(Mr Bayesian, BI)重建長蝽總科的系統(tǒng)發(fā)育關系。首先,使用 IQ-TREE 2.0.6[14]進行最大似然法分析:由軟件自動篩選蛋白質編碼基因的核苷酸序列矩陣最適模型為GTR+F+I+G4模型,蛋白質編碼基因的氨基酸序列矩陣最適模型為mtZOA+F+I+G4模型,系統(tǒng)發(fā)育樹的節(jié)點支持值通過自舉檢驗置信度進行評估(運行次數(shù)設置10 000)。使用PhyloBayes v3.3 進行貝葉斯法分析[15]。首先,執(zhí)行運算兩次,每次運算包含四個鏈。核苷酸數(shù)據(jù)矩陣使用CAT-GTR模型以及氨基酸數(shù)據(jù)矩陣使用CAT-mtART位點異質模型來進行系統(tǒng)發(fā)育分析,利用PhyloBayes軟件包中的bpcomp程序評估鏈間的趨同一致性。以最大差異值小于0.1作為分析達到有效的判定標準。最后,使用bpcomp程序構建50%合意樹。

2 結果與分析

2.1 線粒體基因組分析

茸毛小長蝽的線粒體基因組結構圖利用OGDRAW Version 1.1獲得[16](圖1)。茸毛小長蝽的線粒體全基因組長度為16 760 bp,其中A、T、G、C含量分別為43.10%、33.20%、9.50%和14.20%,其線粒體基因組包含37個基因和1個非編碼控制區(qū),H鏈(heavy strand)編碼9個蛋白質編碼基因和14個轉運RNA基因(表1);L鏈(light strand)編碼4個蛋白質編碼基因和8個轉運RNA基因以及2個核糖體RNA基因。此外,控制區(qū)位于rrnS和trnI之間,長度為2 011 bp,A、T、G、C含量分別為38.69%、36.20%、7.71%和17.40%。茸毛小長蝽的13個蛋白質編碼基因、核糖體RNA和轉運RNA以及整個線粒體基因組的AT含量為73.1%～78.3%(表2),均具有較高的AT含量,其中,位于L鏈的8個轉運RNA基因的AT含量最低為73.1%,位于L鏈的4個蛋白質編碼基因AT含量最高為78.3%。此外,AT偏斜和GC偏斜范圍在-0.326～0.130和-0.198～0.291。

表1 茸毛小長蝽線粒體基因組注釋Table 1 Annotation of the mitochondrial genome of N. graminicola

表2 茸毛小長蝽線粒體基因組的核苷酸組成和偏斜Table 2 Nucleotide composition and skewness of the N. graminicola mitochondrial genome

2.2 蛋白質編碼基因分析

茸毛小長蝽線粒體基因組的13個蛋白質編碼基因的序列在H鏈上的9個蛋白質基因序列總長為6 762 bp,AT含量為75.3%,GC含量為24.7%。L鏈上的4個蛋白質基因序列總長為4 218 bp,AT含量為78.3%,GC含量為21.7%。由上可知,位于H鏈和L鏈上的茸毛小長蝽的蛋白質編碼基因均具有較高的AT含量。茸毛小長蝽的蛋白質編碼基因的起始密碼子除nad2利用GTG和cox1利用TTG作為起始,其余11個蛋白質編碼基因的起始密碼子都是以典型ATN密碼子作為起始,除cox1、cox2、nad3、cob和nad1以不完整的T作為結尾,其余8個蛋白質編碼基因的終止密碼子都是以TAA作為終止密碼子。在茸毛小長蝽線粒體基因組的密碼子使用頻率中,密碼子AAA使用次數(shù)最多為427次,密碼子GCG使用次數(shù)最少,僅使用2次(表3)。

表3 茸毛小長蝽線粒體基因組的相對密碼子使用頻率Table 3 Relative synonymous codon usage of the mitochondrial genome of N. graminicola

2.3 tRNA和rRNA基因分析

茸毛小長蝽線粒體基因組的22個轉運RNA基因的序列長度在61～73 bp范圍內。其中trnA基因最短為61 bp,而trnK基因最長為73 bp。位于H鏈有16個轉運RNA基因,序列全長為2 021 bp。AT含量為77.9%,GC含量為22.1%。位于L鏈6個轉運RNA基因,序列全長為523 bp,AT含量為73.1%,GC含量為26.9%。只有trnS1由于缺少二氫尿嘧啶臂(DHU臂)從而形成一個簡單的環(huán),形成不完整的三葉草結構外,其余21個轉運RNA基因都形成了完整的三葉草結構。此外,trnS1 的反密碼子不是常見的GCU而是UCG(圖2)。茸毛小長蝽的rrnL基因的序列全長為1 256 bp,AT含量為78.4%,GC含量為21.6%。rrnS基因的序列全長為765 bp,AT含量為77%,GC含量為23%。

圖2 茸毛小長蝽 22 個 tRNA基因的二級結構Fig.2 Secondary structure of 22 tRNAs genes of N. graminicola

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

本研究結合長蝽總科的17個科的52個物種作為內群,緣蝽總科Coreoidea的2個物種作為外群,基于它們的蛋白質編碼基因構建了核苷酸數(shù)據(jù)矩陣和氨基酸數(shù)據(jù)矩陣,分別使用最大似然法(ML)和貝葉斯法(BI)重建長蝽總科的系統(tǒng)發(fā)育關系(圖3～圖6)。系統(tǒng)發(fā)育分析結果一致顯示:梭長蝽科、室翅長蝽科、蹺蝽科、桿長蝽科和束長蝽科為單系群。地長蝽科和長蝽科為非單系群。然而,由于象蝽科、尖長蝽科、保長蝽科、皮蝽科、Cryptorhamphidae、Meschiidae、莎長蝽科、Ninidae、束蝽科和大眼長蝽科各僅有一個物種的線粒體基因組數(shù)據(jù)被釋放,系統(tǒng)發(fā)育關系仍是不確定的。

圖3 基于線粒體基因組中13個蛋白編碼基因核苷酸序列構建的最大似然樹Fig.3 Maximum likelihood tree inferred from the nucleotide sequences of 13 protein-coding genes in the mitochondrial genome

圖4 基于線粒體基因組中13個蛋白編碼基因氨基酸序列構建的最大似然樹Fig.4 Maximum likelihood tree inferred from the amino acid sequences of 13 protein-coding genes in the mitochondrial genome

圖5 基于線粒體基因組中13個蛋白編碼基因核苷酸序列構建的貝葉斯樹Fig.5 Bayesian tree inferred from the nucleotide sequences of 13 protein-coding genes the mitochondrial genome

圖6 基于線粒體基因組中13個蛋白編碼基因氨基酸序列構建的貝葉斯樹Fig.6 Bayesian tree inferred from the amino acid sequences of 13 protein-coding genes the mitochondrial genome

利用 Kimura-2-Parameter 模型計算了Nysius屬5個物種(N.cymoides、N.fuscovittatus、茸毛小長蝽N.graminicola、N.plebeius和N.sp.)之間的種間遺傳距離的結果顯示:茸毛小長蝽與N.fuscovittatus在Nysius屬的種間遺傳距離最近為0.143。其次是茸毛小長蝽與N.cymoides和N.sp.的遺傳距離為0.159和0.169,與N.plebeius的遺傳距離結果最遠為0.171(表4)。種間遺傳距離結果與系統(tǒng)發(fā)育樹分析的結果均支持茸毛小長蝽與N.fuscovittatus在Nysius屬中具有相對較近的親緣關系。

3 討論與結論

本研究利用二代測序技術獲得了茸毛小長蝽的線粒體基因組序列,分析了它的核苷酸組成等線粒體基因組結構,詳細對茸毛小長蝽的線粒體基因組進行注釋。茸毛小長蝽的線粒體基因組包含37個基因和一段非編碼控制區(qū)。線粒體基因排列順序與祖先昆蟲的排列順序一致,沒有出現(xiàn)基因重排等特殊現(xiàn)象[3]。茸毛小長蝽的線粒體基因組的核苷酸組成與其它半翅目昆蟲的核苷酸組成相同,都富含AT堿基[17]。此外,預測并繪制了茸毛小長蝽的22個轉運RNA基因的二級結構。在所有轉運RNA基因的二級結構中,只有trnS1基因缺少二氫尿嘧啶臂(DHU臂),其余轉運RNA基因都具有完整的二級結構。這種情況在其它昆蟲的轉運RNA基因也出現(xiàn)過[18-24]。

在蝽次目的系統(tǒng)發(fā)育研究中,Li等[20]使用18S和cox1基因對蝽次目的系統(tǒng)發(fā)育關系進行初步研究的結果對長蝽總科的單系性是模糊的[7]。Yuan等[25]結合蝽次目的5個總科26個物種(其中涉及長蝽總科的5個科的6個物種)以及2個外群的線粒體基因組的蛋白質編碼基因的核苷酸序列的系統(tǒng)發(fā)育結果顯示長蝽總科的5個科(長蝽科、束長蝽科、蹺蝽科、大眼長蝽科和束蝽科)為單系群[25]。此外,Huang等[26]使用長蝽總科的6個科(長蝽科、大眼長蝽科、蹺蝽科、束長蝽科、地長蝽科和束蝽科)12個物種作內群以及2個物種做外群的蛋白質編碼基因與核糖體基因利用最大似然法和Carapelli等[27]基于長蝽總科的6個科(長蝽科、束長蝽科、大眼長蝽科、地長蝽科、蹺蝽科和束蝽科)13個物種作內群和2個物種作為外群的蛋白質編碼基因利用貝葉斯法對長蝽總科內部的系統(tǒng)發(fā)育關系的研究結果一致顯示地長蝽科為非單系群,其余5個科均為單系群。

本研究利用二代測序技術獲得的茸毛小長蝽的線粒體基因組序列,并結合長蝽總科目前已經(jīng)公布的所有線粒體基因組序列(包括17個科的52個物種)作為內群,2個緣蝽總科的物種作為群的線粒體基因組的蛋白質編碼基因構建了蛋白質編碼基因核苷酸序列矩陣和氨基酸數(shù)據(jù)矩陣,分別使用最大似然法和貝葉斯法進行系統(tǒng)發(fā)育分析,結果顯示梭長蝽科、室翅長蝽科、蹺蝽科、桿長蝽科和束長蝽科為單系群。地長蝽科為非單系群,這一結果與Huang等[26]和Carapelli等[27]的系統(tǒng)發(fā)育結果一致。然而,核苷酸數(shù)據(jù)矩陣利用最大似然法和貝葉斯法以及氨基酸數(shù)據(jù)矩陣利用貝葉斯法系統(tǒng)發(fā)育分析中由于大眼長蝽科的嵌入長蝽科而導致長蝽科為非單系群。此外,氨基酸數(shù)據(jù)矩陣利用最大似然法系統(tǒng)發(fā)育結果顯示由于束蝽科的嵌入也導致了長蝽科為非單系群。

Nysius屬5個物種的種間遺傳距離與系統(tǒng)發(fā)育分析結果都顯示茸毛小長蝽與N.fuscovittatus在Nysius屬中具有較近的親緣關系。這可能是由于它們在形態(tài)特征上以及自然界中的生物學習性上擁有較多的共同特點而形成的。

本研究利用二代測序技術首次獲得了茸毛小長蝽的線粒體全基因組,分析了它的線粒體基因組結構特點,并基于數(shù)據(jù)庫已經(jīng)公布的所有的長蝽總科的線粒體基因組數(shù)據(jù)重建了目前較為全面的長蝽總科的系統(tǒng)發(fā)育關系。使用Kimura-2-Parameter 模型和兩種系統(tǒng)發(fā)育分析方法明確了茸毛小長蝽的系統(tǒng)發(fā)育位置,為后續(xù)長蝽總科的系統(tǒng)發(fā)育研究以及茸毛小長蝽的種群遺傳學等研究提供線粒體基因組數(shù)據(jù)。