黑翅土白蚁工蚁和兵蚁转录组测序比较及表达分析

黄双

摘    要：基于Illumina高通量测序RNA-seq技术对黑翅土白蚁Odontotermes formosanus （Shiraki，1909）的工蚁和兵蚁进行了转录组测序分析，了解不同品级的转录组特征，丰富白蚁转录组数据信息，为探明黑翅土白蚁品级分化的分子机制打下基础。经过测序获得质量值高于20的碱基比例（Q20）均不低于98%的数据，所有reads组装成65 604个Unigenes，平均长度337.78 bp;将所得序列注释到NR、NT、KOG、GO和KEGG等数据库进行比对，所得Unigenes均被注释。兵蚁与工蚁差异基因的上调基因有与能量代谢及蛋白质翻译有关，下调相关基因有与能量运输相关，差异基因丰富，有助于黑翅土白蚁品级分化相关基因挖掘。本研究扩充了黑翅土白蚁的基因资源，为科学研究白蚁种群的基因构造、遗传社会性、获取有效分子标记、综合防控等提供珍贵的生物资源。

关键词：黑翅土白蚁;转录组;工蚁;兵蚁;比较分析

中图分类号：Q963         文献标识码：A            DOI 编码：10.3969/j.issn.1006-6500.2022.06.001

Comparative Transcriptome Sequencing and Expression Analysis of Worker and Soldier Termites of Odontotermes Formosanus

HUANG Shuang

（Key Laboratory of Southwest China Wildlife Resources Conservation （ Ministry of Education） ，China West Normal University， Nanchong ， Sichuan 637002， China）

Abstract： Based on Illumina high-throughput sequencing RNA-seq technology， the transcriptome of worker and soldier termites of Odontotermes formosanus （Shiraki， 1909） was sequenced to understand the transcriptome characteristics of different grades， to enrich the information of termite transcriptome data， and to lay the foundation for exploring the molecular mechanism of grade differentiation in the black-winged soil termite. After sequenced， we obtained data that the proportion of bases with quality values higher than 20 （Q20） was not less than 98%， and all reads were assembled into 65 604 Unigenes with an average length of 337.78 bp. The obtained sequences were annotated to NR， NT， KOG， GO and KEGG databases for comparison， and all Unigenes were annotated. The up-regulated genes were related to energy metabolism and protein translation， and the down-regulated genes were related to energy transport， which were rich in differential genes and could help to further develop their work to find their class differentiation-related genes. The study expands the genetic resources of black-winged soil termites and provides valuable biological resources for scientific research on the genetic structure of termite populations， genetic sociality， acquisition of effective molecular markers， and integrated prevention and control.

Key words： Odontotermes formosanus （Shiraki， 1909）; transcriptome; worker ants; solider ants; comparative analysis

黑翅土白蚁Odontotermes formosanus（Shiraki，1909）为白蚁科（Termitidae）大白蚁亚科（Macroter-mitinae）土白蚁属（Odontotermes），屬土栖性白蚁，主要分布在我国北纬35°以南的地区，分布范围广，食性杂，主要以树根、皮层、幼苗、幼树、庄稼秸杆、甘蔗、竹根、草根、水稻等为取食对象，也能取食牛粪[1-2]。黑翅土白蚁品级包括蚁王与蚁后、有翅成虫、兵蚁、工蚁4个品级[3]，各品级白蚁具有不同的行为模式和社会分工，其中蚁后与有翅成虫为繁殖蚁，而工蚁与兵蚁为不育蚁，其形态结构及社会职能为不同属性[4]。根据巢体大小其卵分化为不同品级，其分化时间及龄期会不同。黑翅土白蚁工蚁数量全巢最多，可达到200万个，从卵孵化到工蚁需要17 d;兵蚁数量次之，兵蚁有雌雄之分，但无交配生殖能力，从卵孵化到兵蚁分化大约需要24 d[5]。工蚁社会属性为作工，主要有筑巢、供食、修路、照顾幼蚁、觅食等;兵蚁社会属性为保卫巢穴，用发达的上颚和能分泌毒液的额腺作为防御武器[6]。卵发育为工蚁或兵蚁受环境影响，且在最后1个龄期前其品级完成转化，通过工蚁与兵蚁的高通量转录组测序比较分析，可为影响分化品级基因发生的分子机制提供研究信息[7]。本研究采用Illumina Hiseq 2000高通量测序技术对黑翅土白蚁工蚁与兵蚁进行测序及生物信息学分析，建立转录组数据库，分析兵蚁和工蚁表达差异显著的基因，为其分化品级相关分子机制的研究提供依据[8]。

1 材料和方法

1.1 试验材料

2020年6月在四川南充西山进行了黑翅土白蚁样品的采集，采集后立刻将其保存在RNA保存液中。

1.2 试验方法

1.2.1 总RNA提取 选取样品中的兵蚁20头和工蚁30头，分别取下头部，成功地从工蚁和兵蚁头部提取总RNA。

1.2.2 转录组测序 通过Illumina Hiseq 2000平台测序，构建cDNA文库，在无参考基因组的条件下，通过样本间比较对差异基因进行分析。其中，对照组为工蚁，处理组为兵蚁。测序由上海生工生物有限公司完成。

1.2.3 数据处理 在Hiseq平台上经CASAVA碱基识别 （Base Calling） 分析转化得到原始测序序列（Sequenced Reads），黑翅土白蚁的兵蚁和工蚁的转录组原始数据通过FastQC进行原始数据质量评估，使用Trimmomatic[9]进行质量剪切，得到clean 数据。再使用Trinity软件[10]组装成转录本，得到转录组拼接序列（Transcript）和冗余后的拼接转录本序列（Unigene）。使用Salmon[11]计算基因的表达量，Unigene表达量的计算使用TPM（Transcripts Per Million）。计算方式为：

TPMi=××106

Xi=total exon fragment/reads Li=

1.2.4 Unigene的注释、功能分类和生物学通路分析 通过BLAST[12]将拼接后的Unigene序列与CDD、KOG、NR、NT、PFAM、Swissprot、TrEMBL等多个数据库进行比对，摘取重要数据库比对信息得到其重要功能注释信息。比对到NCBI非冗余核酸数据库（NCBI non-redundant protein sequences，NR）中，分析物种相似性。根据转录本与Swissprot、TrEMBL的注释结果得到基因功能分类数据库（Gene Ontology，GO）功能注释信息，从宏观上认识黑翅土白蚁的工蚁及兵蚁的基因和基因产物的属性。对Unigene进行真核蛋白质直系同源数据库（Clusters of orthologous

groups for eukaryotic complete genomes，KOG）比对，对其可能的功能分类统计[13]。利用KAAS得到转录本基因组百科全书代谢通路库（Kyoto encyclopedia of genes and genomes，KEGG）注释信息，了解其代谢途径、生物学功能以及基因间的相互作用等[14]。

1.2.5 工蚁和兵蚁差异基因表达数据分析处理  利用聚类分析软件对表达模式相似的不同基因进行聚类，并用数据分析比较基因差异。

2 结果与分析

2.1 工蚁和兵蚁RNA-seq转录组数据拼接及注释结果

2.1.1 工蚁和兵蚁RNA-seq转录组数据拼接结果 黑翅土白蚁兵蚁共获得24 138 602条原始序列，Q20含量98.93%，GC含量62.28%。工蚁共获得39 376 474条原始序列，Q20含量98.97% ，GC含量57.40%，具体见表1。转录组测序结果较好，能够满足后续数据分析。

将转录组进行de novo拼接，这些干净读数组装成 99 639 条Transctipt，这些转录本最后被组装成 65 604条Unigenes。对已组装的转录本和Unigenes进行了长度分布统计，具体见表2。

2.1.2 工蚁和兵蚁Unigene序列的比对及NR注释  用NCBI Blast+将工蚁和兵蚁所有Unigenes序列比对到PFAM、Swissprot、TrEMBL、CDD、KOG、NR、NT7个数据库，其注释基因分别有8 738，18 326，

23 458，13 235，11 562，23 646，13 256条。NR数据库可查阅同源序列的功能信息及与相似物种的近似值，注释匹配较高的相似物种分别为内华达州湿木白蚁（Zootermopsis nevadensis）、小鼠（Mus muscuius）、家白蚁（Coptotermes formosanus），分别有8 549，

2 477，789条相似序列。

2.1.3 工蚁和兵蚁Unigene的GO注释及其分类  根据NR数据库注释所得，对21 285条Unigenes（占总Unigenes的32.44%）进行了GO功能注释，共注释206 751个GO功能，可了解其基因功能分布特征。在GO注释中有3大类别，分别是生物过程、细胞组分和分子功能。如图1所示，生物过程共注释了92 356个基因，数量最多，占总注释的44.67%，其包含的26个功能亚类中，以细胞进程（Cellular Processes）和代谢过程（Metabolism）、生物调节（biological regulation）占主导地位，分别占该类型15.95%，12.76%，8.87%，以生物相（biological phase）最低，仅占生物过程全部的0.000 3%。其次是細胞组分，共注释了82 799个基因，占40.04%，其包含的22个功能亚类中，细胞（cell）、细胞组分（cell part）、细胞器（organelle）占主导地位，分别占该类别的18.88%，18.81%，13.97%，而以其他细胞器（other organism）和其他细胞器组成（other organism part）较低，均占该类别的0.02%。分子功能最少，共31 596个基因被注释，占15.28%，其包含的21个功能亚类中，以结合（binding）、催化活性（catalytic activity）居多，分别占该类型的42.41%，32.18%，而以形态发生素活性（morphogen activity）、化学排斥物活性（chemorepellent activity）、营养库活性（nutrient reservoir activity）、金属伴侣蛋白活性（metallochaperone activity）较低，分别占该类别的0.001%，0.019%，0.022%，0.025%。这些结果显示了黑翅土白蚁头部基因表达情况，以生物活性及代谢调节相关的基因量较多，表明黑翅土白蚁代谢能力较强。

2.1.4 黑翅土白蚁转录组Unigenes的KOG注释及其分类 在KOG注释成功的基因有11 562条Unigenes，占总Unigene的17.62%。按其数据库功能分类为25个。其比对结果如图2所示，信号传导机制（Signal transduction mechanisms）最多，占12.8%，其次是翻译后修饰、蛋白质折叠和分子伴侣（Posttranslational modification，protein turnover，chaper-

ones，11.85%），一般功能基因（General function prediction only，9.89%），翻译、核糖体结构与生物起源（Transiation，ribosomal structure and biogenesis，7.95%）;而防御机制（Defense machanisms， 0.72%）、核结构（Nuclear structure，0.36%）和细胞运动性（Cell motility，0.10%）是最小的类群。结果表明，黑翅土白蚁在信号传递和转录、翻译等基因表达丰度高。

2.1.5 工蚁和兵蚁差异表达基因Unigene的KEGG代谢通路分析 根据KO与Pathway的关联性对KEGG代谢通路进行分析。结果如图3所示，共有9 539 Unigenes获得注释（14.54%），共获得6 657个KEGG功能注释。黑翅土白蚁Unigenes参与或涉及的代谢通路可归为5个类别，即环境信息处理（Environmental Information Processing）、代谢（Metabolism）、细胞过程（Cellular Processes）、生物系统（Organismal Systems）、遗传信息处理（Genetic Information Processing），共有33个亚类293个信号通路（Pathway）。大多数Unigenes都被注释到了代谢、有机系统。参与代谢的共有2 271条Unigenes，主要集中在参与糖代谢过程（Carbohydrate metabolism）的Unigenes数量为419条，比例为18.45%;其次参与能量代谢过程（Energy metabolism）的为328条（14.44%）;参与全局整体映射（Overview）和氨基酸代谢（Amino acid metabolism）的Unigenes分别有301条（13.25%）和279条（12.29%）。参与生物系统的共有1 415条Unigenes，主要集中在内分泌系统（Endocrine system，21.98%）和免疫系统（Immune system，18.80%）。参与细胞过程的共有 895条Unigenes，主要集中在运输与代谢（Transport and catabolism，38.77%）和細胞群落（Cellular community，28.38%）。环境信息处理涉及912条Unigenes，包括信号转导（Signal transduction，79.50%），信号分子和交互作用（Signaling molecules and interaction，13.38%），膜转运（Membrane transport，7.13%）。基因信息处理涉及1 164条Unigenes，主要涉及翻译（Translation，46.48%），折叠、分类和降解（Folding sorting and degradation，31.27%）。

2.2 工蚁和兵蚁差异基因表达的分析

2.2.1 工蚁和兵蚁差异基因Unigenes的GO富集分析 所有Unigenes比对到数据库中所有差异基因数为53 092个，其中工蚁和兵蚁差异表达基因GO生物过程比对基因数为16 957个，如图4所示，最显著富集基因数为3 472个，主要用于发育过程（developmental process）、多细胞生物过程（multicellular organismal process）、生物调节（biological regulation）、生物过程调节（regulation of biological process）、行为（behavior）、信号（signaling）等;细胞组成中比对基因数有3 908个，其中最显著富集基因数为3 908个，主要用于膜部分（membrane part）、突触（synapse）、突触部分（synapse part）;分子功能比对基因数为18 062个，其中最显著富集基因数为3 810个，主要功能是结构分子活性（structural molecule activity）、受体活性（receptor activity）、运输活性（transporter activity）等。

2.2.2 工蚁和兵蚁差异基因Unigenes的KEGG富集分析 KEGG通路中所有基因数3 591个，其中差异表达基因数890个，差异表达基因显著富集通路20个（表3），其中5个是与营养因子合成及代谢相关的通路，5个是与蛋白质、激素合成相关的通路，3个是与免疫应答相关的通路，7个是与信号传导相关通路。

2.2.3 工蚁和兵蚁差异基因表达的分析 在兵蚁和工蚁的处理组中，有20 486个差异基因，筛选出10 503个具有显著表达差异基因，其中上调基因有5 099个，下调基因有5 404个（图5）。具有功能注释的上调基因有2 291个，具有功能注释的下调基因有4 063个。工蚁和兵蚁一共有27 167个相同的差异表达基因，12 279条有功能注释。工蚁有19 889个特有差异基因，7 969条有功能注释;兵蚁有17 868个特有差异基因，9 670条有功能注释。

2.2.4 工蚁和兵蚁相比差异幅度最大且前18位的基因 与兵蚁头部基因表达相比，工蚁头部中上调或下调幅度前18位基因见表4。统计学差异显著性检验指标P<0.01;q为校正后的P值，q<0.01。从表4可以看出，兵蚁上调幅度最大前7位基因中有主要功能为碳水化合物运输和代谢及翻译、核糖体结构和生物生成，其余功能未知;在工蚁中下调幅度最大的前11位基因主要功能为能量运输相关，其余功能未知。

3 结论与讨论

本研究中，采用Illumina hisiqTM高通量测序对黑翅土白蚁的兵蚁和工蚁的转录组差异进行了分析，获得黑翅土白蚁转录组信息。测序结果发现，共获得0.58 Gb的clean reads，每个样品的Q20含量均在98%以上，说明测序质量较好。经过组装拼接构建Unigenes数据库，共获得65 604个Unigenes，平均长度为337.78 bp，N50为319 bp，N90为216 bp，且多数Unigenes能比对到亲缘关系较近的物种，说明其测序组装质量较好。采用生物信息学分析方法，对黑翅土白蚁兵蚁及工蚁转录组的序列进行功能注释及其功能分类。通过NR数据库的比对，有23 646条序列比对上，占总Unigene的36.04%。在KOG分类中，共得到11 562个功能注释，涉及25个KOG类别，其中信号传导机制比例最大（1 632个，12.80%），其次是翻译后修饰、蛋白质转换、分子伴侣（1 511个，11.85%），防御机制比例较小（92个，0.72%），最小是核结构（46个，0.36%），这些基因反映了兵蚁与工蚁在一定时期的表达情况。在NR数据库中发现剩余的未被注释的基因在无参考基因组的条件下，可能是黑翅土白蚁自身特有的基因，或是含有蛋白质编码未翻译或非保守的一段区域[15]。

在GO分类中，获得206 751个功能性注释，占生物过程注释总数的44.67%，其次是细胞成分（40.04%）和分子功能（15.28%）。在KEGG分类中，有3 591条Unigenes具有具体的定义，共获得6 657个KEGG功能注释，在5大分支中，基因主要参与信号传导、翻译、糖代谢、折叠排序和降解。对黑翅土白蚁兵蚁及工蚁差异表达基因进行分析，发现细胞胚胎发育、生物调节、膜部分、突触等基因表达均存在差异。黑翅土白蚁中兵蚁及工蚁的差异基因表达分析结果可根据qPCR试验进行进一步分析。兵蚁失去自身取食能力，需由工蚁将食物吞入自身体内然后将已经消化或半消化的营养从口吐出或肠管下侧排出予以喂饲才能存活。因此，转录组分析可进一步探索和发现纤维素利用、消化酶的调控基因[16]。

通过转录组测序发现有65 604条Unigenes具有注释信息，占总Unigenes的100%。其中，GO分析得到69个功能分类，以生物学过程、细胞组分占大多数，各分亚类中以细胞组分、细胞、细胞进程等类别的基因含量丰富，表明黑翅土白蚁具有很强的细胞活性机制;KOG分析可分为25个不同功能分类，以信号传导机制和转录功能、一般功能基因的较多;KEGG的293个代谢通路，以代谢、有机系统较为丰富。兵蚁与工蚁差异基因的上调基因有与能量代谢及蛋白质翻译有关，下调相关基因有与能量运输相关，差异基因丰富，有助于黑翅土白蚁品级分化相关基因发掘。黑翅土白蚁基因丰富，分析其基因分布特征，確定基因代谢通路的分布，为黑翅土白蚁基因提供了转录组水平上的支持，也为后续黑翅土白蚁基因构造，遗传社会性，获取有效分子标记，综合防控等提供珍贵的生物资源。

参考文献：

[1] 黄复生， 朱世模， 平正明， 等. 中国动物志 昆虫纲 第十七卷 等翅目[M]. 北京： 科学出版社， 2000.

[2] 蔡邦华， 陈宁生. 中国白蚁分类和区系问题[J]. 昆虫学报， 1964， 13（1）： 25-37.

[3] 徐志德， 李德运， 周贵清， 等. 黑翅土白蚁的生物学特性及综合防治技术[J]. 昆虫知识， 2007， 44（5）： 763-769.

[4] 陈艳. 黑翅土白蚁巢群间竞争和种间竞争的研究[D]. 武汉： 华中农业大学， 2007.

[5] 黄求应. 黑翅土白蚁觅食行为学基础及诱杀系统的研究[D]. 武汉： 华中农业大学， 2006.

[6] PRESTWICH G D. Chemical systematics of termite exocrine secretions[J]. Annual Review of Ecology and Systematics， 1983， 14（1）： 287-311.

[7] MIURA T， KAMIKOUCHI A， SAWATA M， et al. Soldier caste-specific gene expression in the mandibular glands of Hodotermopsis japonica （Isoptera： Termopsidae）[J]. Proceedings of the National Academy of Sciences of the United States of America， 1999， 96（24）： 13874-13879.

[8] LEONARDO F C， DA CUNHA A F， DA SILVA M J， et al. Analysis of the workers head transcriptome of the Asian subterranean termite， Coptotermes gestroi[J]. Bulletin of Entomological Research， 2011， 101（4）： 383-391.

[9] BOLGER A M， LOHSE M， USADEL B. Trimmomatic： a flexible trimmer for Illumina sequence data[J]. Bioinformatics， 2014， 30（15）： 2114-2120.

[10] HAAS B J， PAPANICOLAOU A， YASSOUR M， et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis[J]. Nature Protocols， 2013， 8（8）： 1494-1512.

[11] PATRO R， DUGGAL G， LOVE M I， et al. Salmon provides fast and bias-aware quantification of transcript expression[J]. Nature Methods， 2017， 14（4）： 417-419.

[12] ALTSCHUL S F， MADDEN T L， SCHÄFFER AA， et al. Gapped BLAST and PSI-BLAST： a new generation of protein database search programs[J]. Nucleic Acids Research， 1997， 25（17）： 3389-3402.

[13] TATUSOV R L， GALPERIN M Y， NATALE D A， et al. The COG database： a tool for genome-scale analysis of protein functions and evolution[J]. Nucleic Acids Research， 2000， 28（1）： 33-36.

[14] KANEHISA M， GOTO S. KEGG： kyoto encyclopedia of genes and genomes[J]. Nucleic Acids Research， 2000， 28（1）： 27-30.

[15] HOWARD R， HAVERTY M I. Seasonal variation in caste proportions of field colonies of Reticulitermes flavipes （Kollar）[J]. Environmental Entomology， 1981， 10（4）： 546-549.

[16] MCMAHAN E A. Studies of termite wood-feeding preferences[J]. Proceedings， Hawaiian Entomological Society， 1966， 19（2）： 239-250.