茶网蝽线粒体基因组全序列测定及系统发育分析

杨金宏，谢满超，文欣茹，陈蕊茹，孔卫青

1. 安康学院现代农业与生物科技学院/陕西省茶叶省市共建重点实验室，陕西 安康 725000；2. 安康学院陕西省蚕桑重点实验室，陕西 安康 725000

茶网蝽（Stephanitis chinensis）为半翅目（Hemiptera）异翅亚目（Heteroptera）网蝽科（Tingidae）昆虫，不完全变态，以我国西南地区的油茶和茶树等经济作物为主要寄主。茶网蝽若虫和成虫以刺吸式口器从茶树树冠中、下层成熟叶片背面吸食汁液，其排泄物在叶片上形成黑色黏质物[1]，受害茶树树势衰弱，生长缓慢，发芽率低且芽包小，严重影响茶树生长和茶叶产量[2]。

茶网蝽在陕南地区的茶叶产区一年发生2～3代，繁殖力强，群体数量庞大，种群扩散能力强，新旧茶园均受危害。2010年陕西省汉中市首次报道茶园发生茶网蝽危害[3]，2012年扩散至陕西省安康市，2017年蔓延至汉江流域和嘉陵江流域[4]。茶网蝽成虫体长 3.1～4.2 mm，体小扁平，无趋性，初孵若虫肉眼无法识别，亟待分子生物学鉴定方法[5]。

线粒体基因组（Mitochondrial genome，mitogenome）不易发生重组、倒位、易位等突变，结构保守。线粒体基因进化速率快，为核基因的 2～6倍，被广泛应用于种群多样性鉴定、物种系统发育、物种地理学谱系和系统发生学等方面的研究[6-7]。全世界已知网蝽科昆虫有260属2 124种，但目前NCBI公布的线粒体基因组仅27个（截至2022年3月31号），网蝽科昆虫的线粒体基因组大小差异较大，主要是由于控制区（Control region）数量、长度不等的重复序列造成[8]。本研究利用Illumina和Sanger测序，获得了陕西安康茶区茶网蝽的线粒体基因组全序列及特征，并进行了网蝽科的系统发育分析，旨在为进一步对茶网蝽扩散传播途径、防治及异翅亚目的分类研究提供依据。

1 材料与方法

1.1 试验材料

茶网蝽成虫于2021年5月17日采集自陕西省安康市紫阳县和平茶厂茶园，采集的样本存放在无水乙醇中，置于-20℃冰箱保存备用。

1.2 基因组DNA提取与高通量测序

取乙醇浸泡的茶网蝽单个个体，ddH2O冲洗2遍，然后使用快速DNA提取检测试剂盒[KG203，天根生化科技（北京）有限公司]按照说明书提取基因组 DNA，最后用100 μL ddH2O洗脱。使用 Multiskan FC微量酶标仪（美国Thermo Fisher Scientific公司）对提取的基因组DNA进行检测。选取部分检测合格的基因组DNA委托生工生物工程（上海）股份有限公司进行建库和Illumina Hiseq 2500平台双向高通量测序，测序长度为150 bp，目标数据量4 Gb。

1.3 序列拼接组装与控制区的PCR扩增

对Illumina高通量测序数据进行质控，去除接头和带 N的序列，滑窗法去除得分低于20和长度小于 35 bp的序列。最后数据利用GetOrganelle[9]（K=31、55、85）和NOVOPlasty[10]（K=31）软件进行组装，组装结果使用Bandage[11]进行可视化分析，手动去除可信度低和不确定的部分。在获得的序列缺口处设计引物rrnS-F 和 trnM-R（表1），使用 Gloria Nova HS 2×高保真Taq DNA聚合酶（RK20715，武汉爱博泰克生物科技有限公司）进行 PCR扩增，PCR反应程序：94℃预变性 2 min；94℃变性 30 s，55℃退火 30 s，72℃延伸 5 min，30个循环；72℃延伸10 min。使用1%琼脂糖凝胶电泳检测PCR产物，之后对产物进行 Sanger双向测序，并根据获得序列设计测序引物（表1）进行Gap补齐。

表1 引物序列Table 1 Primers sequence

1.4 基因组注释和分析

使用在线软件MITOS2[12]（http://mitos2.bioinf.uni-leipzig.de/index.py）对茶网蝽线粒体基因组进行注释，遗传密码选择5（Invertebrate）。同时利用blastn和blastp比对网蝽科其他昆虫的线粒体基因组数据确定蛋白质编码基因（Protein-coding genes，PCGs）、tRNA 和 rRNA的边界位置。最后的注释文件提交至NCBI数据库，并利用在线 OGDRAW 软件[13]（https://chlorobox.mpimp-golm.mpg.de/OGDr aw.html）绘制结构图。tRNA的二级结构预测使用在线软件tRNAscan-SE[14]（http://lowelab.ucsc.edu/tRNAscan-SE），未成功预测的二级结构进行手工绘制。

使用在线工具 Tandem Repeats Finder[15]（https://tandem.bu.edu/trf/trf.html）查找控制区中的串联重复序列，使用MISA[16]检测其中的简单重复序列（Simple sequence repeats，SSRs），参数设置为1～6碱基单元的最少重复数分别为10、5、4、3、3、3个。茎环结构预测使用在线工具MFold[17]（http://www.unafold.org/mfold/applications/rna-folding-form.php）。

使用 Mega 11软件[18]统计基因组及各区域 4种碱基的组成和蛋白质编码基因的密码子情况，并通过公式AT skew=(A-T)/(A＋T)和 GC skew=(G-C)/(G＋C)[19]计算核苷酸组成的偏向性。

1.5 系统发育分析

从 NCBI数据库下载异翅亚目昆虫线粒体基因组，分别提取13个PCGs序列，然后利用 macse_v2.05软件[20]基于编码框架模式对序列进行比对，并把对齐的DNA合并成一个序列。利用 raxmlHPC[21]构建极大似然（Maximum likelihood，ML）发育树，GTRGAMMA模式，抽样 1 000次计算bootstrap值，展示高于50%的位点。

2 结果与分析

2.1 茶网蝽线粒体基因组的测序与组装

Illumina测序数据经质控，获得2 025 131条序列，GetOrganelle和NOVOPlasty组装软件获得的线粒体基因组序列完全一致的部分长度为14 561 bp，平均深度为415.04倍。控制区序列的组装得到多种结果，可能存在重复结构。进一步利用引物rrnS-F和trnM-R进行PCR扩增，获得1条4 000 bp左右的条带（图1）。使用设计的多条测序引物对控制区进行Sanger测序，拼接测序结果，获得序列长度为4 040 bp，与Illumina测序组装数据拼接成环，获得茶网蝽线粒体基因组全长为 18 085 bp（GenBank登录号：OM397399）。

图1 rrnS-F和trnM-R引物的PCR扩增Fig. 1 The PCR amplification of rrnS-F and trnM-R

2.2 茶网蝽线粒体基因组注释与结构

分析发现，茶网蝽线粒体基因组为典型的双链闭合环状结构，包含37个基因，分别为13个PCGs、22个tRNA基因和2个rRNA基因（图2，表2）。基因组控制区为3 678 bp，其后为trnI(GAU)、trnQ(UUG)、trnM(CAU)和nad2基因，与昆虫线粒体祖先基因顺序（Ancestral gene order）相同[22]。

图2 茶网蝽线粒体基因组图谱Fig. 2 Map of S. chinensis mitogenome

表2 茶网蝽线粒体基因组结构Table 2 The mitogenome organization of S. chinensis

茶网蝽的线粒体基因组的37个基因共存在18处基因重叠，最长的两处位于nad1和trnL1(UAG)、trnF(GAA)和nad5之间，重叠长度分别为-36 bp和-32 bp，其他位置的重叠长度均在-8 bp以内。基因间隔现象少于基因重叠，37个基因仅存在8处间隔，其中最长两处位于trnV(UAC)和rrnS、trnL1(UAG)和rrnL之间，长度分别为47 bp和25 bp，其他间隔最长8 bp，另有无间隔和重叠的区域10处（表2）。

2.3 茶网蝽线粒体基因组核苷酸组成

茶网蝽线粒体基因组的 AT含量为78.09%，具有明显的AT偏好性（表3）。rRNA基因的AT含量最高，为80.84%，AT富集区为79.99%。PCGs的AT含量最低，为76.76%，其中atp8基因的AT含量最高，达84.28%，其次是nad2和nad4l，分别为 82.21%和80.70%；AT含量最低的是cox1和cob，分别为69.08%和72.21%。核苷酸组成在线粒体基因组不同链间也不对称，多数编码基因所在的J链（Majority strand）的PCGs多为AT偏移和CG偏移，而N链（Minority strand）的PCGs和rRNA基因，为TA偏移和GC偏移（表2）。

表3 茶网蝽线粒体基因组核苷酸组成Table 3 Nucleotide composition of the mitogenome of S. chinensis

2.4 蛋白质编码基因

茶网蝽线粒体基因组的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之间存在相互重叠的现象，atp8和atp6基因间重叠了7 bp，nad4和nad4l基因、nad6和cob基因分别重叠了1 bp和2 bp。

进一步分析发现，茶网蝽线粒体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 茶网蝽线粒体基因组蛋白编码基因相对同义密码子使用度Table 4 Relative synonymous codon usage (RSCU) in the mitogenome of S. chinensis

2.5 tRNA和rRNA基因

茶网蝽线粒体基因组包含22个tRNA基因（图3），长度在63～70 bp。8个tRNA基因（trnC，trnF，trnH，trnL1，trnP，trnQ，trnV，trnY）位于基因组N链，其余14个位于J链（表2）。预测 tRNA基因的二级结构发现，trnS1(GCU)缺少DHU臂，其他tRNA均能形成典型的三叶草结构（图3），该现象存在于多数后生动物中[23-24]。同时，部分tRNA的二级结构还存在G-U和U-U等非经典配对，该现象在其他昆虫类群 tRNA二级结构中也较为常见[25-26]。

图3 茶网蝽线粒体基因组tRNA预测二级结构Fig. 3 Predicted secondary structures of tRNA genes in mitogenome of S. chinensis

茶网蝽线粒体基因组含有 2个 rRNA基因，rrnL位于trnL1(UAG)和trnV(UAC)之间，rrnS位于trnV(UAC)和控制区之间，长度分别为1 237 bp和767 bp，AT含量分别为80.22%和81.88%，具有明显的AT偏向性。

2.6 控制区分析

茶网蝽线粒体基因组的控制区位于rrnS和trnI(GAU)基因之间，长度为 3 678 bp，根据结构将其分为两部分（图4）。第一部分位于序列的前端（1～2 447 bp，front-end region，FER），存在以下结构：（1）3组非串联重复序列。第一组长度为1 540 bp（R1），包含2个498 bp的完整重复单元和 2个该单元的部分序列（477 bp和 63 bp）；第二组长度为 154 bp（R2），包含2个77 bp的完整重复单元；第三组长度为582 bp（R3），由3个194 bp的完整重复单元组成。（2）4处(TTAG)n，其中(TTAG)7和(TTAG)9分别连接非串联重复 R2的 2个单元，另 2处均为(TTAG)5，位于 R2重复的2个单元中。第二部分位于序列的后端（2 448～3 678 bp，back-end region，BER），为串联重复区（R4），包含5个长度为237 bp的完整重复单元和1个该单元的部分序列。同时序列存在5种茎环结构，其中2种5处为单一环结构，另3种11处为多环结构。

图4 茶网蝽线粒体基因组控制区结构Fig. 4 Structure of mitogenome control region of S. chinensis

进一步分析该区域的AT含量为79.99%，高于Pseudacysta perseae、Agramma hupehanum和杨柳网蝽（Metasalis populi）该区的AT含量，低于其他网蝽科昆虫。FER和BER的AT含量分别为79.85%和80.26%（表3），两段区域对不同核苷酸的使用偏好也有所不同，FER的 AT偏度为-0.12，发生 TA偏移，BER的AT使用量相当，有GC偏移倾向。

2.7 系统发育分析

以双线铜翅沫蝉（Aeneolamia contigua）为外群，构建了基于异翅亚目41种昆虫的线粒体PCGs的系统发育树（图5）。结果表明，茶网蝽与直脊冠网蝽（Stephanitis mendica）的亲缘关系最近，其次为杨柳网蝽和Pseudacysta perseae。网蝽科聚为一簇，位于发育树的根部，是臭虫次目（Cimicomorpha）盲蝽科（Miridae）、花蝽科（Anthocoridae）、姬蝽科（Nabidae），以及蝽次目（Pentatomomorpha）、黾蝽次目（Gerromorpha）、歇蝽次目（Nepomorpha）、奇蝽次目（Enicocephalomorpha）的姐妹群。奇蝽次目和黾蝽次目距离最近，与歇蝽次目共同聚合在一簇，蝽次目单独聚合为一个分支。系统发育树支持了歇蝽次目、奇蝽次目、黾蝽次目和蝽次目的单系性，而臭虫次目为并系群。

图5 基于41种昆虫线粒体蛋白质编码基因序列构建的ML系统发育树Fig. 5 Maximum likelihood phylogenetic tree of 41 insects based on the PCG sequences of mitogenome

3 讨论

本研究结合Illumina测序和Sanger测序，获得安康市茶网蝽线粒体基因组全长18 085 bp，是目前 NCBI登录的网蝽科昆虫线粒体全序列（15 208～16 667 bp）中长度最长的。控制区是昆虫线粒体基因组大小变异的主要来源，一般与线粒体基因组的大小呈正相关，控制区的长度差异可以发生在不同物种甚至不同个体之间[27]。目前昆虫中已报道的最短的控制区仅为70 bp，而象鼻虫的控制区长度达 9～13 kb[28]。已报道的网蝽科昆虫线粒体基因组的控制区为733～2 215 bp，而本研究中安康市茶网蝽的控制区为3 678 bp。

控制区是线粒体基因组复制和转录的主要调控区，包含线粒体基因组复制相关信息[29]。黾蝽科圆臀大黾蝽（Aquarius paludum）线粒体基因组控制区存在 TATA motif、poly T stretch、茎环结构等保守结构元件[30]，地长蝽科大黑毛肩长蝽（Neolethaeus assamensis）的线粒体基因组控制区存在poly A结构、茎环结构、串联重复区等[31]。本研究发现，茶网蝽基因组中有多种重复序列和茎环结构，而重庆市茶网蝽[32]仅有 R4重复序列，没有 FER区，对于不同地域茶网蝽基因组之间的差异，有待进一步研究。

本研究茶网蝽的线粒体基因组所有区域均有明显的 AT偏好性，AT含量最高的是rRNA，这与网蝽科昆虫杨柳网蝽、菊方翅网蝽（Corythucha marmorata）、古无孔网蝽（Dictyla platyoma）中的情况相同。蛋白质编码基因的AT含量最低，这有利于物种编码其所需蛋白质的功能[33]。cox1基因是常用的DNA条形码基因，应用于昆虫的识别鉴定[34]。相对较高的 GC含量有利于降低基因的突变率，保持基因稳定性。Yang等[8]研究发现，网蝽科昆虫线粒体cox1基因的进化速率较低，atp8基因最高，本研究中安康市茶网蝽cox1基因的GC含量最高，atp8基因的GC含量最低，与其结论一致。与多数动物一样，茶网蝽线粒体基因组不同链的核苷酸组成也不对称[19]。

网蝽科昆虫线粒体基因组中绝大部分蛋白基因以ATN为起始密码子，使用最多的是ATG，其中cox1和atp6基因的起始密码子均为 ATG，ATA的使用频次略低于 ATT。终止密码子以TAA为主，cox1、atp8、nad6和nad1基因的终止密码子都是TAA；其次为T，cox2基因均以单个T作为终止密码子，cox3和nad5基因使用 T作为终止密码子的现象较多，较少使用 TA或 TAG。昆虫线粒体基因组以 T或TA作为终止密码子现象比较常见，推测其通过转录后 mRNA多腺苷酸化形成完全终止密码子TAA[35]。

臭虫次目是异翅亚目最大的次目，有 17个科超过20 000种昆虫，其系统发育研究目前没有一致的结论。利用形态性状、转录组和核基因等研究发现，臭虫次目昆虫为单系[36-38]，而使用臭虫次目昆虫线粒体基因组数据进行系统发育研究，鉴定其为并系群，尤其是网蝽科昆虫的数据，该现象更显著[8,39-40]。Liu等[39]基于线粒体基因组数据分析异翅亚目的系统发育发现，网蝽科和奇蝽次目、鞭蝽次目、黾蝽次目聚为一簇且接近发育树的根部，是其他异翅亚目昆虫的姐妹群。本研究利用线粒体PCGs进行的系统发育重建结果，能很好地支持奇蝽次目、黾蝽次目、歇蝽次目和蝽次目的单系性。但臭虫次目的情况比较复杂，其姬蝽科和花蝽科聚合在一起，形成与奇蝽次目、黾蝽次目和歇蝽次目并列的分支，而其盲蝽科和网蝽科的位置接近发育树的根部，这与Liu等[39]的研究结果相似。但本研究接近发育树根部的几个分支的支持率均低于50%，这可能与网蝽科昆虫线粒体基因组序列的变异率较大有关[8,39]，较大碱基差异不容易区分系统发育关系，更适合于亲缘关系近的物种的分化检测。