野桑蠶蛹滯育時期的轉錄組測序分析與蛹滯育關聯基因篩選

摘 要 旨在篩選野桑蠶蛹滯育的調控關聯基因，為在分子水平揭示該昆蟲蛹滯育調控機制提供理論基礎。以野桑蠶蛹滯育不同時期的全蛹樣本為研究對象，采用高通量測序平臺Illumina HiSeqTM 2000進行轉錄組測序分析；利用KAAS在線pathway比對分析工具對篩選出的滿足padj lt; 0.05且fold change ≥ 32或padj lt; 0.05且fold change ≤ 1/32條件的差異表達基因進行KEGG通路富集分析；利用stem軟件對差異表達基因進行時間序列分析。選取前滯育期和滯育解除期低表達，滯育時期高表達的基因為蛹滯育關聯基因，并進行GO和KEGG富集分析。結果表明，在滯育前期、滯育中期、滯育晚期和滯育解除期分別檢測到13 764、" 13 765、13 765和13 764個表達轉錄本。樣本組的兩兩比較發現差異表達基因4 621 個。與滯育前期相比，滯育中期、滯育晚期和解除滯育期的上調表達基因分別是1 921、451和38 個，下調表達的基因分別是1 143、171和59個。利用時間序列分析方法在差異表達基因中發現290 個蛹滯育關聯基因，GO富集分析將202個基因映射到52項二級GO類目，KEGG通路富集分析將67 個基因映射到126 個三級代謝通路。

關鍵詞 野桑蠶；轉錄組；時間序列分析；蛹滯育關聯基因

昆蟲滯育指昆蟲在特定環境刺激下發育停滯，在另一種特定條件刺激后繼續發育的現象，是昆蟲適應不良自然環境而獲得的一種生存本領[1]。滯育按發生所處的發育階段，可分為胚胎滯育、幼蟲滯育、蛹滯育和成蟲滯育等。蛹滯育在鱗翅目和雙翅目昆蟲中普遍存在。蛹滯育調控機制復雜，目前分析鑒定的蛹滯育調控相關基因包括熱激蛋白[2-3]、能量儲存與代謝[4-6]、激素代謝及其調控[7-9]、表觀遺傳[10-11]、生物鐘系統[12]、光信號傳導[13]、關鍵性的細胞信號轉導通路[14-16]等方面的基因或作用因子。由于物種和滯育蟲態及其生理途徑不盡相同，滯育控制機制差別較大，因此對于昆蟲滯育機制的解析或對不同昆蟲的蛹滯育調控異同尚需深入了解[17]。

野桑蠶（Bombyx mandarina）是鱗翅目（Lepidoptera）蠶蛾科（Bombycidae）昆蟲，在東亞、南亞廣泛分布。野桑蠶具有豐富的遺傳多樣性，可以和家蠶雜交并產生正常子代，因而在家蠶新品種創制和分子改良方面具有重要的應用價值。Shen等[18]比較了中國（蘇州）和日本（東京）野桑蠶蛹發育情況，結果表明中國野桑蠶蛹期為10至40 d，而日本野桑蠶在第 1 代中存在夏滯育（summer diapause）現象，即少量個體可維持蟲蛹狀態達3個月而不化蛾。進一步的研究表明，日本野桑蠶在蛹滯育期間卵巢發育停滯，呼吸耗氧量、蛻皮激素含量均維持在低水平，深入探究發現野桑蠶蛹滯育的形成可能與腦組織中的促前胸腺激素（PTTH）合成受阻有關[19-21]。筆者此前的研究表明，秦巴山區野桑蠶是一化、二化和多化性的混合群體，第一代于五月中下旬化蛹，單蛾圈中蛹期時長大于30 d的個體數量可占總數的90%以上，少量蛹個體可遲至9-10月羽化，期間即便改變環境條件也不會即刻恢復變態發育，具有蛹滯育現象，但其關鍵性的控制因子及調控機制尚不明確。

滯育關聯蛋白（Diapause Associated Protein，DAP）是在昆蟲滯育期間特異性表達、在非滯育階段不表達或微量表達的一類蛋白，通常是儲藏蛋白、抗凍蛋白、熱休克蛋白、分子伴侶或重要的代謝調節酶類物質，在昆蟲滯育適應、抵抗逆境脅迫等生化途徑中發揮重要作用[22-23]。目前，柞蠶[24]、棉鈴蟲[25]、七星瓢蟲[26]、美國白蛾[23]等的滯育關聯蛋白有所報道。有研究認為，滯育關聯蛋白的調控與中樞神經內分泌系統以及相應的代謝變化有關。在滯育期間，較高的保幼激素和較低的蛻皮激素水平，可促使滯育關聯蛋白在脂肪體中大量累積，并在整個滯育過程中維持在高水平，維系著蟲蛹的生存。隨著滯育的結束，保幼激素水平下降，滯育關聯蛋白隨之逐漸消失[27]。與之相對應的滯育關聯基因，即在滯育期間特異性上調、在滯育解除后下調的一類基因，在解析滯育調控機制方面有著重要作用。

近年來，組學技術在蔥蠅[28]、吉氏蚜蟲[29]、柑橘大實蠅[30]、菜粉蝶[31]等昆蟲蛹滯育的研究中得到廣泛運用，為發現滯育關聯基因提供了快捷途徑。分析野桑蠶蛹滯育調控機制可以為創制長蛹期家蠶新品種提供基因資源和理論依據。本研究對野桑蠶蛹滯育不同時期的材料開展轉錄組測序，分析不同階段的全基因組基因表達情況及差異表達基因，利用時間序列分析篩選蛹滯育關聯基因，為進一步研究野桑蠶蛹滯育的調控基因以及分子機理奠定基礎。

1 材料與方法

1.1 昆蟲材料

野桑蠶為陜西省蠶桑重點實驗室提供。越年卵經低溫冷藏后于四月上旬催青，收蟻后將野桑蠶飼養于自然條件下的飼育大棚中，五月中下旬營繭化蛹，收集野桑蠶蛹于室內25" ℃保護，野桑蠶蛹期20 d～120 d。按照體表、眼點、翅脈等器官組織的著色特征，選取化蛹后第3天、第20天、第30天和第60天等4個時期的蠶蛹樣本，分別標記為前滯育期（preD）、滯育中期（EarD）、滯育晚期（LatD）和解除滯育期（PosD）。每個時期組設置3個重復，每個重復中包含6只雌蛹。蛹體表面經超純水清潔之后液氮速凍，-80 ℃超低溫冰箱中保存，備用。

1.2 總RNA提取

使用DEPC處理后的1×PBS緩沖液清洗樣本，液氮中研磨至粉末狀。使用TRIzol Reagent試劑盒（Invitrogen，美國）提取總RNA。RNA質量檢測的主要方法：（1）瓊脂糖凝膠電泳分析RNA降解程度以及是否有污染；（2）NanoDrop檢測RNA的純度（OD260/OD280比值）；（3）Qubit對RNA濃度進行精確定量；（4）Agilent 2100精確檢測RNA的完整性。經檢測，提取的RNA純度均滿足OD值（260/280）為1.8～2.2，濃度大于200 ng/μL，總含量大于3.5 μg的測序建庫要求。將每組的6個RNA樣品等濃度混合后用于構建測序文庫，共形成4組12個樣本。

1.3 構建RNA-seq文庫和Illumina高通量測序

各樣品總RNA檢測合格后，用帶有Oligo（dT）的磁珠富集真核mRNA，然后加入片段緩沖液隨機打斷mRNA。以mRNA為模板，用隨機引物合成cDNA第一鏈，然后加入緩沖液、dNTPs、RNase H和DNA polymerase I合成cDNA第二鏈，最后用AMPure XP beads純化cDNA。純化的雙鏈cDNA進行末端修復、加A尾并連接測序接頭，然后用AMPure XP beads進行片段大小選擇，最后PCR擴增富集得到cDNA文庫。文庫質檢合格之后使用Illumina Hiseq"" 2 500測序平臺進行高通量測序。構建文庫和測序由北京百邁克生物科技有限公司完成。轉錄組測序原始數據已上傳至GenBank SRA數據庫，接收編號為SRR7811946、SRR7811953、SRR7811952、SRR7811949、SRR7811948、SRR7811951、SRR7811950、SRR7811945、SRR7811944、SRR7811955、SRR7811954、SRR7811947。

1.4 RNA-seq測序數據處理和樣本相關性分析

原始數據（raw data）經過濾去除帶接頭、低質量的讀序（reads）、核糖體序列、N率較高序列和長度過短序列后獲得高質量測序數據。利用TopHat2軟件將高質量測序數據比對至野桑蠶參考基因組，并進行轉錄組注釋和表達量的計算。參考基因組來源于筆者前期利用Oxford Nanopore、Illumina NovaSeq和Illumina HiSeq X測序組裝獲得的秦巴野桑蠶基因組，基因組大小為457.85 Mb（數據待發表）。基于各樣本表達量結果，利用R（http：//www.r-project.org）開展主成分分析（Principal Component Analysis，PCA），計算各樣本之間的聚類關系。

1.5 差異表達基因篩選與功能富集分析

采用FPKM方法計算基因的表達量[32]。利用DEGseq軟件包篩選差異表達基因（DEG），篩選閾值為Fold change gt; 2和Q-value lt; 0.05，Q-value越小，基因表達量差異越顯著[33]。使用Blast2GO軟件（version 2.5）（https：//www.geneontology.org/）進行GO基因本體功能注釋分析。經校正，Q-valuelt;0.05的GO條目被認為在DEGs中顯著富集。使用KAAS軟件（r140224）進行KEGG代謝通路富集分析。經校正，Q-valuelt;0.05的KEGG通路被認為在DEGs中顯著富集。

1.6 時間序列分析和滯育關聯基因篩選

利用STEM[34]軟件進行時間序列分析。輸入按生物學邏輯順序排列的分組樣本的全部差異基因表達量文件，選擇參數（-pro 20 -ratio 1.000 0）進行時間序列分析，然后對獲得的各趨勢中的基因進行GO和KEGG功能富集分析，并通過假設檢驗計算p value。通過FDR校正p value后，以Q-value ≤ 0.05為閾值，滿足此條件的GO和Pathway分類，被定義為在該趨勢中顯著富集的GO term和Pathway。依據文獻報道并通過對比篩選，將與非蛹滯育期樣品相比，滯育期樣品中上調且在滯育解除期樣品中下調的基因認定為滯育關聯基因[24，26，31]。

2 結果與分析

2.1 測序結果的組裝與序列比對分析

利用高通量測序平臺對12個樣本進行RNA-seq測序，共獲得103.3 Gb的原始數據（raw data），經質量控制后獲得clean data 102.03 Gb，各樣本clean data為7.49 Gb～9.17 Gb，GC含量為44.19%～47.53%，Q30均在" 91.84%以上，說明測序質量可靠。將高質量讀序與參考基因組進行比對，結果顯示各樣本讀序與基因組的匹配度為75.82%～83.31%，其中唯一的匹配占比為74.62%～79.92%，多重匹配占比為" 3.03%～4.20%（表1）。

2.2 表達量統計和樣本關系分析

歸一化處理后，將FPKM值小于0.5的轉錄組本去除，獲得各時期樣本中存在表達的轉錄本數量。結果表明，野桑蠶在蛹滯育前期（preD）、滯育早期（EarD）、滯育末期（LatD）和滯育解除期（PosD）可檢測到的表達轉錄本分別是13 764、137 65、13 765和13 764個，其中無表達（0～" 0.5）、低表達（≥0.5-1）、中表達（≥1-10）和高表達（≥10）轉錄本分別為1 827～2 470、652～703、3 439～3 778和7 185～7 501個（圖1）。

基于各樣本表達量結果，通過PCA分析和計算樣本間皮爾斯相關系數，了解樣本之間的重復性情況。在PCA分析中，樣品組分越相似，反映在PCA圖中的距離越近，來自不同有效處理的樣品往往表現出各自聚集的分布。結果表明，滯育早期和滯育晚期的樣品組各自形成聚類，且與其他樣本組保持著較遠距離，但滯育前（preD）和滯育解除后（posD）的樣品組聚類距離極為接近（圖2）。

2.3 差異表達基因篩選和功能富集分析

在蛹發育4個時期的樣本中共發現差異表達基因4 621個（圖3）。與滯育前期相比，滯育中期、滯育晚期和解除滯育期的上調表達基因分別是1 921、451和38個，下調表達的基因分別是" 1 143、171和59個；與滯育中期相比，滯育后期、解除滯育期上調表達的基因分別是1 230和" 1 378個，下調表達的基因為1 786和2 141個；與滯育后期相比，滯育解除期上調和下調表達的基因分別是214和520個。差異表達基因的功能富集分析表明，3 321個基因獲得GO注釋，其中分子功能（Molecular function）、生物過程（Biologic process）和細胞組分（Cellular component）分別富集3 162、2 921和 2 341個基因（圖4）。

2.4 差異表達基因的時間序列分析與蛹滯育關聯基因篩選

將4 621個差異表達基因進行時間序列分析，獲得20組基因表達模式，其中前6組（FDRlt;0.05）富集基因分別為1 477、1 061、394、309、290和278個。模式16中的290個基因表現為在前滯育期、滯育解除期處于較低表達水平，滯育前期、滯育后期處于高表達水平，是本研究中篩選獲得的蛹滯育關聯基因（圖5）。

202個蛹滯育關聯基因獲得GO功能注釋（圖6）。生物學過程中，富集基因數量較多的是單細胞過程（Single-organism process）、細胞過程（Cellular process）和代謝過程（Metabolic process）等，而生物黏連（Biological adhesion）、節律過程（Rhythmic process）和突觸傳遞中的突觸前過程（Presynaptic process involved in synaptic transmission）等富集基因數量較少；細胞組分中，基因富集數量較高的包括細胞（Cell）、細胞組分（Cell part）和細胞膜（Membrane）等，超分子纖維（Supramolecular fiber）、病毒粒子（Virion）和病毒粒子組分（Virion part）等富集數量較少；分子功能中，催化活性（Catalytic activity）、偶聯（Binding）和轉運活性（Transporter activity）等類目富集的基因數量較高，電子載體活性（Electron carrier activity）、轉錄因子活性（Transcription factor activity），蛋白偶聯（Protein binding）和翻譯調控活性（Translation regulator activity）等基因富集數量較少。

2.5 滯育關聯基因的KEGG通路分析

對野桑蠶蛹滯育關聯基因進行KEGG富集分析，90個蛹滯育關聯基因獲得KEGG功能注釋（圖7，表2），新陳代謝（Metabolism）、遺傳信息處理（Genetic information processing）、環境信息處理（Environmental information processing）、細胞過程（Cellular processes）、有機體系統（Organismal Systems）等5個系統分別富集39、3、13、16和24個基因。

新陳代謝系統中，滯育關聯基因共參與9個代謝途徑、42個通路，分別有15、12、7、6、6、4、4、2和1個基因參與脂質、碳水化合物、氨基酸、甘氨酸合成與分解、輔因子和維生素代謝、其他氨基酸代謝、外源物生物降解和代謝、核苷酸、能量等代謝通路。遺傳信息處理系統中，滯育關聯基因涉及2個途徑，包含核糖體、泛素介導的蛋白質降解2個代謝通路，其中前者富集了核糖體蛋白大亞基基因（rpl1），后者富集了泛素結合酶E2R（UBE2R）基因和泛素結合酶E2Q（UBE2Q）基因。環境信息處理系統中，13個滯育關聯基因共參與3個途徑，涉及13個信號通路，其中參與FOXO、AMPK、PI3K-Akt、mTOR信號通路的基因有磷酸烯醇丙酮酸羧激酶（pepck）和3-磷酸肌醇依賴性蛋白激酶1（ pdk1）。細胞過程系統中，16個滯育關聯基因共參與運輸和分解代謝、細胞生長和死亡、真核生物細胞通信和細胞遷移等4個途徑，涉及溶酶體、過氧化物酶體、吞噬體、動物體中的自噬、壞死樣凋亡、細胞凋亡、緊密連接、黏著斑、細胞骨架調節等9個通路。有機體系統中， 24個滯育關聯基因參與33個通路，其中一氧化氮合酶基因（nsl）參與晝夜節律調控，生醇脂肪酰輔酶A還原酶基因（wat）參與長壽調控途徑，磷酸烯醇丙酮酸羧激酶（pepck）和3-磷酸肌醇依賴性蛋白激酶1參與胰島素信號通路，PPAR信號通路則富集了酰基輔酶A氧化酶基因（CG5009）、磷酸烯醇丙酮酸羧激酶（pepck）、3-磷酸肌醇依賴性蛋白激酶1（ pdk1）、脂肪轉運蛋白1/4（ slc27a4）和脂肪酸結合蛋白3（ fabp3）等基因。

3 討" 論

野桑蠶是家蠶的祖先種。蛹期是蠶的馴化性狀之一，現行家蠶品種蛹期均在10 d左右，而在中國、日本的野桑蠶中均有蛹滯育導致的長蛹期現象，可達20多天直至數月不化蛾。解析滯育分子調控機制，可以為病蟲害防治以及家蠶新種創制和家蠶遺傳改良提供線索。本研究利用轉錄組測序技術分析了野桑蠶蛹滯育不同時期的基因表達情況，并對野桑蠶蛹滯育關聯基因進行了篩選與分析。

時間序列分析是利用多個連續型樣本對基因的表達模式進行聚類分析的方法。有研究表明，蛹滯育不同階段存在著特異性的調控機制，時間序列分析可以帶來階段特異性的生物學信息。在野桑蠶蛹4個不同時期中共發現4 621個差異表達基因，通過時間序列分析獲得20個表達模式聚類，表明野桑蠶蛹滯育存在復雜的調控機制，其中模式16中的基因表達模式為滯育期間表達上調而滯育解除后迅速下調，是本研究中篩選出的蛹滯育關聯基因。值得注意的是，模式3中的278個基因表達情況與模式16的情況相反。比較二者KEGG代謝注釋結果異同，結果發現模式3中有111個基因富集于遺傳信息處理系統，顯著高于模式16中的情形，這些基因涉及翻譯、轉錄、DNA修復、蛋白折疊和運輸等途徑，表明在野桑蠶蛹滯育過程中，DNA、RNA和蛋白質的合成、降解過程均受到顯著抑制，蟲體生理活動處于較低水平，這與其他昆蟲中的研究及野桑蠶滯育蛹的表現形態是一致的，相關基因可能受到激素的調控，而其調控機制值得深入研究。此外，模式0和模式19的基因蔟分別呈持續下調或上調，可能與滯育維持和重新開始發育有關，而模式14/5、13/6、11/8、17/2等分別在滯育的前期或晚期出現基因的上調或下調，可能與蛹滯育不同階段的發育狀態轉換相關聯。

研究表明，昆蟲滯育受內分泌激素調控，促前胸腺激素（PTTH）受體信號轉導途徑、保幼激素和蛻皮激素合成途徑均與滯育密切相關[35]。此外，研究還表明多種細胞信號轉導途徑，包括MAPK信號通路、Wnt信號通路、mTOR信號通路、鈣信號通路、類固醇激素合成途徑、FoxO信號途徑以及胰島素信號等通路相關基因[35]，以及光周期[36]、晝夜節律[37]等生態因素參與昆蟲蛹滯育調控，野桑蠶蛹滯育關聯基因中也發現與此相關的多個功能基因，一方面表明野桑蠶蛹滯育調控的復雜性，另一方面也暗示著野桑蠶中存在與棉鈴蟲、柑橘鳳蝶等昆蟲相似的調控機制，可以為昆蟲蛹滯育機制研究提供參考。值得注意的是，根據文獻報道，蛹滯育進入、蛹滯育解除等不同階段均存在不同的調控機制[31]，亦有研究表明，野桑蠶蛹滯育的形成與其頭部的激素分泌有關[8]。本研究對全蛹進行轉錄組測序分析，今后還需要對特定器官特定發育時期開展基因表達分析，可能有助于獲得蛹滯育形成與調控的新發現。滯育不同時期的基因差異表達比較發現，前滯育期和滯育解除期之間的差異基因數量顯著少于其他組別，這與已有的研究存在一定的一致性[31]，可能是與野桑蠶夏滯育期間發育相對緩慢甚至停滯、基因表達受到抑制有關。

參考文獻 Reference：

［1］TAUBER M，TAUBER C.Insect seasonality：diapause maintenance，termination，and postdiapause [J].Development，1976，21（1）：81-107.

[2]ZHANG J，MIANO F N，JIANG T，et al.Characterization of three heat shock protein genes in Pieris melete and their expression patterns in response to temperature stress and pupal diapause [J]. Insects，2022，13（5）：430.

[3]ZHAO J，HUANG Q，ZHANG G，et al.Characterization of two small heat shock protein genes （Hsp17.4 and Hs20.3） from Sitodiplosis mosellana，and their expression regulation during diapause [J].Insects，2021，12（2）：119.

[4]ZHAO L，WANG X，LIU Z，et al.Energy consumption and cold hardiness of diapausing fall webworm pupae[J].Insects， 2022，13（9）：853.

[5]LI Y N，LIU Y B，XIE X Q，et al.The modulation of trehalose metabolism by 20-hydroxyecdysone in antheraea pernyi （lepidoptera：saturniidae） during its diapause termination and post-termination period [J].Journal of Insect Science，2020，20（5）：18.

[6]CHEN C，MAHAR R，MERRITT M E，et al.ROS and hypoxia signaling regulate periodic metabolic arousal during insect dormancy to coordinate glucose，amino acid，and lipid metabolism [J].Proceedings of the National Academy of Sciences of the United States of America，2021，118（1）：e2017603118.

[7]CHEN Z，DONG Y，WANG Y，et al.Pupal diapause termination in Bactrocera minax：an insight on 20-hydroxyecdysone induced phenotypic and genotypic expressions [J].Scientific Reports，2016，6：27440.

[8]WANG Q，LUO YT，WANG Y，et al.Expression patterns of three important hormone genes and respiratory metabolism in antheraea pernyi during pupal diapause under a long photoperiod [J].Insects，2021，12（8）：699.

[9]XU QY，DENG P，ZHANG Q，et al.Ecdysone receptor isoforms play distinct roles in larval-pupal-adult transition in Leptinotarsa decemlineata [J].Insect Science，2020，" 27（3）：487-99.

[10] REYNOLDS J A，BAUTISTA-JIMENEZ R，DENLINGER D L.Changes in histone acetylation as potential mediators of pupal diapause in the flesh fly，Sarcophaga bullata [J].Insect Biochemistry and Molecular Biology，2016，76：29-37.

[11] PEGORARO M，BAFNA A，DAVIES N J，et al.DNA methylation changes induced by long and short photoperiods in Nasonia [J].Genome research，2016，26（2）：203-210.

[12]ZHU L，TIAN Z，GUO S，et al.Circadian clock genes link photoperiodic signals to lipid accumulation during diapause preparation in the diapause-destined female cabbage beetles Colaphellus bowringi [J].Insect Biochemistry and Molecular Biology， 2019，104：1-10.

[13]ZHU L，LIU W，TAN Q Q，et al.Differential expression of circadian clock genes in two strains of beetles reveals candidates related to photoperiodic induction of summer diapause [J].Gene，2017，603：9-14.

[14]LEE I H，DUVALL L B.Maternally instigated diapause in Aedes albopictus：coordinating experience and internal state for survival in variable environments [J].Frontiers in Behavioral Neuroscience，2022，16：778264.

[15]WANG T，GENG S L，GUAN Y M，et al.Deacetylation of metabolic enzymes by Sirt2 modulates pyruvate homeostasis to extend insect lifespan [J].Aging （Albany NY），2018，10（5）：1053-1072.

[16]ZHANG X S，WANG T，LIN X W，et al.Reactive oxygen species extend insect life span using components of the insulin-signaling pathway [J].Proceedings of the National Academy of Sciences of the United States of America，2017，114（37）：E7832-E40.

[17]HUTFILZ C.Endocrine regulation of lifespan in insect diapause [J].Frontiers in Physiology，2022，13：825057.

[18]SHEN W，HAMANO K，FUGO H.Role of the brain on the resumption of development in summer diapausing pupae of Bombyx mandarina [J].Proceedings of the Japan Academy，Series B，1993，69（6）：149-53.

[19]SHEN W，HAMANO K，FUGO H.The relationships between pupal development and ecdysteroid levels in summer diapausing silkworm，Bombyx mandarina [J].Proceedings of the Japan Academy Ser B Physical and Biological Sciences，1993，69（6）：139-143.

[20]SHEN W，HAMANO K.The relationships between oxygen consumption and duration of pupal-adult development in the silkworm Bombyx mandarina [J].Proceedings of the Japan Academy Series B，Physical and Biological Sciences，1993，69（6）：134-138.

[21]SHEN W，HAMANO K，FUGO H.Relationships between an application of KK-42 and pupal development in the Silkworm，Bombyx mandarina [J].Proceedings of the Japan Academy，Series B，1993，69（6）：144-148.

[22]BI Z，YANG X，YU W，et al.Diapause-associated protein3 functions as Cu/Zn superoxide dismutase in the Chinese oak silkworm （Antheraea pernyi） [J].PLoS One，2014，" 9（3）：e90435.

[23]DENG Y，LI F，RIESKE L K，et al.Transcriptome sequencing for identification of diapause-associated genes in fall webworm，Hyphantria cunea Drury [J].Gene，2018，668：229-36.

[24]DU J，ZHAO P，WANG J，et al.Pupal diapause termination and transcriptional response of antheraea pernyi （Lepidoptera：saturniidae） triggered by 20-hydroxyecdysone [J].Frontiers in Physiology，2022，13：888643.

[25]GENG S L，ZHANG X S，XU W H.COXIV and SIRT2-mediated G6PD deacetylation modulate ROS homeostasis to extend pupal lifespan [J].The FEBS Journal，2021，288（7）：2437-2453.

[26]QI X，ZHANG L，HAN Y，et al.De novo transcriptome sequencing and analysis of Coccinella septempunctata L. in non-diapause，diapause and diapause-terminated states to identify diapause-associated genes [J].BMC Genomics，2015，16：1086.

[27]FAN L，LIN J，ZHONG Y，et al.Shotgun proteomic analysis on the diapause and non-diapause eggs of domesticated silkworm Bombyx mori [J].PLoS One，2013，8（4）：e60386.

[28]REN S，HAO Y J，CHEN B，et al.Global transcriptome sequencing reveals molecular profiles of summer diapause induction stage of onion maggot，Delia antiqua （Diptera：Anthomyiidae） [J].G3-Genes Genomes Genetics，2018，" 8（1）：207-217.

[29]ZHANG H Z，LI YY，AN T，et al.Comparative transcriptome and iTRAQ proteome analyses reveal the mechanisms of diapause in Aphidius gifuensis Ashmead （Hymenoptera：Aphidiidae） [J].Frontiers in Physiology，2018，9：1697.

[30]WANG J，RAN L L，LI Y，et al.Comparative proteomics provides insights into diapause program of Bactrocera minax （Diptera：Tephritidae [J]PLoS One，2020，15（12）：e0244493.

[31]PRUISSCHER P，LEHMANN P，NYLIN S，et al.Extensive transcriptomic profiling of pupal diapause in a butterfly reveals a dynamic phenotype [J].Molecular Ecology，2022，31（4）：1269-1280.

[32]TRAPNELL C，WILLIAMS B A，PERTEA G，et al.Transcript assembly and quantification by RNA-seq reveals unannotated transcripts and isoform switching during cell differentiation [J].Nature Biotechnology，2010，28（5）：511-515.

[33]ANDERS S，HUBER W.Differential expression analysis for sequence count data [J].Genome Biology，2010，" 11（10）：R106.

[34]ERNST J，BAR-JOSEPH Z.STEM：a tool for the analysis of short time series gene expression data [J].BMC Bioinformatics，2006，7：191.

[35]KARP X.Hormonal regulation of diapause and development in nematodes，insects，and fishes [J].Frontiers in Ecology and Evolution，2021（9）：735924.

[36]易傳輝，陳曉鳴，史軍義，等.光周期對柑橘鳳蝶幼蟲生長發育及存活的影響[J].西北農業學報，2009，18（1）：162-164，169.

[37]TAKEDA M，SUZUKI T.Circadian and neuroendocrine basis of photoperiodism controlling diapause in insects and mites：a review [J].Frontiers in Physiology，2022，13：867621.

Transcriptome" Sequencing Analysis during Pupal Diapause in ""Bombyx mandarina and Screening" of Genes Associated to Pupal Diapause

MENG Gang，CHU Qu，WANG Ruixian，LIU Qiang，PENG Yunwu， CHEN Anli and" LING Jun

（Shaanxi Key Laboratory of Sericulture，Ankang University，Ankang" Shaanxi 725000，China）

Abstract This study aims to screen the regulatory genes associated with pupal diapause in wild mulberry silkworms，providing a theoretical basis for revealing the regulatory mechanism of pupal diapause at the molecular level. Pupae at different developmen stages were utilized as materials for transcriptome sequencing using the Illumina HiSeqTM 2000 platform.KEGG enrichment analysis for differentially expressed genes was conducted by the KAAS online toolkit with the conditions of padjlt;0.05 and fold change ≥32 or padjlt;0.05 and fold change≤1/32. Temporal series analysis was conducted on total DEGs using stem toolkit. The genes with low expression in the early and late stages of diapause，high expression in the diapause stage were selected as the pupal diapause associated genes，then GO and KEGG enrichment analysis were performed. A total of 13 764，13 765，13 765，and 13 764 transcripts were detected in the pre-，mid-，late-，and post-diapause stages，respectively.Pairwise comparisons revealed 4 621 DEGs，with 1 921，451 and 38" genes upregulated，1 143，171 and 59 genes downregulated" in middle-，late-，and release pupa diapause stages，respectively，compared to pre-diapause.Temporal series analysis identified 290 pupal diapause associated genes，with GO enrichment"" analysis mapping 202 genes to 52 secondary GO categories and KEGG enrichment analysis mapping 67 genes to 126 tertiary metabolic pathways. The study provides a theoretical reference for screening key genes involved in pupal diapause regulation in wild mulberry silkworms.

Key words Bombyx mandarina；Transcriptome；Temporal-series analysis；Pupal diapause-associated genes

Received 2023-12-01

Returned 2024-03-12

Foundation item Key Ramp;D" Plan of Shaanxi Province （No. 2020NY-014，No.2023-JC-YB-206）; Open Foundation of State Key Laboratory of Silkworm Genome Biology （No. SKLSQB1819-4）; Education Department of Shaanxi Province in China （No.20JS002）; the Sci-tech Project of Ankang City （No. AK2022-NY-06）.

First author MENG Gang，male，associate research fellow. Research area：silkworm germplasm resources and genetic breeding. E-mail：nsymg@aku.edu.cn

（責任編輯：顧玉蘭 Responsible editor：GU Yulan）

收稿日期：2023-12-01

修回日期：2024-03-12

基金項目：陜西省重點研發計劃（2020NY-014，2023-JC-YB-206）；家蠶基因組生物學國家重點實驗室開放課題（SKLSQB1819-4）；陜西省教育廳重點科學研究計劃（20JS002）；安康市科技計劃項目（AK2022-NY-06）。

第一作者：孟 剛，男，副研究員，研究方向為蠶種質資源與遺傳育種。E-mail：nsymg@aku.edu.cn