王宇翔，楊斌，吳培福，冷堂健，邱雨，周杰瓏，*

（1. 西南林業(yè)大學生命科學學院，昆明650224；2. 西南林業(yè)大學生物多樣性保護學院，云南省森林災害預警與控制重點實驗室，昆明650224）

思茅松毛蟲Dendrolimus kikuchii（Matsumura，1927）是我國危害最嚴重的6 種松毛蟲之一，隸屬于鱗翅目Lepidoptera 枯葉蛾科Lasiocampidae 松毛蟲屬Dendrolimus，模式標本來自云南思茅（侯陶謙，1987），在我國分布于云南、四川、貴州、廣東、湖南、浙江、江蘇、安徽、福建、臺灣等省，主要危害思茅松Pinus kesiyavar.langbianensis、云南松Pinus yunnanensis、滇 油 杉Keteleeria fortunei、馬 尾 松Pinus massoniana及落葉松Larix gmelinii等針葉類植物（陳昌潔，1990；侯陶謙，1993；曹先聰，2017）。爆發(fā)成災時，針葉短時間被蠶食精光，被害樹林遠看枯黃一片，如同火燒，俗稱“無煙森林火災”；蟲害輕則影響松樹生長、落葉及枝條枯死，降低松果、木材和松脂的產(chǎn)量（侯陶謙，1987；Daiet al.，2012；Menet al.，2017），重則造成樹木死亡，給林業(yè)經(jīng)濟和生態(tài)景觀造成巨大損失（童清，劉宏屏，2009；聶中良等，2019）。松毛蟲與人接觸還會引起“松毛蟲病”（陳昌潔，1990）。

關(guān)于思茅松毛蟲的研究集中在危害調(diào)查及防治措施、腸道微生物、種群動態(tài)監(jiān)測、線粒體系統(tǒng)發(fā)育、組織病理及生理生化等方面（張榮超，2015；曹先聰，2017；Liuet al.，2017；Menet al.，2017；Wuet al.，2017；馮玉元，2018；萬鷹等，2018；李雕益等，2022；邱雨等，2022）。Zhou 等（2021）報道了思茅松毛蟲的全基因組，為該物種的系統(tǒng)進化和防治策略提供了新的視角。但迄今為止有關(guān)思茅松毛蟲基因表達及調(diào)控相關(guān)研究較為匱乏。隨著下一代測序技術(shù)（NGS）和生物信息學的發(fā)展，RNA-Seq分析成為發(fā)育生物學、進化生物學、分類學、寄主互作及毒理學等領(lǐng)域的重要研究手段（Zhanget al.，2014；Tianet al.，2015；姚順，2018；Chenet al.，2019；Norigaet al.，2019），可從整體水平上分析基因結(jié)構(gòu)和功能，有利于發(fā)現(xiàn)生物學過程的基因表達變化規(guī)律（Wang & Kirkness，2005）。昆蟲轉(zhuǎn)錄組分析為深入系統(tǒng)地研究昆蟲生命活動相關(guān)基因功能、系統(tǒng)進化與分類、生長發(fā)育以及昆蟲與相關(guān)寄主互作等提供依據(jù)（張棋麟，袁明龍，2013）?；诖耍狙芯坎捎肈NBSEQ-T7 測序平臺對思茅松毛蟲不同發(fā)育階段的轉(zhuǎn)錄組進行測序，開展各階段基因表達譜差異和功能富集分析。研究結(jié)果能豐富思茅松毛蟲遺傳資源的組學數(shù)據(jù)庫，有助于揭示基因表達與特定階段生理功能的關(guān)系，理解昆蟲全變態(tài)發(fā)育的生物學過程，為進一步闡明思茅松毛蟲生長發(fā)育的分子機制、開展功能基因表達調(diào)控及害蟲生物防治策略相關(guān)研究奠定理論基礎(chǔ)。

1 材料與方法

1.1 樣品采集

思茅松毛蟲蟲繭采自云南省安寧市縣街野外滇油杉林（102°24'E，24°46'N），采回后于室內(nèi)羽化、交配、產(chǎn)卵、孵化及飼養(yǎng)（飼養(yǎng)條件：28 ℃，相對濕度75%，光照L∶D=16∶8，幼蟲飼喂新鮮滇油杉）。收集各發(fā)育階段樣品：卵50粒、1齡幼蟲10頭、4齡幼蟲 10 頭、7 齡幼蟲 6 頭、蛹 4 頭和成蟲 4 頭。將上述樣品用PBS 溶液（pH7.4，濃度0.01 mol·L?1）清洗干凈，?80 ℃保存。

1.2 RNA提取、文庫構(gòu)建及測序

采用以Trizol 法提取總RNA 的試劑盒（TIANGEN，DP431）分別提取不同發(fā)育階段樣品的總RNA。質(zhì)檢后，使用MGIEasy RNA 文庫制備試劑套裝（深圳華大智造科技股份有限公司）建立6個cDNA文庫。通過DNBSEQ-T7（原MGISEQ-T7）測序平臺上對合格文庫進行測序。文庫構(gòu)建及測序由武漢希望組生物科技有限公司完成。轉(zhuǎn)錄組測序數(shù)據(jù)提交至NCBI_SRA 數(shù)據(jù)庫，登錄號為PRJNA732117。

1.3 數(shù)據(jù)質(zhì)控及比對參考基因組

使用 Fastp（https：//github.com/OpenGene/fastp）對原始數(shù)據(jù)進行過濾，然后利用FastQC（http：//www.bioinformatics.babraham.ac.uk/projects/fastqc）對clean data 進行質(zhì)控，質(zhì)控合格后進行后續(xù)分析。使用 Hisat2（Kimet al.，2015）將 clean reads 與參考基因組比對（DDBJ/ENA/GenBank登錄號：JAHHIN?000000000）。

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

通過定位到基因組區(qū)域或基因外顯子區(qū)的測序序列（reads）的計數(shù)來估計基因的表達水平。本研究以每百萬reads 中來自某一基因每千堿基長度的fragments 數(shù)目（FPKM）作為每個基因表達量的衡量標準，以FPKM=0.1 作為判斷基因是否表達的閾值：0.115 為高豐度表達水平。參考Bardou 等（2014）的方法繪制基因維恩圖；各發(fā)育階段基因表達量經(jīng)log2轉(zhuǎn)換后，通過R 語言繪制聚類熱圖。采用 EdgeR（Robinsonet al.，2010）對不同發(fā)育階段的DEGs 進行分析，篩選條件為|Fold Change |≥2 且 P-adjust（FDR）≤0.005。 采 用clusterProfiler（https：//github. com/GuangchuangYu/clusterProfiler；Yuet al.，2012）對 DEGs 進行 GO 和KEGG 富集分析，以FDR<0.05 作為篩選顯著富集結(jié)果的閾值，確定DEGs 主要參與的代謝通路或信號途徑。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組測序數(shù)據(jù)統(tǒng)計與質(zhì)量評估

測序獲總的數(shù)據(jù)量為369 804 176 條序列，各發(fā)育階段樣品的為47 517 160~75 263 362條；經(jīng)數(shù)據(jù)過濾質(zhì)控，約有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）。轉(zhuǎn)錄組測序質(zhì)量較高，數(shù)據(jù)可滿足生物信息學分析相關(guān)要求。

表1 不同發(fā)育階段轉(zhuǎn)錄組測序數(shù)據(jù)統(tǒng)計與質(zhì)量分析Table 1 Statistics and quality analysis of transcriptome data at different developmental stages

2.2 基因表達水平分析

各發(fā)育階段基因數(shù)量分別為11 580 個（卵）、12 128 個（1 齡幼蟲）、11 334個（4齡幼蟲）、11 209個（7 齡幼蟲）、12 160 個（蛹）和11 478 個（成蟲）?；驍?shù)量和高表達基因均為蛹階段最多，4 齡幼蟲期最少；各發(fā)育階段中，除幼蟲期外，其他3個階段均為高表達基因占比最高（表2）。

表2 不同表達水平的基因數(shù)量統(tǒng)計表Table 2 Statistics of the number of genes in different expression levels

從維恩圖可以看出（圖1：A），不同發(fā)育階段均表達的基因有9 251 個；特有基因在卵、幼蟲、蛹、成蟲階段分別有 213 個、261 個、183 個、141 個，其中幼蟲階段特有基因分別為153個（1齡幼蟲）、43個（4 齡幼蟲）、65 個（7 齡幼蟲）。熱圖顯示（圖1：B），不同發(fā)育階段呈現(xiàn)不同的基因表達聚類模式，其中4 齡幼蟲和7 齡幼蟲聚為一支，再與1 齡幼蟲匯聚，蛹和成蟲聚為一支，再與卵匯聚，最后兩大支匯合；不同發(fā)育階段的基因表達譜存在明顯差異，各發(fā)育階段均有幾個主要的高表達基因簇，卵和蛹階段尤為突出。

圖1 基因表達維恩圖（A）和基因表達模式聚類圖（B）Fig. 1 Venn diagram of gene expression（A）and cluster diagram of gene expression patterns（B）

2.3 DEGs分析

DEGs 總數(shù)為 12 927 個，其中，卵與 1 齡幼蟲的最多（3 336 個：上調(diào)1 902 個，下調(diào) 1 434 個），1 齡幼蟲與 4 齡幼蟲的有 2 398 個（上調(diào) 1 362 個，下調(diào)1 036 個），4 齡幼蟲與7 齡幼蟲的有1 758 個（上調(diào)642 個，下調(diào) 1 116 個），7 齡幼蟲與蛹的有 2 034 個（上調(diào)1 166 個，下調(diào)868 個），蛹與成蟲的最少（1 654個：上調(diào)821個，下調(diào)833個）。

不同發(fā)育階段基因表達的比較結(jié)果顯示，從卵至1齡幼蟲，顯著上調(diào)基因有肌球蛋白相關(guān)基因（XP_026493160.1 和AAV91411.1）、化學感受蛋白相關(guān)基因（AII01039.1 和AII01031.1）、視黃醇脫氫酶相關(guān)基因（XP_026737674.1）、絲氨酸蛋白酶相關(guān)基因（ACR15967.2）；而與有絲分裂特異性細胞周期蛋白cyclin-B1 相關(guān)基因（XP_026731829.1）、細胞分裂周期蛋白20（CDC20）相關(guān)基因（XP_026733697.1）明顯下調(diào)。4 齡幼蟲相對1 齡幼蟲，載脂蛋白相關(guān)基因（XP_026760716.1）、細胞壁蛋白相關(guān)基因（XP_004932399.1）、表皮蛋白相關(guān)基因（XP_026734677.1）、保幼激素相關(guān)基因（XP_004925704.1 和 XP_026743329.1）、蛻皮激素相關(guān)基因（XP_026761583.1 和XP_026735381.1）等顯著上調(diào)；7 齡幼蟲相對4 齡幼蟲，解毒酶細胞色素P450 相關(guān)基因（XP_021191965.1 和 AID54859.1）、表皮蛋白相關(guān)基因（XP_026762046.1）、保幼激素相關(guān)基因（XP_026733514.1和XP_026743329.1）、蛻皮相關(guān)基因（XP_026330046.1和XP_026761583.1）等表達上調(diào)，而脂肪酶相關(guān)基因（AEA76291.1）表達下調(diào)。從7齡幼蟲到蛹階段，顯著上調(diào)基因有溶菌酶相關(guān)基因（XP_026321249.1）、脂肪酸合成酶相關(guān)基因（XP_026499944.1）、過氧化物酶相關(guān)基因（XP_026725384）、長鏈酯酰輔酶A 合成酶相關(guān)基因（XP_026724857.1）、儲存蛋白相關(guān)基因（XP_026327069.1）、表皮蛋白相關(guān)基因（XP_026501352.1）、卵黃原蛋白基因（BAB16412.1）等，顯著下調(diào)基因有細胞色素P450 家族相關(guān)基因（AID54859.1 和（XP_021191965.1）。從蛹到成蟲階段，角質(zhì)層蛋白相關(guān)基因（XP_026756401.1）、肌鈣蛋白 C 相關(guān)基因（XP_026757608.1）、脂肪酸合成酶相關(guān)基因（XP_026499944.1）、化學感受蛋白相關(guān)基因（AII01034.1）、飛行有關(guān)基因（XP_026317122.1和XP_026757413.1）等顯著上調(diào)，而與血淋巴蛋白酶相關(guān)基因（AAV91004.1）、溶酶體相關(guān)基因（XP_026734095.1）、Toll樣受體相關(guān)基因（XP_026744547.1）表達下調(diào)。

2.4 不同發(fā)育階段DEGs功能富集分析

GO分析表明，各階段DEGs均能顯著富集到生物過程、細胞組分及分子功能，其中涉及分子功能的條目最多、細胞組分的最少（圖2）。KEGG 分析發(fā)現(xiàn)，各階段DEGs顯著富集到代謝、生物體系統(tǒng)、遺傳信息處理、細胞過程、環(huán)境信息處理5類途徑，涉及101條通路，其中代謝途徑最多（57條通路；圖3）。

圖2 不同發(fā)育階段差異表達基因GO富集分析Fig. 2 GO enrichment analysis of differentially expressed genes at different developmental stages

圖3 不同發(fā)育階段差異表達基因KEGG通路富集分析Fig. 3 KEGG pathway enrichment analysis of differentially expressed genes at different developmental stages

從卵到1 齡幼蟲，DEGs 主要富集于催化活性和代謝過程。相比于1齡幼蟲，卵期上調(diào)基因主要參與細胞分裂過程和蛋白質(zhì)合成過程。隨齡期增加，GO 富集情況與上、下調(diào)基因參與KEGG 通路均有所變化，主要體現(xiàn)在免疫功能和解毒代謝途徑上。在幼蟲變態(tài)為蛹階段，DEGs 顯著富集于催化活性、結(jié)構(gòu)組分和結(jié)合分子功能，氧化還原過程、跨膜運輸和蛋白質(zhì)水解（生物過程），核糖體、胞外區(qū)域（細胞組分）；上調(diào)基因富集于DNA 復制、核苷酸切除修復，下調(diào)通路主要有核糖體和大量代謝途徑。當蛹變?yōu)槌上x，DEGs富集于跨膜運輸、氧化還原過程、催化活性、細胞粘附及蛋白質(zhì)折疊的細胞增殖過程；上調(diào)基因參與氨基酸代謝、脂肪酸代謝和與飛行有關(guān)的晝夜節(jié)律等途徑。

3 討論

通過對思茅松毛蟲發(fā)育過程的轉(zhuǎn)錄組學研究和差異表達基因的分子解析，可為后續(xù)開展發(fā)育生物學研究奠定分子基礎(chǔ)，為今后害蟲綜合防治提供新的方向和策略。本文采用DNBSEQ-T7測序平臺對思茅松毛蟲進行轉(zhuǎn)錄組測序分析，發(fā)現(xiàn)各發(fā)育階段差異表達基因與每個階段特定生理活動存在顯著相關(guān)性。在卵期，參與細胞周期和細胞分裂相關(guān)基因有較高表達。幼蟲階段，與消化酶、代謝、肌肉生長、感覺等相關(guān)基因明顯上調(diào)；隨齡期增加，保幼激素、解毒、表皮蛋白及抗性等相關(guān)基因也顯著上調(diào)。蛹階段，與細胞增殖和分化相關(guān)基因被激活，脂肪酸合成酶、儲存蛋白、表皮蛋白、溶菌酶等相關(guān)基因顯著上調(diào)，這些可能與保證蛹發(fā)育和角質(zhì)層形成提供足夠蛋白質(zhì)和氨基酸（Pick&Burmester，2009；Allenet al.，2018）、由凋亡和自噬介導對幼蟲器官的自我消化（Tashiroet al.，1976；Joneset al.，1995）和免疫防御機制（Chung &Ourth，2000；Xieet al.，2013）等有關(guān)。成蟲階段，高表達基因主要與組織器官形成、外部形態(tài)、繁殖和飛行等過程顯著相關(guān)。此外，基因表達具有發(fā)育時空特異性，卵與1齡幼蟲間差異表達基因數(shù)目最多，說明從卵向幼蟲的轉(zhuǎn)變是最為復雜的過程，這與戚禮興（2016）、Yang 等（2016）及 Han 等（2019）在菜粉蝶Pieris rapae、馬尾松毛蟲Dendrolimus punctatus和云南松毛蟲Dendrolimus houi上的研究結(jié)果一致。

在GO 功能富集中，各發(fā)育階段的DEGs 顯著富集到了氧化還原、蛋白質(zhì)水解、跨膜運輸、碳水化合物代謝等生物過程，這些過程為機體生長發(fā)育、運動提供營養(yǎng)和能量來源（Claudiuset al.，2014；Massoet al.，2016；Gautamet al.，2020）。在分子功能中，較多地富集在活性功能和結(jié)合功能上。其中，催化活性最為典型，生物體各類催化活動都離不開酶的參與，保證了昆蟲物質(zhì)代謝和營養(yǎng)能量轉(zhuǎn)換等過程正常進行（趙星，2021）。結(jié)合功能與一些分子和關(guān)鍵因子/離子需要結(jié)合后才能變?yōu)榛钚?、有功能的調(diào)控有關(guān)。結(jié)合功能和催化活性功能顯著富集情況與姚順（2018）的研究結(jié)果較為一致。幾丁質(zhì)結(jié)合是唯一一個在各發(fā)育階段比較中均被顯著富集的條目，這可能與其全變態(tài)發(fā)育的周期性表皮降解和重構(gòu)等生理過程有關(guān)（Co?hen，2001）。在細胞組分中，顯著富集到胞外區(qū)域、線粒體內(nèi)膜、質(zhì)膜等，這說明相關(guān)細胞活動及細胞膜內(nèi)外分子在昆蟲發(fā)育過程有著重要作用，發(fā)揮各種物質(zhì)轉(zhuǎn)運和信號傳遞交流等重要生理功能（Zhaoet al.，2020）。

氧化磷酸化是昆蟲生理生化過程的一個重要環(huán)節(jié)，是產(chǎn)生ATP 的主要方式（Dewaret al.，2022）；碳水化合物是生物體生長發(fā)育、運動的主要能量來源（Oppertet al.，2015）；脂肪酸是肌肉收縮、加快組織新陳代謝的重要能量來源，可以與磷脂結(jié)合形成細胞膜（Glatzet al.，2010）；檸檬酸（TCA）循環(huán)是氨基酸、糖類和脂類代謝的共同代謝通路，是機體通過物質(zhì)代謝獲取能量的最有效方式（陳牧等，2012）；CYP450藥物代謝及外源物質(zhì)代謝過程對于生物體有毒代謝物排出、植物次生代謝物降解中發(fā)揮重要作用（Wenet al.，2006；Maoet al.，2007；Pottieret al.，2012；Pakharukovaet al.，2015）；溶酶體在多細胞生物生長發(fā)育、結(jié)構(gòu)重建、維持細胞內(nèi)穩(wěn)態(tài)等過程中起著重要作用（Kobayashi & Liang，2015；Vargaet al.，2015）；吞噬體具有清除凋亡細胞和入侵病原體的功能（Stuart & Ezekowitz，2008），上述代謝、解毒及免疫相關(guān)通路顯著富集于幼蟲階段，說明幼蟲在活動旺盛、取食頻繁及接觸環(huán)境增多的情況下，生物體組織代謝活動加強，解毒和免疫功能逐漸顯現(xiàn)，以滿足機體生長發(fā)育所需和免受外界侵害。隨齡期增加，上調(diào)基因主要參與調(diào)控生長通路（內(nèi)質(zhì)網(wǎng)蛋白加工）（Chenet al.，2017）、與肌動蛋白合成相關(guān)的信號途徑和細胞連接、分化及遷移（細胞外基質(zhì)受體互作）（Alldingeret al.，2006；Peschet al.，2016）、耐藥性相關(guān)通路（藥物代謝-其他酶）（Santoset al.，2020）、先天性防御有關(guān)通路（FcγR 介導 的吞噬作用）（Popovicet al.，2021），這些對于幼蟲生長發(fā)育過程中的組織器官功能和結(jié)構(gòu)穩(wěn)定性、生理生化功能及免疫相關(guān)機能完善等方面有重要作用（Liuet al.，2014）。幼蟲變態(tài)為蛹，下調(diào)通路明顯多于上調(diào)通路、且主要富集于代謝通路，推測可能與蛹階段不再取食和活動，相關(guān)生理生化代謝減弱有關(guān)；上調(diào)基因則富集一些與細胞分裂有關(guān)的通路，這可能與蛹階段細胞增殖和分化活動加劇有關(guān)，以利于成蟲器官的形成（Sotillos&Celis，2005；Wonget al.，2012）。從蛹轉(zhuǎn)變到成蟲，上調(diào)基因顯著富集通路多數(shù)與成蟲飛行活動、交配、產(chǎn)卵及其他行為等重要功能發(fā)揮需要較多營養(yǎng)和能量有關(guān)；與飛行有關(guān)的晝夜節(jié)律通路在該階段也被上調(diào)基因富集，這可能與羽化后形成的白天飛行、夜間棲息的生理節(jié)律相關(guān)。