眼优角蚱两性成体转录组差异分析

覃英灿，容万韬，李月妹，陆桂红，李晓东

(河池学院化学与生物工程学院，广西宜州 546300)

蚱类昆虫在分类上属于直翅目Orthoptera蝗亚目Caelifera蚱总科Tetrigoidea，由于其背板呈菱形，旧称菱蝗。蚱类昆虫为不完全变态昆虫，在整个昆虫纲中属于较为原始的类群，约在2.2亿年前，由一个共同祖先进化而来(Songetal.,2015)。大多数种类分布于热带和亚热带，是生物多样性和动物区系的重要组成部分，世界已知 9科 270属1 880多种，我国已知7科57属788种，占世界三分之一多(邓维安,2016)。由于蚱类昆虫既不是益虫也不是害虫，所以国内外对蚱类昆虫的研究不多，主要集中在分类学方面，其它方面较少。但近年来发现，由于蚱类昆虫多以苔藓、地衣或腐殖质为食，属地栖性昆虫，运动能力不强，与环境特征联系紧密，可以作为重要的生态指示物种(邓维安等,2007;Tsuruietal.,2010;Wennerstenetal.,2012)。因此有必要采用新的研究技术，对这一古老的类群进行多方面的研究，以挖掘该类群的科学和应用方面的研究价值。

眼优角蚱Eucriotettixoculatus属于刺翼蚱科Scelimenidae优角蚱属Eucriotettix，广泛分布于中国南部各省及南亚地区，属于蚱类中的广布和常见种类(郑哲民,2005)，具有代表性。本研究基于转录组高通量测序技术，对眼优角蚱E.oculatus雌雄两性虫体之间的转录差异进行研究，以从转录组层面了解蚱总科昆虫雌雄分化的分子基础，为昆虫分子生物学、发育生物学和进化生物学提供一定的基础资料。

1 材料与方法

1.1 材料准备

为了尽可能避免不同环境因素和发育时期引起的雌、雄虫之间转录组的差异，本研究所用的眼优角蚱是2019年6-7月分别在广西河池市宜州区和南丹县不同地点采集的雌性和雄性成虫带回实验室交配产卵后，取其卵块在昆虫饲养室中孵化并在稳定的环境条件下从1龄饲养到成虫后收集的进入成虫期时间相近的雌、雄成虫样品。雌虫样品组(编号AF)由9个生物学重复样品(编号为AF1-AF9)组成，且每个生物学重复样品均由3头雌虫组成。雄虫样品组(编号AM)也是由9个生物学重复样品(编号为AM1-AM9)组成，每个样品也由3头雄虫组成。所有样品均为活虫经48 h饥饿后用液氮冻存。

1.2 转录组测序

使用QIAGEN动物组织RNA提取试剂盒提取各样品总RNA并构建cDNA文库，然后基于Illumina技术测序平台，采用双末端测序(Pair-End)法，完成18个样品二代转录组测序分析。

1.3 基因差异表达分析

本研究测定的所有二代转录组数据与本实验室前期已测定的眼优角蚱三代转录组数据(由卵和各龄期雌、雄虫所组成的混合样品经Pacbio Sequel II技术测序平台利用单分子实时测序法获得的三代数据，编号为EO)联合分析。首先使用二代转录组数据对三代转录组数据进行转录本的校正，然后利用CD-HIT软件(Fuetal.,2012)对校正后的consensus序列进行去冗余，接着对去冗余的序列基于Nr(Lietal.,2002)，Nt，Pfam(Finnetal.,2008)，KOG/COG(Tatusovetal.,2003)，Swiss-prot(Amosetal.,2000)，KEGG(Kanehisaetal.,2004)，GO(Ashburneretal.,2000)共7种数据库进行基因功能注释。同时以去冗余后得到的转录本作为基因的参考序列(ref)，将Illumina测序得到的二代转录组测序数据中每个样品的clean reads比对到ref上。利用RSEM软件(Lietal.,2011)对比对结果进行统计，进一步得到每个样品比对到每个基因上的readcount值，并对其进行FPKM转换。基于基因表达的RPKM值，使用DEseq差异表达基因鉴定程序计算AF和AM样品组之间表达量不同的基因。差异基因筛选标准为|log2 (FoldChange)|>0和qvalue<0.05。

1.4 差异基因表达的GO和KEGG分析

采用GOseq(Youngetal.,2010)和KOBAS软件(Chenetal.,2011)对差异基因集进行GO功能富集分析，KEGG通路富集分析。富集分析基于超几何分布原理，其中差异基因集为差异显著分析所得差异基因并注释到GO或KEGG数据库的基因集，背景基因集为所有进行差异显著分析的基因并注释到GO或KEGG数据库的基因集。富集分析结果是对差异比较组合的所有差异基因集、上调基因集、下调基因集进行富集。

1.5 实时荧光定量PCR验证

从眼优角蚱雌性和雄性转录组中筛选出的差异表达基因中随机挑取9个，以β-actin作为内参基因，采用实时荧光定量PCR (qRT-PCR)的方法对转录组的准确性进行验证(验证样品与转录组测序样品一致)，验证基因和对应的引物序列见表1。使用反转录试剂盒将总RNA反转录合成cDNA，反转录体系为40 μL(冰上配制)：5×PrimeScriptRTMaster Mix (Perfect Real-time) 8 μL和RNA模板2 000 ng，最后加入RNase-free dH2O补全至40 μL。反应条件为：37°C 15 min，85°C 5 s，置于4℃冰箱保存。采用ChamQTM SYBR®qPCR Master Mix染料试剂盒进行实时荧光定量PCR扩增，PCR反应体系20 μL：SYBR Premix Ex Taq (2×) 10 μL，10 μmol/L上、下游引物各1 μL，cDNA模板2 μL，ddH2O补全至20 μL。PCR反应程序为：95℃预变性30 s；95℃变性5 s，60℃退火34 s，40个循环。

表1 实时荧光定量引物Table 1 Primer sequences for Real-time quantitative PCR

2 结果与分析

2.1 转录本校正和去冗余

通过本研究测定的二代转录组数据对三代转录组数据校正后，三代consensus的碱基数目有所提升，其它指标基本无变化(表2)，说明眼优角蚱转录组二代和三代数据整体质量较高。进一步按照序列之间95%的相似性对校正后的转录本序列进行聚类去冗余，最终得到13 432条transcripts，总长度26 977 000 bp，平均长度2 009 bp，最大长度8 297 bp，N50长度远大于1 000 bp(表3)，说明测序效果较好，可以进行后续分析。

表2 转录本校正前后长度分布统计表Table 2 Statistical table of length distribution of transcripts before and after correction

表3 去冗余后序列长度分布统计表Table 3 Statistical table of sequence length distribution after de-redundancy

2.2 基因功能注释

对去冗余之后的序列进行基因功能注释，通过七大数据库共有12 074条transcripts被注释，占整个transcripts的90%，未成功注释的transcripts仅占比10%，说明注释成功率比较高(图1)。对Nr数据库所获得的基因注释与已公布的基因组信息的物种进行匹配，其中被匹配最多的物种有20个(图2)，分别是内华达古白蚁Zootermopsisnevadensis、亚利桑那沫蝉Clastopteraarizonana、飞蝗Locustamigratoria、麦茎蜂Cephuscinctus、赤拟谷盗Triboliumcastaneum等皆为昆虫纲物种，因此可证明本研究转录组注释结果的可靠性。

图1 七大数据库注释结果统计图Fig.1 Statistical chart of annotation results of seven databases

图2 Nr数据库物种注释统计图Fig.2 Statistical chart of species annotation in NR database

2.3 差异表达基因筛选

为了获得雌性和雄性眼优角蚱的表达量不同的基因，运用FPKM法对雌性组(AF)和雄性组(AM)的transcripts进行分析，共获得两者之间的差异表达基因3 597个，其中雌性组(AF)相对于雄性组(AM)，共获得1 295个显著上调基因，获得2 302个显著下调基因(图3)。

图3 眼优角蚱雌性与雄性间差异表达基因的比较Fig.3 Comparison of differential expressed genes between male and female Eucriotettix oculatus

2.4 差异表达基因GO富集结果

对差异表达基因进行GO功能分析可以了解差异表达基因的生物学意义。差异表达基因GO功能分类显示，共有2 344条transcripts被GO注释，并被显著富集到生物学过程(Biological process)、细胞组分(Cellular component)和分子功能(Molecular function)三大类的主要35个亚类中。其中细胞成分组织或合成(Cellular component organization or biogenesis)、细胞成分组织(Cellular component organization)、细胞器组织(Organelle organization)、细胞(Cell)、细胞成分(Cell part)、细胞内(Intracellular)、细胞内成分(Intracellular part)、有机环状化合物结合(Organic cyclic compound binding)、杂环化合物结合(Heterocyclic compound binding)、核酸结合(Nucleic acid binding)分别是三大类中显著富集差异基因数量最多的条目(图4)。

图4 差异表达基因GO功能分类Fig.4 Go functional classification of differentially expressed genes

雌性相对于雄性，表达上调的1 295个基因中，被GO注释的有922个，但没有显著富集的亚类。表达下调的2 302个基因中，被GO注释的有1 422个，显著富集到三大类的47个亚类。其中富集差异基因数目最多的是分子功能(Molecular function)大类中的蛋白质结合(Protein binding)亚类的468个。其次是细胞组分(Cellular component)大类中的细胞器组分(Organelle part)和细胞内细胞器组成(Intracellular organelle part)，分别是157和155个。再次是生物过程(Biological process)中的细胞成分组织或生物发生(Cellular component organization or biogenesis)和细胞成分组织(Cellular component organization)，分别是140个和133个(图5)。其余亚类富集到的基因数据均在100以下。

图5 下调表达基因GO功能分类Fig.5 Go functional classification of down regulated expressed genes

2.5 差异表达基因KEGG富集结果

雌性相对于雄性，表达上调的1 295个基因经过KEGG注释，共有1 188个基因被注释到248条通路，富集显著的前20条通路(图6)，其中神经活性配体受体相互作用(Neuroactive ligand-receptor interaction)、DNA复制(DNA replication)、蛋白质消化吸收(Protein digestion and absorption)、错配修复(Mismatch repair)、嘧啶代谢(Pyrimidine metabolism)和胰腺分泌(Pancreatic secretion)通路是富集最显著的6条通路(qvalue<0.05)。表达下调的2 303个基因中共有1 886个基因被注释到264条通路，富集显著的前20条通路(图7)，其中蛋白酶体(Proteasome)、氮代谢(Nitrogen metabolism)、泛素介导的蛋白水解(Ubiquitin mediated proteolysis)、糖酵解/糖异生(Glycolysis/Gluconeogenesis)、丙酮酸代谢(Pyruvate metabolism)和RNA转运(RNA transport)通路是富集最显著的6条通路(qvalue<0.05)。

图6 雌性相对于雄性上调基因KEGG富集散点图Fig.6 Scatter plot of KEGG enrichment in females relative to males of up-regulated gene

图7 雌性相对于雄性下调基因KEGG富集散点图Fig.7 Scatter plot of KEGG enrichment in females relative to males of down-regulated gene

2.6 实时荧光定量PCR验证结果

qRT-PCR验证结果显示，上调表达的5个基因(VgR、VG2、PRSS1_2_3、TUBA、PLIN2)在转录组测序中的表达均是呈上调趋势，下调表达的4个基因(gdhA、PKLR、UBE2O、CA)在转录组测序中的表达均是呈下调趋势，其中gdhA和PKLR在qRT-PCR验证中均是下调1倍(图8)。

图8 眼优角蚱雌性与雄性差异表达基因qRT-PCR与RNA-Seq分析结果Fig.8 Results of qRT-PCR and RNA-Seq analysis of differential expressed genes between male and female Eucriotettix oculatus

3 结论与讨论

动物的雌雄分化是进化的产物。昆虫作为生物多样性最丰富的动物类群，大多数种类在形态上都表现出了性别二态性，即同种昆虫的雌、雄个体除生殖器官的结构差异外，在大小、颜色、结构等方面也常有明显差异(王孟卿等,2005)。但是，目前对昆虫雌雄性别在转录组层面的差异了解还不充分。

眼优角蚱雌、雄成虫在外部形态上也具有明显的差异，雌虫腹部末端具有细长的产卵瓣，体型较雄虫粗大，前胸背板背面较雄性粗糙，其它构造与体色与雄虫相似(邓维安等,2007)。作为昆虫中较古老的物种，对眼优角蚱雌性和雄性成虫的转录差异开展研究，可为揭示昆虫雌雄差异的分子机制提供一定的基础数据。本研究测定了2组18个样品的二代转录组，并与本实验室已测定的三代转录本进行校正后，最终得到13 432条转录本，平均长度2 009 bp，最大长度8 297 bp，N50长度远大于1 000 bp，说明本研究所得到的转录本质量较较高。本研究仅有10%的转录本没有注释成功，主要是因为数据库中有关直翅目昆虫的基因信息匮乏。

本研究对眼优角蚱雌性和雄性的转录组进行了基因差异表达分析，共获得1 295个显著上调和2 302个显著下调基因。使用qRT-PCR方法验证了9个基因的表达情况，虽然qRT-PCR结果与mRNA测序结果在倍数之间有稍许差异，但是每个基因有同样的上下调趋势，可以说明转录组数据分析的可靠性。在上调的基因中有与雌性发育有关的卵黄蛋白原受体VgR、卵黄蛋白原2VG2、胰蛋白酶PRSS1_2_3、围脂滴蛋白2PLIN2等基因，在下调的基因中有与雄性发育有关的碳酸酐酶CA等基因，说明了差异表达基因与性别紧密联系。为了解差异基因的生物学意义，通过差异基因的GO富集可以发现，差异基因显著富集的GO条目大多是与细胞的形成有关，特别是雄性，其相对于雌性来说，结合(Binding)、细胞组分(Cellular component)和细胞过程(Cellular processes)等条目相关基因显著高表达，说明雌、雄个体在成虫后都有大量的细胞增殖活动，雄性较雌性更加剧烈，这应该是与雄性进入成虫期后要形成更多的精子为交配繁殖做好准备有关(刘杰等,2020)。同时，基于整个基因组背景，应用超几何检验对差异基因进行了的KEGG Pathway富集分析，发现在雌性成虫体内高表达的基因显著富集在主要的6条通路中。其中神经活性配体受体相互作用(Neuroactive ligand-receptor interaction)通路属于生殖相关通路，该通路的基因能影响促性腺素合成和释放，在性腺发育过程中主要是分别作用于脑和卵巢(储明星等,2009;Xuetal.,2015;Zhangetal.,2019)。蛋白质消化吸收(Protein digestion and absorption)和胰腺分泌(Pancreatic secretion)通路属于消化系统，说明雌性眼优角蚱相对雄性消化能力更强，可能是由于雌性个体较大，需要消化吸收更多的营养物质，或者是雌性要储备更多营养物质以备后期生殖需要。而DNA复制(DNA replication)、错配修复(Mismatch repair)和嘧啶代谢(Pyrimidine metabolism)通路均与遗传物质的复制和修复有关。雄性体内高表达的基因也是显著富集在主要的6条通路。其中蛋白酶体(Proteasome)通路可能与雄性眼优角蚱的生精过程有关，因为已有研究发现敲除蛋白酶体REGγ基因可导致精原干细胞数目减少，使雄性小鼠生殖能力降低(Alietal.,2013;李宁等,2016)。氮源是昆虫生长发育和繁殖所必需的营养物质之一。由于食物中氮含量与蛋白质含量直接相关，食物中的氮能直接影响昆虫的交配、取食速率和繁殖等(张碧尧,2019)，研究发现与氮代谢有关的精氨酸对精子发生发挥关键作用(Moralesetal.,2003;O’Flahertyetal.,2004;徐学玉等,2017)，因此雄性眼优角蚱氮代谢通路被富集。泛素介导的蛋白水解(Ubiquitin mediated proteolysis)属于信号通路，参与卵发生、精子发生和性腺成熟等一系列生殖过程的调控活动(Mengetal.,2015;Pengetal.,2015;Yangetal.,2016;董忠典等,2021)。与能量代谢相关的糖酵解(Glycolysis / Gluconeogenesis)、丙酮酸代谢(Pyruvate metabolism)通路在雄性中被显著富集，可能是雄性眼优角蚱在交配前较雌性活动多，需要消耗更多的能量(肖舒晴等,2019)。RNA转运(RNA transport)通路显著富集，可能是由于蛋白质合成与RNA转运密切相关，雄性需要大量合成与精子形成有关的蛋白，以完成精子储备，这与雄性按蚊转录表达相似(刘杰等,2020)。综上所述，眼优角蚱雌性和雄性成虫在转录水平的差异主要还是集中在与生殖细胞形成和生殖活动有关的基因表达方面。本研究对进一步深入解析昆虫雌雄分化的分子调控差异提供了一定的理论基础。