張 博, 韓海斌, 徐林波, 高書晶, 甘 霖, 岳方正, 劉愛萍,*

(1.中國農(nóng)業(yè)科學院草原研究所, 呼和浩特 100010; 2.阿拉善盟科技信息研究所, 內(nèi)蒙古阿拉善 750399;3.國家林業(yè)和草原局森林和草原病蟲害防治總站, 沈陽 110034)

傘裙追寄蠅Exoristacivilis是一種牧草優(yōu)勢寄生性天敵，隸屬于雙翅目(Diptera)寄蠅科(Tachinidae)，能夠?qū)φ诚xMythimnaseparata、草地螟Loxostegesticticalis、草地貪夜蛾Spodopterafrugiperda等多種牧草害蟲起到有效的防控作用(王建梅等, 2013, 2015)。昆蟲滯育是一種生理過程，其特征在于行為修飾、代謝抑制以及增強的脅迫耐受性(King and Macrae, 2015)。近年來，關(guān)于傘裙追寄蠅滯育類型、滯育溫光周期的誘導(dǎo)以及滯育關(guān)聯(lián)蛋白的研究取得了一定的成果(韓海斌等, 2018; 姜珊, 2018)，為揭示其滯育的分子機制奠定了基礎(chǔ)。

新一代測序技術(shù)的出現(xiàn)，極大促進了昆蟲轉(zhuǎn)錄組學研究(張棋麟和袁明龍, 2013)。目前利用高通量技術(shù)對雙翅目昆蟲的滯育研究也取得了一定的進展，利用高通量技術(shù)發(fā)現(xiàn)了在柑橘大實蠅Bactroceraminax滯育調(diào)控的相關(guān)研究中篩選出12個顯著上調(diào)和6個顯著下調(diào)的滯育關(guān)聯(lián)基因以及重要的滯育關(guān)聯(lián)蛋白hsp23, hsp70和hsp90(Luetal., 2014; Wangetal., 2019)。張倩等(2019)利用iTRAQ技術(shù)篩選淡色庫蚊Culexpipienspallens越冬后期和滯育期的差異表達蛋白，共鑒定出244個差異表達蛋白，包括126個上調(diào)蛋白和118個下調(diào)蛋白，這些差異表達蛋白的功能與脂質(zhì)代謝、細胞骨架重塑、碳水化合物代謝等有關(guān)。

目前，調(diào)控傘裙追寄蠅滯育的分子調(diào)控機制仍然不完全明確。為此，本研究對傘裙追寄蠅正常發(fā)育及滯育的蛹進行轉(zhuǎn)錄組測序，對篩選出的滯育關(guān)聯(lián)基因進行KEGG通路富集分析，了解滯育關(guān)聯(lián)基因主要參與的代謝途徑和信號轉(zhuǎn)導(dǎo)機制，為進一步研究傘裙追寄蠅的滯育調(diào)控基因以及分子機理奠定基礎(chǔ)，并為通過遺傳信息研發(fā)滯育劑、滯育解除劑，從而實現(xiàn)人工調(diào)控傘裙追寄蠅等天敵昆蟲的壽命，實現(xiàn)在田間精準釋放，發(fā)揮生物防治的最佳效果提供科學依據(jù)。

1 材料與方法

1.1 供試昆蟲及滯育誘導(dǎo)

傘裙追寄蠅蛹為中國農(nóng)業(yè)科學院草原研究所實驗室所飼養(yǎng)。將正常發(fā)育的非滯育蛹在25℃±1℃室溫條件下培養(yǎng)1 d，經(jīng)用液氮處理后放入-80℃冰箱保存?zhèn)溆?；將正常發(fā)育的非滯育蛹置于人工氣候箱內(nèi)，在溫度11±1℃、相對濕度70%±10%、光周期8L∶16D的條件下，經(jīng)過40 d誘導(dǎo)可獲得滯育態(tài)蛹(姜珊, 2018)，40 d后選取重量與正常發(fā)育非滯育蛹一致的存活滯育蛹，經(jīng)用液氮處理后放入-80℃冰箱保存?zhèn)溆?。然后，將備用的傘裙追寄蠅正常發(fā)育非滯育蛹和滯育蛹進行RNA樣品檢測，主要方法為：(1)瓊脂糖凝膠電泳分析RNA降解程度以及是否有污染；(2)NanoDrop檢測RNA的純度(OD260/OD280比值)；(3)Qubit對RNA濃度進行精確定量；(4)Agilent 2100精確檢測RNA的完整性。

1.2 轉(zhuǎn)錄組測序和組裝

傘裙追寄蠅非滯育和滯育蛹的轉(zhuǎn)錄組測序工作由北京諾禾致源科技股份有限公司完成。測序進行3個生物學重復(fù)，每個生物學重復(fù)以4頭蛹為樣本，兩種蛹共構(gòu)建6個文庫。樣品檢測合格后，用帶有Oligo(dT)的磁珠富集真核生物mRNA。隨后加入fragmentation buffer將mRNA打斷成短片段，以mRNA為模板，用六堿基隨機引物(random hexamers)合成cDNA第1鏈，然后加入緩沖液，再加純化的雙鏈cDNA。純化的雙鏈cDNA先進行末端修復(fù)、加A尾并連接測序接頭，再進行片段大小選擇。最后進行PCR擴增，得到最終的文庫。文庫構(gòu)建完成后，先使用Qubit2.0進行初步定量，稀釋文庫至1.5 ng/μL，隨后使用Agilent 2100對文庫的insert size進行檢測，insert size符合預(yù)期后，使用Q-PCR方法對文庫的有效濃度進行準確定量(文庫有效濃度>2 nmol/L)，以保證文庫質(zhì)量。庫檢合格后，把不同文庫按照有效濃度及目標下機數(shù)據(jù)量的需求pooling后進行高通量測序。測序完成后，高通量測序(Illumina HiSeqTM2000)得到的原始圖像數(shù)據(jù)文件經(jīng)堿基識別(base calling)分析轉(zhuǎn)化為原始測序序列，然后進行測序錯誤率分布檢查，獲得高質(zhì)量的序列(clean reads)后，采用Trinity軟件對clean reads進行拼接，Corset參差聚類后得到最長的序列進行后續(xù)法分析。

1.3 生物信息學分析

為獲得全面的基因功能信息，對傘裙追寄蠅進行了七大數(shù)據(jù)庫的基因功能注釋，包括：NCBI非冗余蛋白序列數(shù)據(jù)庫(NCBI non-redundant protein sequences, Nr)，NCBI非冗余核苷酸序列數(shù)據(jù)庫(NCBI non-redundant nucleotide sequences, Nt)，蛋白家族(protein family, Pfam)，直系同源基因簇(cluster of orthologous groups of proteins, KOG/COG), Swiss-Prot數(shù)據(jù)庫(a manually annotated and reviewed protein sequence database), 京都基因與基因組百科全書(KEGG ortholog database, KO)，基因本體論(gene ontolooy, GO)。拼接得到的傘裙追寄蠅轉(zhuǎn)錄組測序結(jié)果與UniProt 或者Nr數(shù)據(jù)庫進行Blastx比對(將核苷酸序列翻譯為蛋白，再進行比對)，篩選條件E-value<1e-5，得到注釋基因。通過RSEM軟件得到每個樣品比對到每個基因上的readcount數(shù)目，并對其進行FPKM(expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced)轉(zhuǎn)換，分析基因的表達水平；對轉(zhuǎn)錄組測序的非滯育蛹和滯育蛹的樣品進行了樣品間的相關(guān)性檢查，采用R語言進行Pearson相關(guān)系數(shù)的計算；將正常非滯育蛹和滯育蛹的轉(zhuǎn)錄組測序結(jié)果進行比較，采用校正后的P值(p-adjusted, padj)與差異表達倍數(shù)(fold change, FC)方法進行差異表達基因挑選，padj是對P值進行了多重假設(shè)檢驗，能更好地控制假陽性率。對篩選出的差異表達基因進行聚類分析，使用聚類軟件包pheatmap，針對差異表達基因的并集，以基因的相對表達水平值進行聚類。采用相應(yīng)的距離算法計算每個基因之間的距離，然后通過反復(fù)迭代，計算基因之間的相對距離，最后根據(jù)基因的相對距離遠近分成不同的subcluster，從而實現(xiàn)聚類。對傘裙追寄蠅非滯育蛹和滯育蛹差異表達基因的篩選條件為: padj<0.05且fold change≥32 或 padj<0.05 且fold change≤1/32。與非滯育蛹的相比，滯育蛹基因差異表達倍數(shù)在32倍以上的滯育關(guān)聯(lián)基因，并進行重點分析。

1.4 滯育關(guān)聯(lián)基因 KEGG通路分析

將滯育關(guān)聯(lián)基因序列制成FASTA格式文件，應(yīng)用KAAS(KEGG Automatic Annotation Server.http:∥www.genome.jp/kegg/)在線pathway比對分析工具進行KEGG通路富集分析。

2 結(jié)果

2.1 傘裙追寄蠅非滯育和滯育蛹轉(zhuǎn)錄組測序結(jié)果

本研究對傘裙追寄蠅非滯育和滯育蛹進行轉(zhuǎn)錄組測序，高質(zhì)量樣品數(shù)據(jù)超過7.42 Gb，每個樣品的Q20均大于96.08%，Q30均大于90.01%(表1)，通過轉(zhuǎn)錄組數(shù)據(jù)結(jié)果發(fā)現(xiàn)其準確度較高，可用于后續(xù)分析。

表1 傘裙追寄蠅非滯育和滯育蛹轉(zhuǎn)錄組測序數(shù)據(jù)質(zhì)量評估統(tǒng)計

樣品間基因表達水平相關(guān)性是檢驗實驗可靠性和樣本選擇是否合理的重要指標。樣品間基因表達水平的相關(guān)性反映了樣品間的相似程度，即不同處理或組織的樣品在基因表達水平方面的相似度。相關(guān)系數(shù)越接近1，表明樣品之間基因表達模式的相似度越高，樣品間的差異表達基因也越少(張宇鐳等, 2005)。本研究中生物學重復(fù)間的樣品的相關(guān)系數(shù)大于生物學重復(fù)外的樣品的相關(guān)系數(shù)(圖1)，表明滯育蛹樣品與非滯育蛹樣品基因表達模式存在明顯差異。

圖1 傘裙追寄蠅非滯育和滯育蛹轉(zhuǎn)錄組樣品間基因表達水平Pearson相關(guān)系數(shù)熱圖

差異表達基因聚類分析用于判斷不同實驗條件下差異基因表達量的聚類模式(楊春梅等, 2003)。圖2紅色表示高表達，藍色表示低表達。本研究檢測結(jié)果中藍色和紅色部分能夠明顯區(qū)分，說明樣品整體的測序質(zhì)量好。

圖2 傘裙追寄蠅非滯育和滯育蛹轉(zhuǎn)錄組中差異表達基因聚類圖

2.2 滯育關(guān)聯(lián)基因 KEGG通路分析

傘裙追寄蠅非滯育蛹和滯育蛹轉(zhuǎn)錄組測序在KEGG直系同源系統(tǒng)(KEGG orthdogy, KO)中成功注釋了6 643個基因，KO系統(tǒng)將各個KEGG注釋系統(tǒng)聯(lián)系在一起，可完成轉(zhuǎn)錄組的功能注釋。根據(jù)篩選條件：padj<0.05且fold change≥2或padj<0.05且fold change≤0.5，獲得在傘裙追寄蠅非滯育蛹和滯育蛹轉(zhuǎn)錄組間差異2倍以上的差異表達基因共有2 648個，滯育蛹中上調(diào)表達基因1 518個，下調(diào)表達基因1 130個。對這些差異表達基因進行KEGG分類(圖3)，所屬KEGG代謝通路分為5類，包括細胞過程、環(huán)境信息處理、遺傳信息處理、代謝(metabolism)和有機系統(tǒng)，這些基因主要集中在信號轉(zhuǎn)導(dǎo)、內(nèi)分泌系統(tǒng)、細胞群落和碳水化合物代謝等途徑。與傘裙追寄蠅非滯育蛹相比，滯育蛹轉(zhuǎn)錄組中差異在32倍以上的滯育關(guān)聯(lián)基因為454個，可分為代謝、環(huán)境信息處理、遺傳信息處理、細胞過程和有機系統(tǒng)五大類，每類中滯育關(guān)聯(lián)基因參與的通路占比前3的如表2所示。其中上調(diào)基因為406個，下調(diào)基因為48個。通過KEGG數(shù)據(jù)庫分析，發(fā)現(xiàn)454個滯育關(guān)聯(lián)基因共參與到141個通路，其中參與焦點粘連通路的滯育關(guān)聯(lián)基因最多，為13個；其次是參與碳代謝通路的，為12個。

表2 傘裙追寄蠅滯育蛹轉(zhuǎn)錄組中滯育關(guān)聯(lián)基因KEGG通路分析

圖3 傘裙追寄蠅滯育蛹轉(zhuǎn)錄組中與非滯育蛹相比差異表達基因的KEGG功能分類

2.2.1代謝:根據(jù)KO分類結(jié)果，滯育關(guān)聯(lián)基因共參與代謝分類中的碳水化合物代謝(carbohydrate metabolism)、能量代謝(energy metabolism)、脂代謝(lipid metabolism)、核苷酸代謝(nucleotide metabolism)、氨基酸代謝(amino acid metabolism)等11個途徑，涉及47個通路。

分析結(jié)果顯示，在糖酵解/糖異生(glycolysis/gluconeogenesis)通路中，丙酮酸激酶(pyruvate kinase, PK)基因表達上調(diào)101.60倍，在檸檬酸循環(huán)(citrate cycle, TCA cycle)中，丙酮酸脫氫酶的E1組分(pyruvate dehydrogenase E1 component, PDHA)α亞基基因表達上調(diào)151.49倍，PK是糖酵解過程中發(fā)生第二次底物磷酸化的關(guān)鍵酶，催化反應(yīng)生成丙酮酸。在有氧條件下，丙酮酸脫氫酶復(fù)合體(pyruvate dehydrogenase complex, PDHC)催化由糖酵解產(chǎn)生的丙酮酸進行氧化脫羧反應(yīng)，促進乙酰輔酶A的生成，進入檸檬酸循環(huán)。在脂代謝中，脂酰輔酶A還原酶(fatty acyl-CoA reductase, FAR)基因表達上調(diào)，上調(diào)倍數(shù)最高可達146.49倍。在能量代謝中，參與氧化磷酸化的滯育關(guān)聯(lián)基因有9個，氧化磷酸化是一個電子傳遞過程，對線粒體呼吸電子鏈的傳遞十分重要。細胞色素C氧化酶是線粒體中電子傳遞鏈的末端酶，催化電子從細胞色素C轉(zhuǎn)移到氧。細胞色素氧化酶的5個亞基基因表達上調(diào)，其中細胞色素C氧化酶亞基4(cytochrome c oxidase subunit 4, Cox4)是1個關(guān)鍵亞基。在昆蟲激素的生物合成(insect hormone biosynthesis)通路中，保幼激素酯酶(juvenile-hormone esterase, JHE)基因表達下調(diào)，下調(diào)了82.28倍。

2.2.2有機系統(tǒng)：根據(jù)KO注釋結(jié)果，有機系統(tǒng)中的滯育關(guān)聯(lián)基因參與免疫系統(tǒng)、內(nèi)分泌系統(tǒng)、循環(huán)系統(tǒng)、排泄系統(tǒng)、神經(jīng)系統(tǒng)、感覺系統(tǒng)以及消化系統(tǒng)等7個系統(tǒng)，其中有58個滯育關(guān)聯(lián)基因參與內(nèi)分泌系統(tǒng)，涉及內(nèi)分泌系統(tǒng)的15個通路。

鈣調(diào)蛋白(calmodulin, CALM)基因上調(diào)表達主要涉及內(nèi)分泌系統(tǒng)，鈣調(diào)蛋白本身無任何酶活性，可以通過結(jié)合Ca2+形成Ca2+-CaM復(fù)合體，調(diào)控細胞代謝過程以及下游的靶蛋白、靶酶。Hsp70的一種：heat shock 70kDa protein 1/8(HSPA1s)在內(nèi)分泌系統(tǒng)的雌激素信號通路(estrogen signaling pathway)、免疫系統(tǒng)的抗原處理和展示(antigen processing and presentation)通路中既有上調(diào)也有下調(diào)，上調(diào)倍數(shù)最高可達98.53倍，下調(diào)倍數(shù)為2.56倍。

2.2.3細胞過程：根據(jù)KO注釋結(jié)果，參與細胞過程的73個滯育關(guān)聯(lián)基因主要涉及焦點粘連(focal adhesion)、肌動蛋白細胞骨架的調(diào)節(jié)(regulation of actin cytoskeleton)、緊密連接(tight junction)、縫隙連接(adherens junction)、吞噬體(phagosome)、溶酶體(lysosome)、細胞凋亡(apoptosis)、卵母細胞減數(shù)分裂(oocyte meiosis)、間隙連接(gap junction)、過氧化物酶體(peroxisome)、p53信號通路(p53 signaling pathway)、細胞周期(cell cycle)、胞吞作用(endocytosis)這13個通路。

在吞噬體通路中肌動蛋白基因表達占比最高，共有9個基因被注釋為肌動蛋白β/γ1(actin beta/gamma 1, ACTB_G1)，肌動蛋白是真核生物中最豐富的蛋白質(zhì)之一，對維持細胞形態(tài)具有重要的作用。在溶酶體中組織蛋白酶表達占比最高，共有6個基因分別被注釋為組織蛋白酶G(cathepsin G, CTSG)、組織蛋白酶S(cathepsin S, CTSS)、組織蛋白酶K(cathepsin K, CTSK)。

在細胞凋亡通路中，細胞色素C(cytochrome c, CYC)基因表達上調(diào)，CYC介導(dǎo)呼吸鏈中的電子傳遞，影響ATP的產(chǎn)生，其釋放可誘導(dǎo)細胞凋亡。在過氧化氫物酶體通路中，過氧化氫酶(catalase, CAT)基因表達上調(diào)，上調(diào)123.17倍。

2.2.4環(huán)境信息處理：根據(jù)滯育關(guān)聯(lián)基因的KO比對結(jié)果，環(huán)境信息處理共涉及到信號轉(zhuǎn)導(dǎo)(signal transduction)、信號分子及其相互作用(signaling molecules and interaction)、膜轉(zhuǎn)運(membrane transport)3個途徑。60個滯育關(guān)聯(lián)基因共參與環(huán)磷酸腺苷(cAMP signaling pathway)，rap1信號通路(rap1 signaling pathway)，MAPK信號通路(MAPK signaling pathway)，鈣信號通路(calcium signaling pathway)，cGMP-PKG信號通路(cGMP-PKG signaling pathway)，hippo信號通路(hippo signaling pathway)等在內(nèi)的19個通路。在ECM-受體相互作用通路中，肌腱蛋白(tenascin, TN)基因表達下調(diào)，下調(diào)倍數(shù)最高可達51.98倍。

2.2.5遺傳信息處理：根據(jù)KO注釋結(jié)果，在遺傳信息處理第二層次的分類中涉及復(fù)制與修復(fù)(replication and repair)，折疊、分類和降解(folding, sorting and degradation)，翻譯(translation)和轉(zhuǎn)錄(transcription)4個途徑，第3層次分類中涉及內(nèi)質(zhì)網(wǎng)蛋白加工(protein processing in endoplasmic reticulum)、RNA運輸(RNA transport)、泛素介導(dǎo)的蛋白質(zhì)水解(ubiquitin-mediated proteolysis)等7個通路。泛素蛋白連接酶E3 XIAP連鎖凋亡抑制蛋白(E3 ubiquitin-protein ligase XIAP, XIAP)基因，ring結(jié)構(gòu)域蛋白1(ring-box protein 1, RBX1)基因共同參與了泛素介導(dǎo)的蛋白質(zhì)水解這一通路，RBX1基因表達上調(diào)，上調(diào)293.29倍；XIAP基因表達下調(diào)，下調(diào)42.81倍。

3 討論

人工誘導(dǎo)傘裙追寄蠅蛹進入滯育，可以有效延長其貨架期。本研究以傘裙追寄蠅非滯育和滯育蛹轉(zhuǎn)錄組為研究對象，通過對比分析，篩選出454個滯育關(guān)聯(lián)基因。這些基因主要集中在信號轉(zhuǎn)導(dǎo)、內(nèi)分泌系統(tǒng)、碳水化合物代謝等途徑。

在傘裙追寄蠅滯育期間，丙酮酸激酶基因、丙酮酸脫氫酶的E1組分α亞基表達上調(diào)，有利于中間產(chǎn)物——乙酰輔酶A生成，然后進入TCA循環(huán)(Naitoetal., 1998)。細胞色素氧化酶相關(guān)基因的表達上調(diào)，利于呼吸電子鏈的傳遞，為昆蟲的滯育期提供能量。由此可以推測糖酵解、檸檬酸循環(huán)、氧化磷酸化在傘裙追寄蠅滯育期間的能量供給并不是維持在一個低水平。在脂代謝中，F(xiàn)AR催化脂酰輔酶A還原成脂肪醇，這是生成脂肪醇的主要途徑。關(guān)于FAR的研究主要集中在鱗翅目昆蟲信息素，在家蠶Bombyxmori和巢蛾屬Yponomeuta的研究中表明一種FAR催化性腺中多種性信息素，而歐洲玉米螟Ostrinianubilalis兩個亞種的研究則表明2個FAR基因生成不同的性信息素(Motoetal., 2003; Wahlenetal., 2009; Lassanceetal., 2010; 楊璞等, 2012)。JHE由脂肪體細胞和上皮細胞產(chǎn)生，是主要的保幼激素(juvenile hormone, JH)特異性降解酶之一，裂解JH酯鍵產(chǎn)生甲醇和JH酸，在調(diào)節(jié)JH代謝中發(fā)揮重要作用(李勝等, 2004)。JHE基因表達下調(diào)，由此可以推測滯育期間傘裙追寄蠅體內(nèi)保幼激素降解受到阻礙，JH含量可能不是保持在一個低水平。

肌腱蛋白屬于大細胞外基質(zhì)糖蛋白家族，在神經(jīng)系統(tǒng)發(fā)育過程中以特定模式表達，并在成熟后下調(diào)。大量的體外研究表明，肌腱蛋白在神經(jīng)系統(tǒng)中的作用與前體細胞遷移，軸突生長和引導(dǎo)有關(guān)(Joester and Faissner, 2001)?？梢酝茰y,在滯育期間傘裙追寄蠅蛹的肌腱蛋白可能處于發(fā)育成熟階段。

過氧化氫酶基因在過氧化氫物酶體通路發(fā)揮著保護酶的作用，使傘裙追寄蠅在滯育期間抵抗逆境脅迫，避免受到損傷(于德玲和王昌留, 2016)。在菜蛾盤絨繭蜂Cotesiavestalis滯育期間，發(fā)現(xiàn)低濃度的過氧化氫和高活性的CAT有利于增強該蜂的抗逆能力(郝仲萍等, 2012; 胡甘雨等, 2014)。

細胞凋亡是一種細胞程序性死亡過程，由許多過程相互關(guān)聯(lián)，共同調(diào)控。在傘裙追寄蠅滯育期間，XIAP基因、RBX1基因共同參與泛素介導(dǎo)的蛋白質(zhì)水解通路，RBX1是泛素連接酶E3的一個催化亞基，可以控制許多短壽命調(diào)節(jié)蛋白的更新(Jiaetal., 2011; 王步勇等, 2013)。在秀麗隱桿線蟲Caenorhabditiselegans的胚胎發(fā)育，生殖細胞豐度以及成年活力或壽命具有重要的調(diào)控作用(Moore and Boyd, 2004)。XIAP是具有抑制細胞凋亡的蛋白，RBX1與XIAP發(fā)揮了拮抗作用。RBX1具有泛素連接酶E3活性，可以引起XIAP自身泛素化而被降解，同時也能夠引起底物發(fā)生泛素化，抑制細胞凋亡的線粒體途徑，發(fā)揮抗凋亡的作用(李中海等, 2003)。

本研究通過對傘裙追寄蠅非滯育蛹與滯育蛹的轉(zhuǎn)錄組差異表達基因KEGG富集通路的分析比較，揭示了其滯育的關(guān)鍵基因，為闡明其滯育機理奠定了一定的理論基礎(chǔ)。關(guān)鍵滯育關(guān)聯(lián)基因在傘裙追寄蠅滯育過程中的差異性表達及作用，還有待于進一步驗證與探究。