思茅松毛虫不同发育阶段的转录组差异分析

王宇翔，杨斌，吴培福，冷堂健，邱雨，周杰珑，*

（1. 西南林业大学生命科学学院，昆明650224；2. 西南林业大学生物多样性保护学院，云南省森林灾害预警与控制重点实验室，昆明650224）

思茅松毛虫Dendrolimus kikuchii（Matsumura，1927）是我国危害最严重的6 种松毛虫之一，隶属于鳞翅目Lepidoptera 枯叶蛾科Lasiocampidae 松毛虫属Dendrolimus，模式标本来自云南思茅（侯陶谦，1987），在我国分布于云南、四川、贵州、广东、湖南、浙江、江苏、安徽、福建、台湾等省，主要危害思茅松Pinus kesiyavar.langbianensis、云南松Pinus yunnanensis、滇 油 杉Keteleeria fortunei、马 尾 松Pinus massoniana及落叶松Larix gmelinii等针叶类植物（陈昌洁，1990；侯陶谦，1993；曹先聪，2017）。爆发成灾时，针叶短时间被蚕食精光，被害树林远看枯黄一片，如同火烧，俗称“无烟森林火灾”；虫害轻则影响松树生长、落叶及枝条枯死，降低松果、木材和松脂的产量（侯陶谦，1987；Daiet al.，2012；Menet al.，2017），重则造成树木死亡，给林业经济和生态景观造成巨大损失（童清，刘宏屏，2009；聂中良等，2019）。松毛虫与人接触还会引起“松毛虫病”（陈昌洁，1990）。

关于思茅松毛虫的研究集中在危害调查及防治措施、肠道微生物、种群动态监测、线粒体系统发育、组织病理及生理生化等方面（张荣超，2015；曹先聪，2017；Liuet al.，2017；Menet al.，2017；Wuet al.，2017；冯玉元，2018；万鹰等，2018；李雕益等，2022；邱雨等，2022）。Zhou 等（2021）报道了思茅松毛虫的全基因组，为该物种的系统进化和防治策略提供了新的视角。但迄今为止有关思茅松毛虫基因表达及调控相关研究较为匮乏。随着下一代测序技术（NGS）和生物信息学的发展，RNA-Seq分析成为发育生物学、进化生物学、分类学、寄主互作及毒理学等领域的重要研究手段（Zhanget al.，2014；Tianet al.，2015；姚顺，2018；Chenet al.，2019；Norigaet al.，2019），可从整体水平上分析基因结构和功能，有利于发现生物学过程的基因表达变化规律（Wang & Kirkness，2005）。昆虫转录组分析为深入系统地研究昆虫生命活动相关基因功能、系统进化与分类、生长发育以及昆虫与相关寄主互作等提供依据（张棋麟，袁明龙，2013）。基于此，本研究采用DNBSEQ-T7 测序平台对思茅松毛虫不同发育阶段的转录组进行测序，开展各阶段基因表达谱差异和功能富集分析。研究结果能丰富思茅松毛虫遗传资源的组学数据库，有助于揭示基因表达与特定阶段生理功能的关系，理解昆虫全变态发育的生物学过程，为进一步阐明思茅松毛虫生长发育的分子机制、开展功能基因表达调控及害虫生物防治策略相关研究奠定理论基础。

1 材料与方法

1.1 样品采集

思茅松毛虫虫茧采自云南省安宁市县街野外滇油杉林（102°24'E，24°46'N），采回后于室内羽化、交配、产卵、孵化及饲养（饲养条件：28 ℃，相对湿度75%，光照L∶D=16∶8，幼虫饲喂新鲜滇油杉）。收集各发育阶段样品：卵50粒、1龄幼虫10头、4龄幼虫 10 头、7 龄幼虫 6 头、蛹 4 头和成虫 4 头。将上述样品用PBS 溶液（pH7.4，浓度0.01 mol·L−1）清洗干净，−80 ℃保存。

1.2 RNA提取、文库构建及测序

采用以Trizol 法提取总RNA 的试剂盒（TIANGEN，DP431）分别提取不同发育阶段样品的总RNA。质检后，使用MGIEasy RNA 文库制备试剂套装（深圳华大智造科技股份有限公司）建立6个cDNA文库。通过DNBSEQ-T7（原MGISEQ-T7）测序平台上对合格文库进行测序。文库构建及测序由武汉希望组生物科技有限公司完成。转录组测序数据提交至NCBI_SRA 数据库，登录号为PRJNA732117。

1.3 数据质控及比对参考基因组

使用 Fastp（https：//github.com/OpenGene/fastp）对原始数据进行过滤，然后利用FastQC（http：//www.bioinformatics.babraham.ac.uk/projects/fastqc）对clean data 进行质控，质控合格后进行后续分析。使用 Hisat2（Kimet al.，2015）将 clean reads 与参考基因组比对（DDBJ/ENA/GenBank登录号：JAHHIN‑000000000）。

1.4 差异表达基因（DEGs）及其功能富集分析

通过定位到基因组区域或基因外显子区的测序序列（reads）的计数来估计基因的表达水平。本研究以每百万reads 中来自某一基因每千碱基长度的fragments 数目（FPKM）作为每个基因表达量的衡量标准，以FPKM=0.1 作为判断基因是否表达的阈值：0.115 为高丰度表达水平。参考Bardou 等（2014）的方法绘制基因维恩图；各发育阶段基因表达量经log2转换后，通过R 语言绘制聚类热图。采用 EdgeR（Robinsonet al.，2010）对不同发育阶段的DEGs 进行分析，筛选条件为|Fold Change |≥2 且 P-adjust（FDR）≤0.005。 采 用clusterProfiler（https：//github. com/GuangchuangYu/clusterProfiler；Yuet al.，2012）对 DEGs 进行 GO 和KEGG 富集分析，以FDR<0.05 作为筛选显著富集结果的阈值，确定DEGs 主要参与的代谢通路或信号途径。

2 结果与分析

2.1 转录组测序数据统计与质量评估

测序获总的数据量为369 804 176 条序列，各发育阶段样品的为47 517 160~75 263 362条；经数据过滤质控，约有99%被作为clean reads 保留，获得368 027 762 条序列进行下游分析，各阶段clean reads 为 47 433 276~74 726 282 条，其中，成虫期最多，4龄幼虫期最少；各阶段clean reads的比对率为86.29%~91.10%；Q20 均 在 96.46% 以 上 ，Q30 在90.00% 左右；GC 含量为 41.10%~47.60%（表 1）。转录组测序质量较高，数据可满足生物信息学分析相关要求。

表1 不同发育阶段转录组测序数据统计与质量分析Table 1 Statistics and quality analysis of transcriptome data at different developmental stages

2.2 基因表达水平分析

各发育阶段基因数量分别为11 580 个（卵）、12 128 个（1 龄幼虫）、11 334个（4龄幼虫）、11 209个（7 龄幼虫）、12 160 个（蛹）和11 478 个（成虫）。基因数量和高表达基因均为蛹阶段最多，4 龄幼虫期最少；各发育阶段中，除幼虫期外，其他3个阶段均为高表达基因占比最高（表2）。

表2 不同表达水平的基因数量统计表Table 2 Statistics of the number of genes in different expression levels

从维恩图可以看出（图1：A），不同发育阶段均表达的基因有9 251 个；特有基因在卵、幼虫、蛹、成虫阶段分别有 213 个、261 个、183 个、141 个，其中幼虫阶段特有基因分别为153个（1龄幼虫）、43个（4 龄幼虫）、65 个（7 龄幼虫）。热图显示（图1：B），不同发育阶段呈现不同的基因表达聚类模式，其中4 龄幼虫和7 龄幼虫聚为一支，再与1 龄幼虫汇聚，蛹和成虫聚为一支，再与卵汇聚，最后两大支汇合；不同发育阶段的基因表达谱存在明显差异，各发育阶段均有几个主要的高表达基因簇，卵和蛹阶段尤为突出。

图1 基因表达维恩图（A）和基因表达模式聚类图（B）Fig. 1 Venn diagram of gene expression（A）and cluster diagram of gene expression patterns（B）

2.3 DEGs分析

DEGs 总数为 12 927 个，其中，卵与 1 龄幼虫的最多（3 336 个：上调1 902 个，下调 1 434 个），1 龄幼虫与 4 龄幼虫的有 2 398 个（上调 1 362 个，下调1 036 个），4 龄幼虫与7 龄幼虫的有1 758 个（上调642 个，下调 1 116 个），7 龄幼虫与蛹的有 2 034 个（上调1 166 个，下调868 个），蛹与成虫的最少（1 654个：上调821个，下调833个）。

不同发育阶段基因表达的比较结果显示，从卵至1龄幼虫，显著上调基因有肌球蛋白相关基因（XP_026493160.1 和AAV91411.1）、化学感受蛋白相关基因（AII01039.1 和AII01031.1）、视黄醇脱氢酶相关基因（XP_026737674.1）、丝氨酸蛋白酶相关基因（ACR15967.2）；而与有丝分裂特异性细胞周期蛋白cyclin-B1 相关基因（XP_026731829.1）、细胞分裂周期蛋白20（CDC20）相关基因（XP_026733697.1）明显下调。4 龄幼虫相对1 龄幼虫，载脂蛋白相关基因（XP_026760716.1）、细胞壁蛋白相关基因（XP_004932399.1）、表皮蛋白相关基因（XP_026734677.1）、保幼激素相关基因（XP_004925704.1 和 XP_026743329.1）、蜕皮激素相关基因（XP_026761583.1 和XP_026735381.1）等显著上调；7 龄幼虫相对4 龄幼虫，解毒酶细胞色素P450 相关基因（XP_021191965.1 和 AID54859.1）、表皮蛋白相关基因（XP_026762046.1）、保幼激素相关基因（XP_026733514.1和XP_026743329.1）、蜕皮相关基因（XP_026330046.1和XP_026761583.1）等表达上调，而脂肪酶相关基因（AEA76291.1）表达下调。从7龄幼虫到蛹阶段，显著上调基因有溶菌酶相关基因（XP_026321249.1）、脂肪酸合成酶相关基因（XP_026499944.1）、过氧化物酶相关基因（XP_026725384）、长链酯酰辅酶A 合成酶相关基因（XP_026724857.1）、储存蛋白相关基因（XP_026327069.1）、表皮蛋白相关基因（XP_026501352.1）、卵黄原蛋白基因（BAB16412.1）等，显著下调基因有细胞色素P450 家族相关基因（AID54859.1 和（XP_021191965.1）。从蛹到成虫阶段，角质层蛋白相关基因（XP_026756401.1）、肌钙蛋白 C 相关基因（XP_026757608.1）、脂肪酸合成酶相关基因（XP_026499944.1）、化学感受蛋白相关基因（AII01034.1）、飞行有关基因（XP_026317122.1和XP_026757413.1）等显著上调，而与血淋巴蛋白酶相关基因（AAV91004.1）、溶酶体相关基因（XP_026734095.1）、Toll样受体相关基因（XP_026744547.1）表达下调。

2.4 不同发育阶段DEGs功能富集分析

GO分析表明，各阶段DEGs均能显著富集到生物过程、细胞组分及分子功能，其中涉及分子功能的条目最多、细胞组分的最少（图2）。KEGG 分析发现，各阶段DEGs显著富集到代谢、生物体系统、遗传信息处理、细胞过程、环境信息处理5类途径，涉及101条通路，其中代谢途径最多（57条通路；图3）。

图2 不同发育阶段差异表达基因GO富集分析Fig. 2 GO enrichment analysis of differentially expressed genes at different developmental stages

图3 不同发育阶段差异表达基因KEGG通路富集分析Fig. 3 KEGG pathway enrichment analysis of differentially expressed genes at different developmental stages

从卵到1 龄幼虫，DEGs 主要富集于催化活性和代谢过程。相比于1龄幼虫，卵期上调基因主要参与细胞分裂过程和蛋白质合成过程。随龄期增加，GO 富集情况与上、下调基因参与KEGG 通路均有所变化，主要体现在免疫功能和解毒代谢途径上。在幼虫变态为蛹阶段，DEGs 显著富集于催化活性、结构组分和结合分子功能，氧化还原过程、跨膜运输和蛋白质水解（生物过程），核糖体、胞外区域（细胞组分）；上调基因富集于DNA 复制、核苷酸切除修复，下调通路主要有核糖体和大量代谢途径。当蛹变为成虫，DEGs富集于跨膜运输、氧化还原过程、催化活性、细胞粘附及蛋白质折叠的细胞增殖过程；上调基因参与氨基酸代谢、脂肪酸代谢和与飞行有关的昼夜节律等途径。

3 讨论

通过对思茅松毛虫发育过程的转录组学研究和差异表达基因的分子解析，可为后续开展发育生物学研究奠定分子基础，为今后害虫综合防治提供新的方向和策略。本文采用DNBSEQ-T7测序平台对思茅松毛虫进行转录组测序分析，发现各发育阶段差异表达基因与每个阶段特定生理活动存在显著相关性。在卵期，参与细胞周期和细胞分裂相关基因有较高表达。幼虫阶段，与消化酶、代谢、肌肉生长、感觉等相关基因明显上调；随龄期增加，保幼激素、解毒、表皮蛋白及抗性等相关基因也显著上调。蛹阶段，与细胞增殖和分化相关基因被激活，脂肪酸合成酶、储存蛋白、表皮蛋白、溶菌酶等相关基因显著上调，这些可能与保证蛹发育和角质层形成提供足够蛋白质和氨基酸（Pick&Burmester，2009；Allenet al.，2018）、由凋亡和自噬介导对幼虫器官的自我消化（Tashiroet al.，1976；Joneset al.，1995）和免疫防御机制（Chung &Ourth，2000；Xieet al.，2013）等有关。成虫阶段，高表达基因主要与组织器官形成、外部形态、繁殖和飞行等过程显著相关。此外，基因表达具有发育时空特异性，卵与1龄幼虫间差异表达基因数目最多，说明从卵向幼虫的转变是最为复杂的过程，这与戚礼兴（2016）、Yang 等（2016）及 Han 等（2019）在菜粉蝶Pieris rapae、马尾松毛虫Dendrolimus punctatus和云南松毛虫Dendrolimus houi上的研究结果一致。

在GO 功能富集中，各发育阶段的DEGs 显著富集到了氧化还原、蛋白质水解、跨膜运输、碳水化合物代谢等生物过程，这些过程为机体生长发育、运动提供营养和能量来源（Claudiuset al.，2014；Massoet al.，2016；Gautamet al.，2020）。在分子功能中，较多地富集在活性功能和结合功能上。其中，催化活性最为典型，生物体各类催化活动都离不开酶的参与，保证了昆虫物质代谢和营养能量转换等过程正常进行（赵星，2021）。结合功能与一些分子和关键因子/离子需要结合后才能变为活性、有功能的调控有关。结合功能和催化活性功能显著富集情况与姚顺（2018）的研究结果较为一致。几丁质结合是唯一一个在各发育阶段比较中均被显著富集的条目，这可能与其全变态发育的周期性表皮降解和重构等生理过程有关（Co‑hen，2001）。在细胞组分中，显著富集到胞外区域、线粒体内膜、质膜等，这说明相关细胞活动及细胞膜内外分子在昆虫发育过程有着重要作用，发挥各种物质转运和信号传递交流等重要生理功能（Zhaoet al.，2020）。

氧化磷酸化是昆虫生理生化过程的一个重要环节，是产生ATP 的主要方式（Dewaret al.，2022）；碳水化合物是生物体生长发育、运动的主要能量来源（Oppertet al.，2015）；脂肪酸是肌肉收缩、加快组织新陈代谢的重要能量来源，可以与磷脂结合形成细胞膜（Glatzet al.，2010）；柠檬酸（TCA）循环是氨基酸、糖类和脂类代谢的共同代谢通路，是机体通过物质代谢获取能量的最有效方式（陈牧等，2012）；CYP450药物代谢及外源物质代谢过程对于生物体有毒代谢物排出、植物次生代谢物降解中发挥重要作用（Wenet al.，2006；Maoet al.，2007；Pottieret al.，2012；Pakharukovaet al.，2015）；溶酶体在多细胞生物生长发育、结构重建、维持细胞内稳态等过程中起着重要作用（Kobayashi & Liang，2015；Vargaet al.，2015）；吞噬体具有清除凋亡细胞和入侵病原体的功能（Stuart & Ezekowitz，2008），上述代谢、解毒及免疫相关通路显著富集于幼虫阶段，说明幼虫在活动旺盛、取食频繁及接触环境增多的情况下，生物体组织代谢活动加强，解毒和免疫功能逐渐显现，以满足机体生长发育所需和免受外界侵害。随龄期增加，上调基因主要参与调控生长通路（内质网蛋白加工）（Chenet al.，2017）、与肌动蛋白合成相关的信号途径和细胞连接、分化及迁移（细胞外基质受体互作）（Alldingeret al.，2006；Peschet al.，2016）、耐药性相关通路（药物代谢-其他酶）（Santoset al.，2020）、先天性防御有关通路（FcγR 介导 的吞噬作用）（Popovicet al.，2021），这些对于幼虫生长发育过程中的组织器官功能和结构稳定性、生理生化功能及免疫相关机能完善等方面有重要作用（Liuet al.，2014）。幼虫变态为蛹，下调通路明显多于上调通路、且主要富集于代谢通路，推测可能与蛹阶段不再取食和活动，相关生理生化代谢减弱有关；上调基因则富集一些与细胞分裂有关的通路，这可能与蛹阶段细胞增殖和分化活动加剧有关，以利于成虫器官的形成（Sotillos&Celis，2005；Wonget al.，2012）。从蛹转变到成虫，上调基因显著富集通路多数与成虫飞行活动、交配、产卵及其他行为等重要功能发挥需要较多营养和能量有关；与飞行有关的昼夜节律通路在该阶段也被上调基因富集，这可能与羽化后形成的白天飞行、夜间栖息的生理节律相关。