茶網(wǎng)蝽線粒體基因組全序列測定及系統(tǒng)發(fā)育分析

楊金宏，謝滿超，文欣茹，陳蕊茹，孔衛(wèi)青

1. 安康學(xué)院現(xiàn)代農(nóng)業(yè)與生物科技學(xué)院/陜西省茶葉省市共建重點實驗室，陜西 安康 725000；2. 安康學(xué)院陜西省蠶桑重點實驗室，陜西 安康 725000

茶網(wǎng)蝽（Stephanitis chinensis）為半翅目（Hemiptera）異翅亞目（Heteroptera）網(wǎng)蝽科（Tingidae）昆蟲，不完全變態(tài)，以我國西南地區(qū)的油茶和茶樹等經(jīng)濟(jì)作物為主要寄主。茶網(wǎng)蝽若蟲和成蟲以刺吸式口器從茶樹樹冠中、下層成熟葉片背面吸食汁液，其排泄物在葉片上形成黑色黏質(zhì)物[1]，受害茶樹樹勢衰弱，生長緩慢，發(fā)芽率低且芽包小，嚴(yán)重影響茶樹生長和茶葉產(chǎn)量[2]。

茶網(wǎng)蝽在陜南地區(qū)的茶葉產(chǎn)區(qū)一年發(fā)生2～3代，繁殖力強，群體數(shù)量龐大，種群擴(kuò)散能力強，新舊茶園均受危害。2010年陜西省漢中市首次報道茶園發(fā)生茶網(wǎng)蝽危害[3]，2012年擴(kuò)散至陜西省安康市，2017年蔓延至漢江流域和嘉陵江流域[4]。茶網(wǎng)蝽成蟲體長 3.1～4.2 mm，體小扁平，無趨性，初孵若蟲肉眼無法識別，亟待分子生物學(xué)鑒定方法[5]。

線粒體基因組（Mitochondrial genome，mitogenome）不易發(fā)生重組、倒位、易位等突變，結(jié)構(gòu)保守。線粒體基因進(jìn)化速率快，為核基因的 2～6倍，被廣泛應(yīng)用于種群多樣性鑒定、物種系統(tǒng)發(fā)育、物種地理學(xué)譜系和系統(tǒng)發(fā)生學(xué)等方面的研究[6-7]。全世界已知網(wǎng)蝽科昆蟲有260屬2 124種，但目前NCBI公布的線粒體基因組僅27個（截至2022年3月31號），網(wǎng)蝽科昆蟲的線粒體基因組大小差異較大，主要是由于控制區(qū)（Control region）數(shù)量、長度不等的重復(fù)序列造成[8]。本研究利用Illumina和Sanger測序，獲得了陜西安康茶區(qū)茶網(wǎng)蝽的線粒體基因組全序列及特征，并進(jìn)行了網(wǎng)蝽科的系統(tǒng)發(fā)育分析，旨在為進(jìn)一步對茶網(wǎng)蝽擴(kuò)散傳播途徑、防治及異翅亞目的分類研究提供依據(jù)。

1 材料與方法

1.1 試驗材料

茶網(wǎng)蝽成蟲于2021年5月17日采集自陜西省安康市紫陽縣和平茶廠茶園，采集的樣本存放在無水乙醇中，置于-20℃冰箱保存?zhèn)溆谩?/p>

1.2 基因組DNA提取與高通量測序

取乙醇浸泡的茶網(wǎng)蝽單個個體，ddH2O沖洗2遍，然后使用快速DNA提取檢測試劑盒[KG203，天根生化科技（北京）有限公司]按照說明書提取基因組 DNA，最后用100 μL ddH2O洗脫。使用 Multiskan FC微量酶標(biāo)儀（美國Thermo Fisher Scientific公司）對提取的基因組DNA進(jìn)行檢測。選取部分檢測合格的基因組DNA委托生工生物工程（上海）股份有限公司進(jìn)行建庫和Illumina Hiseq 2500平臺雙向高通量測序，測序長度為150 bp，目標(biāo)數(shù)據(jù)量4 Gb。

1.3 序列拼接組裝與控制區(qū)的PCR擴(kuò)增

對Illumina高通量測序數(shù)據(jù)進(jìn)行質(zhì)控，去除接頭和帶 N的序列，滑窗法去除得分低于20和長度小于 35 bp的序列。最后數(shù)據(jù)利用GetOrganelle[9]（K=31、55、85）和NOVOPlasty[10]（K=31）軟件進(jìn)行組裝，組裝結(jié)果使用Bandage[11]進(jìn)行可視化分析，手動去除可信度低和不確定的部分。在獲得的序列缺口處設(shè)計引物rrnS-F 和 trnM-R（表1），使用 Gloria Nova HS 2×高保真Taq DNA聚合酶（RK20715，武漢愛博泰克生物科技有限公司）進(jìn)行 PCR擴(kuò)增，PCR反應(yīng)程序：94℃預(yù)變性 2 min；94℃變性 30 s，55℃退火 30 s，72℃延伸 5 min，30個循環(huán)；72℃延伸10 min。使用1%瓊脂糖凝膠電泳檢測PCR產(chǎn)物，之后對產(chǎn)物進(jìn)行 Sanger雙向測序，并根據(jù)獲得序列設(shè)計測序引物（表1）進(jìn)行Gap補齊。

表1 引物序列Table 1 Primers sequence

1.4 基因組注釋和分析

使用在線軟件MITOS2[12]（http://mitos2.bioinf.uni-leipzig.de/index.py）對茶網(wǎng)蝽線粒體基因組進(jìn)行注釋，遺傳密碼選擇5（Invertebrate）。同時利用blastn和blastp比對網(wǎng)蝽科其他昆蟲的線粒體基因組數(shù)據(jù)確定蛋白質(zhì)編碼基因（Protein-coding genes，PCGs）、tRNA 和 rRNA的邊界位置。最后的注釋文件提交至NCBI數(shù)據(jù)庫，并利用在線 OGDRAW 軟件[13]（https://chlorobox.mpimp-golm.mpg.de/OGDr aw.html）繪制結(jié)構(gòu)圖。tRNA的二級結(jié)構(gòu)預(yù)測使用在線軟件tRNAscan-SE[14]（http://lowelab.ucsc.edu/tRNAscan-SE），未成功預(yù)測的二級結(jié)構(gòu)進(jìn)行手工繪制。

使用在線工具 Tandem Repeats Finder[15]（https://tandem.bu.edu/trf/trf.html）查找控制區(qū)中的串聯(lián)重復(fù)序列，使用MISA[16]檢測其中的簡單重復(fù)序列（Simple sequence repeats，SSRs），參數(shù)設(shè)置為1～6堿基單元的最少重復(fù)數(shù)分別為10、5、4、3、3、3個。莖環(huán)結(jié)構(gòu)預(yù)測使用在線工具M(jìn)Fold[17]（http://www.unafold.org/mfold/applications/rna-folding-form.php）。

使用 Mega 11軟件[18]統(tǒng)計基因組及各區(qū)域 4種堿基的組成和蛋白質(zhì)編碼基因的密碼子情況，并通過公式AT skew=(A-T)/(A＋T)和 GC skew=(G-C)/(G＋C)[19]計算核苷酸組成的偏向性。

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

從 NCBI數(shù)據(jù)庫下載異翅亞目昆蟲線粒體基因組，分別提取13個PCGs序列，然后利用 macse_v2.05軟件[20]基于編碼框架模式對序列進(jìn)行比對，并把對齊的DNA合并成一個序列。利用 raxmlHPC[21]構(gòu)建極大似然（Maximum likelihood，ML）發(fā)育樹，GTRGAMMA模式，抽樣 1 000次計算bootstrap值，展示高于50%的位點。

2 結(jié)果與分析

2.1 茶網(wǎng)蝽線粒體基因組的測序與組裝

Illumina測序數(shù)據(jù)經(jīng)質(zhì)控，獲得2 025 131條序列，GetOrganelle和NOVOPlasty組裝軟件獲得的線粒體基因組序列完全一致的部分長度為14 561 bp，平均深度為415.04倍?？刂茀^(qū)序列的組裝得到多種結(jié)果，可能存在重復(fù)結(jié)構(gòu)。進(jìn)一步利用引物rrnS-F和trnM-R進(jìn)行PCR擴(kuò)增，獲得1條4 000 bp左右的條帶（圖1）。使用設(shè)計的多條測序引物對控制區(qū)進(jìn)行Sanger測序，拼接測序結(jié)果，獲得序列長度為4 040 bp，與Illumina測序組裝數(shù)據(jù)拼接成環(huán)，獲得茶網(wǎng)蝽線粒體基因組全長為 18 085 bp（GenBank登錄號：OM397399）。

圖1 rrnS-F和trnM-R引物的PCR擴(kuò)增Fig. 1 The PCR amplification of rrnS-F and trnM-R

2.2 茶網(wǎng)蝽線粒體基因組注釋與結(jié)構(gòu)

分析發(fā)現(xiàn)，茶網(wǎng)蝽線粒體基因組為典型的雙鏈閉合環(huán)狀結(jié)構(gòu)，包含37個基因，分別為13個PCGs、22個tRNA基因和2個rRNA基因（圖2，表2）?；蚪M控制區(qū)為3 678 bp，其后為trnI(GAU)、trnQ(UUG)、trnM(CAU)和nad2基因，與昆蟲線粒體祖先基因順序（Ancestral gene order）相同[22]。

圖2 茶網(wǎng)蝽線粒體基因組圖譜Fig. 2 Map of S. chinensis mitogenome

表2 茶網(wǎng)蝽線粒體基因組結(jié)構(gòu)Table 2 The mitogenome organization of S. chinensis

茶網(wǎng)蝽的線粒體基因組的37個基因共存在18處基因重疊，最長的兩處位于nad1和trnL1(UAG)、trnF(GAA)和nad5之間，重疊長度分別為-36 bp和-32 bp，其他位置的重疊長度均在-8 bp以內(nèi)?；蜷g隔現(xiàn)象少于基因重疊，37個基因僅存在8處間隔，其中最長兩處位于trnV(UAC)和rrnS、trnL1(UAG)和rrnL之間，長度分別為47 bp和25 bp，其他間隔最長8 bp，另有無間隔和重疊的區(qū)域10處（表2）。

2.3 茶網(wǎng)蝽線粒體基因組核苷酸組成

茶網(wǎng)蝽線粒體基因組的 AT含量為78.09%，具有明顯的AT偏好性（表3）。rRNA基因的AT含量最高，為80.84%，AT富集區(qū)為79.99%。PCGs的AT含量最低，為76.76%，其中atp8基因的AT含量最高，達(dá)84.28%，其次是nad2和nad4l，分別為 82.21%和80.70%；AT含量最低的是cox1和cob，分別為69.08%和72.21%。核苷酸組成在線粒體基因組不同鏈間也不對稱，多數(shù)編碼基因所在的J鏈（Majority strand）的PCGs多為AT偏移和CG偏移，而N鏈（Minority strand）的PCGs和rRNA基因，為TA偏移和GC偏移（表2）。

表3 茶網(wǎng)蝽線粒體基因組核苷酸組成Table 3 Nucleotide composition of the mitogenome of S. chinensis

2.4 蛋白質(zhì)編碼基因

茶網(wǎng)蝽線粒體基因組的13個PCGs的長度在 159～1 536 bp，其中nad1、nad4、nad4l和nad5位于基因組N鏈，其他9個均位于J鏈（表2）。13個基因中，cox1、cox3、atp6、nad5、cob和nad1以ATG起始，atp8、nad3、nad4l和nad2以 ATT起始，cox2、nad4和nad6以ATA起始；以TAA為終止密碼子的基因最多，共有9個，cox2、atp6和cox3以不完全的T終止，cob以TAG終止。同時，PCGs之間存在相互重疊的現(xiàn)象，atp8和atp6基因間重疊了7 bp，nad4和nad4l基因、nad6和cob基因分別重疊了1 bp和2 bp。

進(jìn)一步分析發(fā)現(xiàn)，茶網(wǎng)蝽線粒體PCGs使用最多的 5個密碼子分別是 UUA、AUU、UUU、AUA和 AAU，均由 A和 T組成（表4）。編碼蛋白的氨基酸組成及含量也有較大差異，使用頻率最高的為亮氨酸（Leu），占比約14.00%，其次為異亮氨酸（Ile）、苯丙氨酸（Phe）、絲氨酸（Ser）、甲硫氨酸（Met）、天冬酰胺（Asn）。其他氨基酸占比均低于5%，其中精氨酸（Arg）、半胱氨酸（Cys）、谷氨酰胺（Gln）、組氨酸（His）、天冬氨酸（Asp）的占比低于2%。

表4 茶網(wǎng)蝽線粒體基因組蛋白編碼基因相對同義密碼子使用度Table 4 Relative synonymous codon usage (RSCU) in the mitogenome of S. chinensis

2.5 tRNA和rRNA基因

茶網(wǎng)蝽線粒體基因組包含22個tRNA基因（圖3），長度在63～70 bp。8個tRNA基因（trnC，trnF，trnH，trnL1，trnP，trnQ，trnV，trnY）位于基因組N鏈，其余14個位于J鏈（表2）。預(yù)測 tRNA基因的二級結(jié)構(gòu)發(fā)現(xiàn)，trnS1(GCU)缺少DHU臂，其他tRNA均能形成典型的三葉草結(jié)構(gòu)（圖3），該現(xiàn)象存在于多數(shù)后生動物中[23-24]。同時，部分tRNA的二級結(jié)構(gòu)還存在G-U和U-U等非經(jīng)典配對，該現(xiàn)象在其他昆蟲類群 tRNA二級結(jié)構(gòu)中也較為常見[25-26]。

圖3 茶網(wǎng)蝽線粒體基因組tRNA預(yù)測二級結(jié)構(gòu)Fig. 3 Predicted secondary structures of tRNA genes in mitogenome of S. chinensis

茶網(wǎng)蝽線粒體基因組含有 2個 rRNA基因，rrnL位于trnL1(UAG)和trnV(UAC)之間，rrnS位于trnV(UAC)和控制區(qū)之間，長度分別為1 237 bp和767 bp，AT含量分別為80.22%和81.88%，具有明顯的AT偏向性。

2.6 控制區(qū)分析

茶網(wǎng)蝽線粒體基因組的控制區(qū)位于rrnS和trnI(GAU)基因之間，長度為 3 678 bp，根據(jù)結(jié)構(gòu)將其分為兩部分（圖4）。第一部分位于序列的前端（1～2 447 bp，front-end region，F(xiàn)ER），存在以下結(jié)構(gòu)：（1）3組非串聯(lián)重復(fù)序列。第一組長度為1 540 bp（R1），包含2個498 bp的完整重復(fù)單元和 2個該單元的部分序列（477 bp和 63 bp）；第二組長度為 154 bp（R2），包含2個77 bp的完整重復(fù)單元；第三組長度為582 bp（R3），由3個194 bp的完整重復(fù)單元組成。（2）4處(TTAG)n，其中(TTAG)7和(TTAG)9分別連接非串聯(lián)重復(fù) R2的 2個單元，另 2處均為(TTAG)5，位于 R2重復(fù)的2個單元中。第二部分位于序列的后端（2 448～3 678 bp，back-end region，BER），為串聯(lián)重復(fù)區(qū)（R4），包含5個長度為237 bp的完整重復(fù)單元和1個該單元的部分序列。同時序列存在5種莖環(huán)結(jié)構(gòu)，其中2種5處為單一環(huán)結(jié)構(gòu)，另3種11處為多環(huán)結(jié)構(gòu)。

圖4 茶網(wǎng)蝽線粒體基因組控制區(qū)結(jié)構(gòu)Fig. 4 Structure of mitogenome control region of S. chinensis

進(jìn)一步分析該區(qū)域的AT含量為79.99%，高于Pseudacysta perseae、Agramma hupehanum和楊柳網(wǎng)蝽（Metasalis populi）該區(qū)的AT含量，低于其他網(wǎng)蝽科昆蟲。FER和BER的AT含量分別為79.85%和80.26%（表3），兩段區(qū)域?qū)Σ煌塑账岬氖褂闷靡灿兴煌現(xiàn)ER的 AT偏度為-0.12，發(fā)生 TA偏移，BER的AT使用量相當(dāng)，有GC偏移傾向。

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

以雙線銅翅沫蟬（Aeneolamia contigua）為外群，構(gòu)建了基于異翅亞目41種昆蟲的線粒體PCGs的系統(tǒng)發(fā)育樹（圖5）。結(jié)果表明，茶網(wǎng)蝽與直脊冠網(wǎng)蝽（Stephanitis mendica）的親緣關(guān)系最近，其次為楊柳網(wǎng)蝽和Pseudacysta perseae。網(wǎng)蝽科聚為一簇，位于發(fā)育樹的根部，是臭蟲次目（Cimicomorpha）盲蝽科（Miridae）、花蝽科（Anthocoridae）、姬蝽科（Nabidae），以及蝽次目（Pentatomomorpha）、黽蝽次目（Gerromorpha）、歇蝽次目（Nepomorpha）、奇蝽次目（Enicocephalomorpha）的姐妹群。奇蝽次目和黽蝽次目距離最近，與歇蝽次目共同聚合在一簇，蝽次目單獨聚合為一個分支。系統(tǒng)發(fā)育樹支持了歇蝽次目、奇蝽次目、黽蝽次目和蝽次目的單系性，而臭蟲次目為并系群。

圖5 基于41種昆蟲線粒體蛋白質(zhì)編碼基因序列構(gòu)建的ML系統(tǒng)發(fā)育樹Fig. 5 Maximum likelihood phylogenetic tree of 41 insects based on the PCG sequences of mitogenome

3 討論

本研究結(jié)合Illumina測序和Sanger測序，獲得安康市茶網(wǎng)蝽線粒體基因組全長18 085 bp，是目前 NCBI登錄的網(wǎng)蝽科昆蟲線粒體全序列（15 208～16 667 bp）中長度最長的?？刂茀^(qū)是昆蟲線粒體基因組大小變異的主要來源，一般與線粒體基因組的大小呈正相關(guān)，控制區(qū)的長度差異可以發(fā)生在不同物種甚至不同個體之間[27]。目前昆蟲中已報道的最短的控制區(qū)僅為70 bp，而象鼻蟲的控制區(qū)長度達(dá) 9～13 kb[28]。已報道的網(wǎng)蝽科昆蟲線粒體基因組的控制區(qū)為733～2 215 bp，而本研究中安康市茶網(wǎng)蝽的控制區(qū)為3 678 bp。

控制區(qū)是線粒體基因組復(fù)制和轉(zhuǎn)錄的主要調(diào)控區(qū)，包含線粒體基因組復(fù)制相關(guān)信息[29]。黽蝽科圓臀大黽蝽（Aquarius paludum）線粒體基因組控制區(qū)存在 TATA motif、poly T stretch、莖環(huán)結(jié)構(gòu)等保守結(jié)構(gòu)元件[30]，地長蝽科大黑毛肩長蝽（Neolethaeus assamensis）的線粒體基因組控制區(qū)存在poly A結(jié)構(gòu)、莖環(huán)結(jié)構(gòu)、串聯(lián)重復(fù)區(qū)等[31]。本研究發(fā)現(xiàn)，茶網(wǎng)蝽基因組中有多種重復(fù)序列和莖環(huán)結(jié)構(gòu)，而重慶市茶網(wǎng)蝽[32]僅有 R4重復(fù)序列，沒有 FER區(qū)，對于不同地域茶網(wǎng)蝽基因組之間的差異，有待進(jìn)一步研究。

本研究茶網(wǎng)蝽的線粒體基因組所有區(qū)域均有明顯的 AT偏好性，AT含量最高的是rRNA，這與網(wǎng)蝽科昆蟲楊柳網(wǎng)蝽、菊方翅網(wǎng)蝽（Corythucha marmorata）、古無孔網(wǎng)蝽（Dictyla platyoma）中的情況相同。蛋白質(zhì)編碼基因的AT含量最低，這有利于物種編碼其所需蛋白質(zhì)的功能[33]。cox1基因是常用的DNA條形碼基因，應(yīng)用于昆蟲的識別鑒定[34]。相對較高的 GC含量有利于降低基因的突變率，保持基因穩(wěn)定性。Yang等[8]研究發(fā)現(xiàn)，網(wǎng)蝽科昆蟲線粒體cox1基因的進(jìn)化速率較低，atp8基因最高，本研究中安康市茶網(wǎng)蝽cox1基因的GC含量最高，atp8基因的GC含量最低，與其結(jié)論一致。與多數(shù)動物一樣，茶網(wǎng)蝽線粒體基因組不同鏈的核苷酸組成也不對稱[19]。

網(wǎng)蝽科昆蟲線粒體基因組中絕大部分蛋白基因以ATN為起始密碼子，使用最多的是ATG，其中cox1和atp6基因的起始密碼子均為 ATG，ATA的使用頻次略低于 ATT。終止密碼子以TAA為主，cox1、atp8、nad6和nad1基因的終止密碼子都是TAA；其次為T，cox2基因均以單個T作為終止密碼子，cox3和nad5基因使用 T作為終止密碼子的現(xiàn)象較多，較少使用 TA或 TAG。昆蟲線粒體基因組以 T或TA作為終止密碼子現(xiàn)象比較常見，推測其通過轉(zhuǎn)錄后 mRNA多腺苷酸化形成完全終止密碼子TAA[35]。

臭蟲次目是異翅亞目最大的次目，有 17個科超過20 000種昆蟲，其系統(tǒng)發(fā)育研究目前沒有一致的結(jié)論。利用形態(tài)性狀、轉(zhuǎn)錄組和核基因等研究發(fā)現(xiàn)，臭蟲次目昆蟲為單系[36-38]，而使用臭蟲次目昆蟲線粒體基因組數(shù)據(jù)進(jìn)行系統(tǒng)發(fā)育研究，鑒定其為并系群，尤其是網(wǎng)蝽科昆蟲的數(shù)據(jù)，該現(xiàn)象更顯著[8,39-40]。Liu等[39]基于線粒體基因組數(shù)據(jù)分析異翅亞目的系統(tǒng)發(fā)育發(fā)現(xiàn)，網(wǎng)蝽科和奇蝽次目、鞭蝽次目、黽蝽次目聚為一簇且接近發(fā)育樹的根部，是其他異翅亞目昆蟲的姐妹群。本研究利用線粒體PCGs進(jìn)行的系統(tǒng)發(fā)育重建結(jié)果，能很好地支持奇蝽次目、黽蝽次目、歇蝽次目和蝽次目的單系性。但臭蟲次目的情況比較復(fù)雜，其姬蝽科和花蝽科聚合在一起，形成與奇蝽次目、黽蝽次目和歇蝽次目并列的分支，而其盲蝽科和網(wǎng)蝽科的位置接近發(fā)育樹的根部，這與Liu等[39]的研究結(jié)果相似。但本研究接近發(fā)育樹根部的幾個分支的支持率均低于50%，這可能與網(wǎng)蝽科昆蟲線粒體基因組序列的變異率較大有關(guān)[8,39]，較大堿基差異不容易區(qū)分系統(tǒng)發(fā)育關(guān)系，更適合于親緣關(guān)系近的物種的分化檢測。