高桥仁蚧线粒体全基因组测序与分析

邓 鋆, 王 刚, 卢聪聪, 张江涛, 黄晓磊,*

(1.福建农林大学植物保护学院, 闽台作物有害生物生态防控国家重点实验室, 福州 350002;2.江西农业大学林学院, 鄱阳湖流域森林生态系统保护与修复国家林业和草原局重点实验室, 南昌 330045)

线粒体是真核生物细胞中具有独立遗传物质的细胞器，是细胞制造能量的场所，参与细胞分化、细胞信息传递和细胞凋亡等过程(潘宝平和卜文俊, 2005; Cameron, 2014)。截至2020年，在GenBank数据库中已发布489种半翅目(Hemiptera)昆虫完整线粒体基因组。作为半翅目一员，蚧总科(Coccoidea)现已知35科8 310种蚧虫(García Moralesetal., 2016)，但在GenBank数据库中仅有3个蚧虫完整的线粒体基因组序列，均为蚧科(Coccidae)物种(Dengetal., 2019; Luetal., 2020; Xuetal., 2021)。相比现有的蚧虫物种数量，蚧虫线粒体基因组的测序工作明显滞后，GenBank数据库中严重缺少除cox1基因外的其他蚧虫线粒体标记基因数据，这阻碍了分子进化、系统发育和生物地理学等相关工作在蚧虫中的开展，也阻碍了线粒体基因组在半翅目昆虫研究中进一步的应用。Deng等(2019)公布首个蚧虫即日本龟蜡蚧Ceroplastesjaponicus完整的线粒体基因组序列，并于2020年完整注释了蚧科咖啡黑盔蚧Saissetiacoffeae线粒体基因组，弥补了GenBank数据库因长期缺失蚧虫线粒体基因组数据而造成的空白(Luetal., 2020)。研究发现蚧虫线粒体基因组中存在高度的基因重排现象、高AT含量以及非典型的tRNAs二级结构(例如缩减tRNAs，主要表现在缺失DHU臂或TΨC臂)(Luetal., 2020)。事实上，GenBank数据库缺失的不仅仅是完整的蚧虫线粒体基因组序列，数据库中只有两条超过5 000 bp的部分线粒体基因序列。这些现象反映想要获得完整的蚧虫线粒体基因组可能具有较高的难度。

高桥仁蚧Aclerdatakahashii隶属于半翅目蚧总科仁蚧科(Aclerdidae)，分布在巴西、中国、埃及、印度、印度尼西亚、马来西亚、美国等14个国家和地区，涵盖古北界、新北界、东洋界、新热带界(García Moralesetal., 2016)。高桥仁蚧各龄期均生活在植物的叶鞘下，寄主主要为禾本科芒属Miscanthus、甘蔗属Saccharum、粽叶芦属Thysanolaena植物(García Moralesetal., 2016)。本研究对高桥仁蚧粒体基因组进行二代高通量测序，获得其完整的线粒体基因组序列，对其碱基组成、密码子使用情况、tRNA结构进行了预测分析，并基于线粒体基因组13个蛋白编码基因构建半翅目系统发育树。研究结果将为蚧总科分类及仁蚧科线粒体基因组研究提供科学依据。

1 材料与方法

1.1 昆虫样本采集与DNA的提取

本研究所用的高桥仁蚧成虫样本于2019年8月5日采集自福建省建瓯市云际山公园内(27°1′19″N, 118°18′36″E)，利用《中国动物志》(王子清, 2001)进行形态鉴定。测序样本置于95%的酒精中-80℃保存。使用DNeasy DNA Extraction Kit试剂盒(Qiagen, Hilden, 德国)提取10头成虫混合样本的总DNA，经NanoDrop 2000分光光度计和1%的琼脂糖凝胶电泳检测总DNA的质量和浓度。

1.2 高通量测序及线粒体基因组的拼接

将1.1节质量合格浓度大于2 μg的DNA样本送至北京诺禾致源科技股份有限公司进行高通量测序，构建350 bp的短片段文库，利用Illumina HiSeq测序平台进行双末端(paired-end, PE)测序。 将下机的原始数据去接头和质控。利用NOVOPlasty 3.7.1(Dierckxsensetal., 2017)软件以该物种线粒体cox1基因片段为种子序列进行线粒体基因组的拼接，在K-mer=23, 27, 33下均获一条长为16 599 bp的成环序列；为保证拼接的准确性，使用MEGAHIT 1.0(Lietal., 2016)对clean data进行从头组装，建立本地Blast数据库，提取一条长为16 735 bp的线粒体基因组序列；随后将两种软件所拼接的序列进行比对验证，未发现有碱基不同。最终使用Pilon 1.3.2(Walkeretal., 2014)对16 599 bp拼接序列进行校正，获得高桥仁蚧线粒体基因组全长序列。

1.3 线粒体基因组的注释与序列分析

利用MITOS2在线服务器(Berntetal., 2013)对获得的线粒体基因组序列进行初步注释。剩余未发现的tRNA由ARWEN 1.2注释(Laslett and Canbäck, 2008)。使用TRF(Tandem Repeats Finder)在线服务器(https:∥tandem.bu.edu/trf/trf.html)查找串联重复序列(Benson, 1999)。运用Phylosuite 1.2.1(Zhangetal., 2020)分析核苷酸密码子使用频率和核酸组成。利用公式AT-skew=(A-T)/(A+T)和GC-skew=(G-C)/(G+C)分别计算AT偏斜和GC偏斜，使用OGDRAM在线可视化软件绘制高桥仁蚧线粒体全基因组(Greineretal., 2019)。

1.4 系统发育分析

从GenBank数据库中下载半翅目15科31种已报道的线粒体全基因组序列，基于13个蛋白编码基因，在Phylosuite软件中使用MACSE 2.03(Ranwezetal., 2018)进行序列多重比对，并进行串联，获得蛋白编码基因核苷酸序列数据集。以蜚蠊目维州散白蚁Reticulitermesvirginicus线粒体全基因组为外群，应用IQ-TREE1.6.8(Nguyenetal., 2015)和MrBayes 3.2.6(Ronquistetal., 2012)分别构建最大似然法(maximum likelihood, ML)和贝叶斯推断法(Bayesian inference, BI)系统发育树。IQ-TREE系统树分支节点的置信水平用ultrafast bootstrap approximation approach (UFBoot)(重复抽样10 000次)估计。利用ModelFinder(Kalyaanamoorthyetal., 2017)筛选出最佳核苷酸替换模型为GTR+I+G+F，MrBayes在GTR+I+G+F模型进行贝叶斯推断，运算2 000 000代，每1 000代抽样一次，burn-in参数设置为25%。

2 结果

2.1 高桥仁蚧线粒体基因组结构和碱基组成

高桥仁蚧线粒体基因组序列全长16 599 bp(GenBank登录号: MW839575)，由一段控制区和37个基因(13个蛋白编码基因、22个tRNA基因和2个rRNA基因)组成，呈闭合双链环状的DNA分子结构(图1)。通常将多数编码基因所在的链定义为主要编码链J链(majority strand)，则另外一条链定义为次要编码链N链(minority strand)。线粒体基因组中9个蛋白编码基因(cox1,cox2,cox3,cob,atp8,atp6,nad2,nad3和nad6)和14个tRNA基因(trnM,trnW,trnL2,trnK,trnD,trnG,trnR,trnS1,trnE,trnN,trnA,trnI,trnT和trnS2)在J链上，剩余14个编码基因在N链上。线粒体基因组中A, T, C和G碱基含量分别为47.54%, 36.97%, 11.05%和4.44%。总A+T含量为84.51%，AT偏斜为0.125，呈现明显的AT偏向性(表1)。线粒体基因组中共有9处基因重叠(表2)，共76 bp，其中最长的重叠区域位于trnE-trnF之间，为37 bp；其次位于trnI-nad2之间，长为15 bp。存在基因间隔区24个，共1 250 bp；最长间隔区在trnT-cob处，长度为184 bp。最短间隔区位于trnN-trnQ处长度为2 bp。无间隔无重复区有4处(表2)。

图1 高桥仁蚧线粒体基因组结构Fig.1 Structure of the mitochondrial genomeof Aclerda takahashii

表1 高桥仁蚧线粒体基因组核苷酸组成Table 1 Nucleotide composition of the mitochondrial genome of Aclerda takahashi

续表1 Table 1 continued

表2 高桥仁蚧线粒体基因组结构Table 2 Organization of the mitochondrial genome of Aclerda takahashii

2.2 高桥仁蚧线粒体蛋白质编码基因和相对同义密码子使用频率(RSCU)

高桥仁蚧线粒体基因组中13个蛋白编码基因序列总长为10 632 bp，约占线粒体基因组全序列的64.05%。nad5的序列最长，为1 620 bp。atp8的序列最短，长141 bp；均采用标准的起始密码子(ATN)和终止密码子(TAA),其中，3个蛋白编码基因(cox1,atp8和atp6)的起始密码子为ATA， 5个蛋白编码基因(cox2,cox3,cob,nad1和nad4)的起始密码子为ATG，其余5个蛋白编码基因的起始密码子为ATT(表2)。共有62个不同的密码子，其中有5个密码子使用次数最多：AUA(Met)439次，AUU(Ile)408次，UUU(Phe)390次，UUA(Leu2)366次和AAU(Asn)260次。氨基酸占比由高到低的排序为Met(13.62%)，Ile(12.89%)，Phe(12.21%)，Leu2(11.95%)和Asn(9.23%)(图2)。

图2 高桥仁蚧线粒体基因组的相对同义密码子使用频率(RSCU)Fig.2 Relative synonymous codon usage (RSCU) in the mitochondrial genome of Aclerda takahashii

2.3 高桥仁蚧线粒体基因组tRNA基因和rRNA基因

22个tRNA基因中，MITOS2成功注释了15个，剩余7个由ARWEN注释；长度变化范围在48(trnM)～79 bp(trnE)之间；只有10个tRNA基因(trnR,trnN,trnC,trnE,trnI,trnL1,trnL2,trnK,trnF和trnM)的二级结构是典型的三叶草结构，其余12个皆为非典型的结构：trnY和trnW缺失DHU臂，trnS1和trnS2的DHU臂和TΨC臂缺失，剩下8个tRNA基因都缺失TΨC臂。此外在一些tRNA基因的氨基酸接受臂上出现碱基错配现象，例如trnR(U-U),trnE(U-U),trnN(A-A)以及trnF(U-C)。除了trnN的反密码子臂上出现1处碱基错配现象(C-A)，其余都能形成完全配对的臂(图3)。

图3 高桥仁蚧线粒体基因组tRNA基因二级结构Fig.3 Secondary structure of tRNA genes in the mitochondrial genome of Aclerda takahashii

2.4 系统发育

运用最大似然法和贝叶斯推断法构建的系统发育树均展示了相同的拓扑结构(图4)，并且每个分支节点上都有高的支持率(自举值>65，后验概率>0.90)。二者结果充分展示了这些总科的如下系统发育关系：(((跳蝽总科(Saldoidea)+细蝽总科(Leptopodoidea))+(臭虫总科(Cimicoidea)+猎蝽总科(Reduvioidea)))+(角蝉总科(Membracoidea)+(蜡蝉总科(Fulgoroidea)+(木虱总科(Psylloidea)+(粉虱总科(Aleyrodoidea)+(蚜总科(Aphidoidea)+蚧总科(Coccoidea)))))))。两个系统发育树展示蚧科的日本龟蜡蚧C.japonicus、咖啡黑盔蚧S.coffeae、朝鲜球坚蚧Didesmococcuskoreanus均分布在一个分支上，高桥仁蚧作为仁蚧科代表，单独形成一个分支。

图4 最大似然法构建的基于线粒体基因组的13个蛋白质编码基因核苷酸序列的半翅目系统发育树Fig.4 Phylogenetic tree of Hemiptera based on the sequences of 13 protein-coding genesof mitochondrial genome using maximum likelihood method节点旁左侧数字代表Bootstrap支持率，右侧数字代表贝叶斯后验概率。Numbers on nodes refer to bootstrap support values (left) and Bayesian posterior probabilities (right).

3 讨论

高桥仁蚧tRNA二级结构具有特殊性，一部分tRNA无法通过MITOS2得到注释(表2)。主要是由于这些tRNA是非典型的三叶草结构，缺失DHU臂或TΨC臂，其中tRNAser(S1)和tRNAser(S2)更是同时缺失DHU臂和TΨC臂(图3)。同时在部分tRNA的氨基酸臂上，有个别碱基不配对，这些原因导致蚧虫的tRNA无法通过线粒体注释软件MITOS2准确注释。一般情况下，动物的trnS1(AGN)通常会缺失DHU臂，而不具有经典的三叶草结构(Boore, 1999)。这种trnS1缺失DHU臂的现象在半翅目中均有发现(郭仲龙和袁明龙, 2016)，白背飞虱Sogatellafurcifera的线粒体基因组的2个trnS基因和1个trnG基因也有缺臂的现象存在(Zhangetal., 2014)。在高桥仁蚧基因组中存在严重的tRNA缺臂现象，22个tRNA中的12个均存在缺臂现象(图3)，这一现象在蚧科的线粒体基因组中也有发现(Luetal., 2020)。然而在其他类群的昆虫中，大量tRNA的缺臂现象鲜有发生(Luetal., 2020)。类似的报道主要来自于线虫(Jühlingetal., 2012; Lorenzetal., 2017)、螨类(Xueetal., 2016; 赵亚男和李朝品, 2020)和蛛形纲动物(Masta and Boore, 2008)。在后生动物中，线粒体tRNA通过编译可以修复这些缺失的臂，形成典型的三叶草结构的tRNA(Lavrovetal., 2000; Segoviaetal., 2011)。这种转录后的编译可能广泛存在于仁蚧科的类群中。同时，我们也发现高桥仁蚧线粒体蛋白编码基因具有很强的密码子偏好性(图2)，这与半翅目其他昆虫的密码子使用情况类似(Luetal., 2020)。

半翅目昆虫的AT含量一般在70%～80%之间，最低AT含量的半翅目线粒基因组为粉虱科Bemisiaafer(65.67%)，最高的是粉虱科Aleurodicusdugesii(86.33%)(Wangetal., 2015)。其中胸喙亚目(Sternorrhyncha)是半翅目4个亚目中AT含量最高的亚目，平均AT含量达到80.97%(郭仲龙和袁明龙, 2016)。目前GenBank中最高AT含量(88.02%)的昆虫线粒体基因组为膜翅目锤角细蜂科Trichopriadrosophilae线粒体基因组(GenBank登录号: NC_048491)。作为胸喙亚目的一员，3个已知的蚧虫线粒体基因组AT含量均超过了80.00%(Dengetal., 2019; Luetal., 2020)，高桥仁蚧的AT含量达到84.51%(表1)。超高的AT含量造成测序上的困难，AT重复片段以及poly (A)的特殊结构经常导致一代测序的质量不佳，出现套峰或测序信号的衰减中断。在高桥仁蚧的线粒体基因组中，一个非常极端的例子是在nad5基因的位置上，连续出现了31个A碱基，这或可以解释为何没有一个蚧虫线粒体基因组是通过一代测序获得的原因。二代测序技术，也就是高通量测序技术，本身由于对于GC有偏好性，导致在AT含量高的地方往往测序覆盖度不一致，甚至缺失覆盖(Chenetal., 2013; Browneetal., 2020)。这导致利用二代测序技术拼接蚧虫线粒体基因组时，获得的线粒体基因组常碎片化。再加上由于蚧虫本身个体小，采集困难，易被寄生蜂寄生(Dengetal., 2012)，这些原因都导致很难获得完整蚧虫基因组，也成为GenBank数据库中严重缺少蚧虫线粒体基因组数据的主要原因。

在构建的半翅目代表物种系统发育树上(图4)，仁蚧科与蚧科物种亲缘关系最近，这与之前基于18S rRNA, 28S rRNA和EF-1α的3个核基因片段的蚧总科系统发育研究的结果(Vea and Grimaldi, 2016)一致。由于蚧虫线粒体基因组数据有限，在蚧虫中基于线粒体多基因的系统发育研究还有待开展。已知的蚧虫线粒体基因组构建的系统发育树也无法真实反映蚧虫的进化的格局和系统发育关系。半翅目蚧总科中线粒体基因组数据的缺乏，严重阻碍了对蚧虫本身以及整个半翅目昆虫利用线粒体基因组开展进化、系统发育和生物地理等问题的研究，获得更多不同科的蚧虫线粒体基因组序列将成为解决这一问题的关键。

本研究获得了高桥仁蚧线粒体基因组全序列，是仁蚧科首个线粒体基因组序列，提供了仁蚧科线粒体基因结构和组成的基础数据，为后续获得仁蚧科更多物种的线粒体基因组提供了参考。