林興雨,宋 南

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

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

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

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

1 材料與方法

1.1 標本采集和DNA提取

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

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

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

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

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

1.4 多序列比對

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

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

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

2 結(jié)果與分析

2.1 線粒體基因組分析

茸毛小長蝽的線粒體基因組結(jié)構(gòu)圖利用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個蛋白質(zhì)編碼基因和14個轉(zhuǎn)運RNA基因(表1);L鏈(light strand)編碼4個蛋白質(zhì)編碼基因和8個轉(zhuǎn)運RNA基因以及2個核糖體RNA基因。此外,控制區(qū)位于rrnS和trnI之間,長度為2 011 bp,A、T、G、C含量分別為38.69%、36.20%、7.71%和17.40%。茸毛小長蝽的13個蛋白質(zhì)編碼基因、核糖體RNA和轉(zhuǎn)運RNA以及整個線粒體基因組的AT含量為73.1%～78.3%(表2),均具有較高的AT含量,其中,位于L鏈的8個轉(zhuǎn)運RNA基因的AT含量最低為73.1%,位于L鏈的4個蛋白質(zhì)編碼基因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 蛋白質(zhì)編碼基因分析

茸毛小長蝽線粒體基因組的13個蛋白質(zhì)編碼基因的序列在H鏈上的9個蛋白質(zhì)基因序列總長為6 762 bp,AT含量為75.3%,GC含量為24.7%。L鏈上的4個蛋白質(zhì)基因序列總長為4 218 bp,AT含量為78.3%,GC含量為21.7%。由上可知,位于H鏈和L鏈上的茸毛小長蝽的蛋白質(zhì)編碼基因均具有較高的AT含量。茸毛小長蝽的蛋白質(zhì)編碼基因的起始密碼子除nad2利用GTG和cox1利用TTG作為起始,其余11個蛋白質(zhì)編碼基因的起始密碼子都是以典型ATN密碼子作為起始,除cox1、cox2、nad3、cob和nad1以不完整的T作為結(jié)尾,其余8個蛋白質(zhì)編碼基因的終止密碼子都是以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個轉(zhuǎn)運RNA基因的序列長度在61～73 bp范圍內(nèi)。其中trnA基因最短為61 bp,而trnK基因最長為73 bp。位于H鏈有16個轉(zhuǎn)運RNA基因,序列全長為2 021 bp。AT含量為77.9%,GC含量為22.1%。位于L鏈6個轉(zhuǎn)運RNA基因,序列全長為523 bp,AT含量為73.1%,GC含量為26.9%。只有trnS1由于缺少二氫尿嘧啶臂(DHU臂)從而形成一個簡單的環(huán),形成不完整的三葉草結(jié)構(gòu)外,其余21個轉(zhuǎn)運RNA基因都形成了完整的三葉草結(jié)構(gòu)。此外,trnS1 的反密碼子不是常見的GCU而是UCG(圖2)。茸毛小長蝽的rrnL基因的序列全長為1 256 bp,AT含量為78.4%,GC含量為21.6%。rrnS基因的序列全長為765 bp,AT含量為77%,GC含量為23%。

圖2 茸毛小長蝽 22 個 tRNA基因的二級結(jié)構(gòu)Fig.2 Secondary structure of 22 tRNAs genes of N. graminicola

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

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

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

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

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

圖6 基于線粒體基因組中13個蛋白編碼基因氨基酸序列構(gòu)建的貝葉斯樹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.)之間的種間遺傳距離的結(jié)果顯示:茸毛小長蝽與N.fuscovittatus在Nysius屬的種間遺傳距離最近為0.143。其次是茸毛小長蝽與N.cymoides和N.sp.的遺傳距離為0.159和0.169,與N.plebeius的遺傳距離結(jié)果最遠為0.171(表4)。種間遺傳距離結(jié)果與系統(tǒng)發(fā)育樹分析的結(jié)果均支持茸毛小長蝽與N.fuscovittatus在Nysius屬中具有相對較近的親緣關(guān)系。

3 討論與結(jié)論

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

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

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

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

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