Atteva charopis 的比较线粒体基因组研究及系统发育分析

林兴雨 宋南

（河南农业大学植物保护学院，郑州 450046）

巢蛾总科（Yponomeutoidea）是鳞翅目（Lepidop‑tera）早期衍生的一个类群，全世界范围内共分布11个科约1 800 个种，其中温带地区分布较多［1-2］。近年来，由于巢蛾总科是鳞翅目中最早进化取食草本和灌木等植物的一个类群，常被用于鳞翅目昆虫与植物互作的相关研究［3-4］。巢蛾总科中的某些种群会严重危害农业生产，属于重要的农业害虫。例如蔬菜上危害较为严重的小菜蛾（Plutella xylostella），其爆发期会对农业生产造成不可估量的经济损失［5］。此外，巢蛾总科中一些物种的幼虫以当地植被甚至景观植物为食，导致树木提前落叶甚至衰亡［6-7］。然而，目前对于巢蛾总科的系统学的研究相对较少，对其生物多样性了解不足。因此，对巢蛾总科进行深入的系统发育研究具有重要的科学意义。

随着高通量测序技术的不断发展与成熟，越来越多的昆虫线粒体基因组被成功测序并获得。迄今为止，共公布了约19 万条昆虫部分或完整的线粒体基因组数据。线粒体基因组是一条环状、双链、闭合的DNA 分子，平均长度约为15 000-18 000 bp，通常由13 个蛋白质编码基因（13 PCGs）、24 个RNA 编码基因（rrnL、rrnS和22 个tRNA 基因）和一段非编码控制区组成［8］。线粒体基因组数据具有母系遗传、基因组组织结构保守以及进化速率快等独特优势，常作为一种分子遗传标记广泛应用于昆虫分子系统学、系统地理学和种群遗传学等研究［9-12］。截至2023 年2 月，GenBank 中仅注释和释放了1 个Atteva属的完整线粒体基因组数据。有限的线粒体基因组数据在一定程度上制约了从线粒体基因组角度深入研究巢蛾总科的系统发育关系。

Atteva charopis在1903 年被Turner 首次命名（https://www.biodiversitylibrary.org/item/30143#page/96/mode/2up）。A. charopis成虫的前翅呈深橙色，分布有白色斑点；后翅为黄色，在翅尖的部位有暗到深棕色的阴影；翅展大约2.5 cm。目前暂未有对A.charopis生物学习性和系统发育等的研究报道。本研究利用高通量测序技术首次获得并注释了巢蛾总科中斑巢蛾科（Attevidae）的一个物种A. charopis完整的线粒体基因组数据，分析预测了线粒体基因组核苷酸碱基组成、密码子使用频率的情况以及转运RNA 的二级结构等，同时详细描述了巢蛾总科的蛋白质编码基因的进化速率、核苷酸多样性、线粒体基因重排等。此外，以GenBank 中已经释放的巢蛾总科中5 个科的10 个物种的线粒体基因组数据（Acrolepiopsis assectella、Prays oleae、Atteva aurea、Scythropia crataegella、Leucoptera malifolie、Lyonetia clerkella、Plutella armoraciae、Plutella porrectella、Plutella australiana和Plutella xylostella）为内群，以2 个谷蛾总科（Tineoidea）（Phyllocnistis citrella和Caloptilia theivora）物种的线粒体基因组数据为外群进行系统发育分析（表1），并分别利用线粒体基因组中13 个蛋白质编码基因的核苷酸序列和氨基酸序列重建巢蛾总科内部系统发育关系以及A. charopis的系统发育位置。这些结果为深入了解巢蛾总科和A. charopis的系统发育提供了宝贵的线粒体基因组数据。

表1 本研究中用于构建系统进化树的物种样本信息Table 1 Species sample information used in the phylogenetic tree construction in this study

1 材料与方法

1.1 材料

2022 年8 月，从河南省平顶山市（33°47'N,112°19'E）采集A. charopis成虫标本。随后将采集后的标本浸泡在100%无水乙醇中并保存在河南农业大学植物保护学院标本馆的‑20℃超低温冰箱内备用。

1.2 方法

1.2.1 DNA 提取 根据动物基因组DNA 提取试剂盒（TIANamp Genomic DNA Kit, China）中的操作步骤，提取A. charopis成虫胸部肌肉组织的总DNA。使用NanoDrop 2000 分光光度计和1.5%的琼脂糖凝胶电泳检测提取的总DNA 的质量和浓度。

1.2.2 高通量测序与线粒体基因组组装 首先将检测质量合格且浓度大于20 ng/μL 的DNA 送至河南微智因科技股份有限公司构建350 bp 的小片段文库；然后使用MGI 2000 测序平台对所构建的小片段测序文库进行双末端（2 × 150 bp PE）测序；最后将下机的原始数据去除接头，删除低质量的reads，并利用NGS QC Toolkit v2.3.3 软件的默认参数进行质控［13］，获得clean reads。接着利用GetOrganelle v1.7.5.2 软件［14］中‑F animal_mt 参数对质控后的clean reads 进行线粒体基因组的组装和拼接，最终获得高质量线粒体基因组全序列。

1.2.3 线粒体基因组的注释和分析 将组装得到的A. charopis的线粒体基因组数据在MITOS［15］在线网站上进行初步的功能注释，具体参数设置为：Genetic Code：05‑inverterbrate；Reference：RefSeq 89Metazoa；其余参数设置均为默认。通过MITOS 预测22 个转运RNA 基因的二级结构。根据对近缘种线粒体基因组的序列比对验证后，手动矫正13 个蛋白质编码基因、rrnL和rrnS核糖体基因和22 个转运RNA 的基因边界，获得注释好的线粒体基因组序列。然后，将注释好的线粒体基因组序列上传至NCBI并获得GenBank 登陆号（OQ130005）。利用MEGA 11.0［16］软件分别进行分析和计算A. charopis线粒体基因组所有的蛋白质编码基因的碱基组成以及密码子的使用频率情况，并利用Kimura‑2‑Parameter 模型计算A. charopis与巢蛾总科的其他10 个物种之间的种间遗传距离。此外，根据AT‑skew =（A-T）/（A+T）和GC‑skew =（G-C）/（G+C）公式分别计算位于正负链上的蛋白质编码基因、核糖体RNA 基因、转运RNA 基因以及整个线粒体基因组的AT 偏倚和GC 偏倚［17］。使用DnaSP v5.10.01［18］计算巢蛾总科的蛋白质编码基因的核苷酸多样性（nucleotide diversity）以及非同义进化率non‑synonymous（Ka）和同义进化率synonymous（Ks）的比值。

1.2.4 多序列比对 利用MAFFT 软件中的E‑INS‑I策略分别对蛋白质编码基因的核苷酸和氨基酸序列进行序列比对［19］。使用trimAl v1.4 对序列中比对质量差的序列片段进行修剪、删除和对齐［20］。最后，使用FASconCAT‑G_v1.04 进行串联得到建树所需的数据矩阵［21］。本研究用于构建系统发育分析的数据矩阵包括13 个蛋白质编码基因的核苷酸序列以及氨基酸序列构建的数据矩阵。

1.2.5 系统发育分析 采用IQ‑TREE 2.0.6［22］中最大似然法（maximum likelihood）对13 个蛋白质编码基因的核苷酸数据矩阵进行系统发育分析，其中根据贝叶斯信息算法（BIC）筛选的最适系统发育模型为GTR+F+I+G4，另外使用自举检验置信度评估系统发育树的节点支持率（运行次数设置为10 000）。接着运用MrBayes v3.2 中的贝叶斯法对13 个蛋白质编码基因的氨基酸数据矩阵做系统发育分析［23］。其中设置4 条蒙特卡罗马尔柯夫链（Markov Chain Monte Carlo, MCMC），同时运行1×108代，每运行完1 000 代随机进行抽样一次，丢弃25%的老化树。最后获得多数一致树，使用后验概率值评估分枝支持度。

2 结果

2.1 线粒体基因组分析

采用OGDRAW Version 1.1［24］在线绘制A.charopis的线粒体基因组图（图1）。A. charopis的线粒体基因序列全长为15 163 bp，其中A、T、G 和C含量分别为39.60%、41.55%、7.58%和11.27%。该线粒体基因组由37 个基因（13 个蛋白质编码基因、2 个核糖体RNA 基因和22 个转运RNA 基因）和1个非编码控制区组成（表2）。

图1 A. charopis 线粒体基因组结构Fig. 1 Structure of A. charopis mitochondrial genome

表2 A. charopis 线粒体基因组注释Table 2 Annotation of A. charopis mitochondrial genome

在A. charopis的线粒体基因组结构中，4 个蛋白质编码基因和8 个转运RNA 基因以及2 个核糖体RNA 基因位于L 链（light strand）；剩余的9 个蛋白质编码基因和14 个转运RNA 基因位于H 链（heavy strand）。非编码控制区位于rrnS和trnI之间，长度为426 bp，其A、T、G 和C 含量分别为44.37%、52.11%、1.41%和2.11%。A. charopis位于正链的蛋白质编码基因、负链的蛋白质编码基因、负链的核糖体RNA 基因、正链的转运RNA 基因、负链的转运RNA 基因和整个线粒体基因组均具有较高的AT含量，AT 含量在78.3%-85.1%范围内，具有较高的AT 含量。其中位于负链的核糖体RNA 基因最高为85.1%，位于正链的蛋白质编码基因最低为78.3%。AT 偏倚和GC 偏倚分别在‑0.171-0.04 和‑0.196-0.425范围内（表3）。

表3 A. charopis 的核苷酸组成和偏倚Table 3 Nucleotide composition and skewness of the A.charopis mitochondrial genome

2.2 蛋白质编码基因分析

A. charopis的13 个蛋白质编码基因的起始密码子中，除cox1利用CGA 开头外，其余12 个蛋白质编码基因都是以ATN（ATT、ATA 和ATG）作为起始密码子。4 个蛋白质编码基因（cox2、nad5、nad4和nad6）使用不完整的终止密码子T 或TA 结尾，剩余9 个蛋白质编码基因都是以完整的终止密码子TAA 或TAG 结尾。在所有蛋白质编码基因的氨基酸使用中，苏氨酸（Thr）使用频率最高为41.54%，其次为丙氨酸（Ala）和半胱氨酸（Cys），使用频率分别为39.59%和11.27%，甘氨酸（Gly）使用频率最低，为7.58%，其余氨基酸均未使用。A. charopis线粒体基因组的相对密码子使用频率见表4，密码子使用次数最多的是UUU 共使用444 次，密码子GCG 和GGG 使用次数最少，各使用了一次。

表4 A. charopis 线粒体基因组的相对密码子使用频率Table 4 Relative synonymous codon usage of the mitochondrial genome of A. charopis

巢蛾总科的核苷酸多样性和进化率如表5 所示，其中13 个蛋白质编码基因的核苷酸多样性的范围在0.127（cox2）-0.247（nad6）之间。nad6基因（Pi=0.247）在所有蛋白质编码基因中具有最多样化的核苷酸，其次是atp8基因（Pi= 0.205）和nad2基因（Pi= 0.204）。然而，cox2基因（Pi= 0.127）、cox1基因（Pi= 0.129）和cob基因（Pi= 0.159）的核苷酸多样性值相对较低，是相对保守的基因。

表5 巢蛾总科线粒体基因组13 个蛋白质编码基因的核苷酸多样性和进化率分析Table 5 Nucleotide diversity and evolutionary rate analyses of 13 protein-coding genes of the mitochondrial genome of Yponomeutoidea

巢蛾总科的13 个蛋白质编码基因的进化率显示，atp8基因（0.640）、nad6基因（0.500）和nad4基因（0.496）进化率相对较高，而cox1基因（0.133）、cox2基因（0.186）和cob基因（0.207）的进化率相对较低。

2.3 tRNA和rRNA基因分析

A. charopis线粒体基因组的22 个转运RNA 基因的序列长度在60-71 bp 之间。其中trnS1基因最短，为60 bp，trnK基因最长，为71 bp。22 个转运RNA 基因中，H 链上的RNA 基因较多，有14 个转运RNA 基因位于H 链，L 链上有8 个转运RNA 基因。rrnL和rrnS基因全部位于L 链。其中trnS1基因缺少二氢尿嘧啶（DHU）臂，不能形成完整的三叶草结构，只能形成一个简单的环。而其他21 个转运RNA 基因都能形成完整的三叶草结构（图2）。此外，trnS1的反密码子不是常见的GCU，而是UCG。A. charopis线粒体基因组的rrnL的全长为1 330 bp，A+T 含量为84.74%。rrnS的全长为697 bp，A+T 含量为85.79%。综上可知，两个核糖体基因均具有较高的AT 含量。

图2 A. charopis 22 个 tRNA 基因的二级结构Fig. 2 Secondary structure of 22 tRNAs genes of A. charopis

2.4 巢蛾总科线粒体基因重排分析

在现有的巢蛾总科的线粒体基因组数据中，各个物种均出现了不同程度的基因重排事件（图3）。其中潜叶蛾科（Lyonetiidae）的L. malifoliella的线粒体基因组分别出现了基因转位和颠倒事件，巢蛾总科的其他物种的线粒体基因组数据均出现了不同程度的基因转位事件。例如，潜叶蛾科的L. clerkella与Praydidae（Prays oleae）、 菜蛾科（Plutellidae）（Acrolepiopsis assectella、Plutella armoraciae、Plutella porrectella、Plutella australiana和Plutella xylostella）都出现了nad2‑rrnS与trnM基因之间出现基因转位。然而，在斑巢蛾科中A. aurea的trnM与trnQ基因之间以及tRNA与nad3基因之间分别出现了基因转位。本研究测序获得的A. charopis物种的线粒体基因组数据中只有trnM与trnQ基因之间出现了转位。

图3 巢蛾总科和模式昆虫的线粒体基因排列Fig. 3 Gene order of Yponomeutoidea mitogenome and putative insects

2.5 系统发育分析

基于最大似然法和贝叶斯法这两种推断方法的系统发育分析结果分别如图4 和图5 所示。两种系统发育分析方法构建的树拓扑结构都支持潜叶蛾科、Praydidae、斑巢蛾科为单系群，菜蛾科为非单系群。两种系统发育分析方法强烈支持斑巢蛾科为单系群且有较高的节点支持值（BS = 100，PP = 1）。13 个蛋白质编码基因的核苷酸系统发育分析结果显示，潜叶蛾科为单系群，节点支持值相对低，BS=83。然而，13 个蛋白质编码基因的氨基酸系统发育分析结果强烈支持潜叶蛾科为单系群，节点支持值PP=1。

图4 基于13 个蛋白编码基因的核苷酸序列构建的最大似然树Fig. 3 Maximum likelihood tree inferred from the nucleotide sequences of 13 protein-coding genes

图5 基于13 个蛋白编码基因的氨基酸序列构建的贝叶斯树Fig. 5 Bayesian tree inferred from the amino acid sequences of 13 protein-coding genes

此外，基于最大似然法和贝叶斯法的两种系统发育分析结果都强烈支持A. charopis与A. aurea（BS=100，PP=1）的姐妹群关系。

巢蛾总科的其他10 个物种（潜叶蛾科的L.malifoliella和L. clerkella，Praydidae 的P. oleae，斑巢蛾科的A. aurea，菜蛾科的A. assectella、P.armoraciae、P. porrectella、P. australiana和P.xylostella以及Scythropiidae 的S. crataegella）之间的种间遗传距离的结果显示，A. charopis与A. aurea的遗传距离最近为0.122。然而，A. charopi与L.malifoliella遗传距离结果显示最远为0.235（表6）。种间遗传距离与两种系统发育结果所显示的亲缘关系一致支持A. charopis与A. aurea有较近的亲缘关系。

表6 巢蛾总科昆虫的线粒体蛋白质编码基因基于Kimura-2-Parameter 模型的遗传距离Table 6 Genetic distances of mitochondrial protein-coding gene sequences in Yponomeutoidea based on Kimura-2-Parameters

3 讨论

本文利用高通量测序技术首次获得了A. charopis的线粒体基因组全序列，分别分析和预测了它的线粒体基因组的特点以及22 个转运RNA 基因的二级结构。此外，基于GenBank 释放的巢蛾总科的10条线粒体基因组数据，利用13 个蛋白质编码基因分析了巢蛾总科的核苷酸多样性和进化速率，这些结果有助于我们更好地理解它们的系统发育关系。A.charopis的线粒体基因组长度为15 163 bp，而多数线粒体基因组长度范围在15-18 kb 之间［8］，这表明A. charopis的线粒体基因组在多数线粒体基因组长度范围内。巢蛾总科的线粒体基因组的基因排列顺序与模式昆虫不一致［25］，它们出现几种不同的基因重排事件。

计算巢蛾总科的13 个蛋白质编码基因的非同义进化率与同一进化率的比值结果显示，atp8基因的进化率明显高于其他蛋白质编码基因，进化速率高达0.640。然而，cox1基因的进化速率相较于其他蛋白质编码基因进化速率最慢，为0.133。在大多数昆虫中的cox1基因都有相对较低的进化率［26］。此外，Xiao 等［27］的研究结果认为不仅是昆虫，地球上几乎所有动物的cox1基因的非同义进化率与同一进化率的比值结果都相对较低。这可能是由于cox1基因经过了最高的纯化选择造成的［28］。

迄今为止，关于巢蛾总科的系统发育研究的相关报道较少。Sohn 等［1］首次用8 个核基因片段构建了巢蛾总科的系统发育关系。系统发育结果支持巢蛾总科内部的10 个科（巢蛾科Yponomeutidae、茎蛾科Ypsolophidae、菜蛾科Plutellidae、雕蛾科Glyphipterigidae、 银蛾科Argyresthiidae、 斑巢蛾科Attevidae、Praydidae、 举肢蛾科Heliodinidae、Bedelliidae 和Scythropiidae）为单系群。此外，Jeong等［29］基于巢蛾总科的13 个蛋白质编码基因和2 个核糖体基因（rrnL和rrnS）构建了数据矩阵，并利用最大似然法重建了系统发育关系，结果显示菜蛾科、茎蛾科、巢蛾科、斑巢蛾科、Praydidae 和潜叶蛾科Lyonetiidae 为单系群。然而，斑巢蛾科与Praydidae 有较近的亲缘关系但节点支持值较低，BS=64。

本文基于13 个蛋白质编码基因的核苷酸序列利用最大似然法，13 个蛋白质编码基因的氨基酸序列利用贝叶斯法分别重建了巢蛾总科的系统发育关系。两种系统发育推断方法与Sohn 等［1］和Jeong 等［29］结果一致，均支持潜叶蛾科、Praydidae 和斑巢蛾科为单系群。然而，在最大似然法中，菜蛾科的A.assectella物种与Scythropiidae 的S. crataegella物种形成一枝以及贝叶斯法中A. assectella物种单独形成一枝，都导致系统发育结果显示菜蛾科为非单系群。此外，由于Scythropiidae 科的线粒体基因组仅有（S.crataegella）一条被释放，对于它的系统发育关系仍然是不确定的。利用核苷酸数据矩阵在最大似然法分析中，斑巢蛾科和Praydidae 为姐妹群但有相对低的节点支持值BS=87，这一结果与Jeong 等［29］的研究结果一致。然而，利用氨基酸数据矩阵在贝叶斯法系统发育分析强烈支持斑巢蛾科和Praydidae 形成姐妹群（PP=1）。

两种系统发育推断方法与种间遗传距离的结果都支持A. charopis与A. aurea有相对较近的亲缘关系。这可能是由于两个物种的形态特征以及生活习性上有较多的相同之处，促使线粒体基因间更多的基因相同形成的。

4 结论

本研究首次获得了A. charopis的线粒体全基因组数据，分别预测和分析了A. charopis的转运RNA的二级结构和线粒体基因组的结构等特点，比较分析了巢蛾总科的线粒体基因重排，计算了巢蛾总科13 个蛋白质编码基因的核苷酸多样性和进化率。两种系统发育分析结果显示潜叶蛾科、Praydidae 和斑巢蛾科为单系群，斑巢蛾科和Praydidae 为姐妹群。确定了A. charopis与A. aurea亲缘关系较近。本研究为后续巢蛾总科的系统发育、谱系地理学以及种群遗传学等相关研究提供线粒体基因组数据。