二化螟幼蟲熱脅迫響應(yīng)的轉(zhuǎn)錄組分析

朱 宇，劉 洋

(溫州科技職業(yè)學(xué)院，浙江 溫州 325006)

溫度是影響昆蟲生長、發(fā)育、繁殖、存活、分布范圍和多樣性的重要非生物因子[1-2]。昆蟲是變溫動物，對溫度的變化敏感[3]。在高溫條件下，昆蟲發(fā)育速度變得更快，溫度升高會使昆蟲向高緯度地區(qū)擴(kuò)張遷移[1]。高溫會導(dǎo)致昆蟲種群動態(tài)的改變、物種豐富度的下降，以及昆蟲存活率的降低[3-5]。昆蟲也可以通過調(diào)整自身的生理活動來應(yīng)對高溫環(huán)境，包括熱激蛋白的表達(dá)、丙醛濃度的增高以及過氧化氫酶、谷胱甘肽S-轉(zhuǎn)移酶和超氧化物歧化酶等抗氧化物質(zhì)的提高[6-8]。隨著全球氣候變暖，越來越多的研究者開始關(guān)注昆蟲對高溫或熱脅迫的適應(yīng)機(jī)制[9]。

二化螟是水稻和茭白等作物上的重要害蟲，每年給農(nóng)業(yè)生產(chǎn)造成巨大的經(jīng)濟(jì)損失[2]。二化螟對溫度的適應(yīng)性強(qiáng)，其越冬幼蟲能在-20 ℃下存活[10]。羅舉等[11]發(fā)現(xiàn)，每天在40 ℃下4 h的二化螟幼蟲，也能發(fā)育為成蟲。二化螟發(fā)育的最適溫度在27 ℃左右，發(fā)育溫度的最低閾值在10～15 ℃，最高閾值在35～40 ℃，在40 ℃溫度下孵出的幼蟲不能存活[11-12]。在高溫脅迫下，二化螟幼蟲體內(nèi)會產(chǎn)生一系列的生理反應(yīng)來消除氧化應(yīng)激反應(yīng)對身體的傷害，包括熱激蛋白的產(chǎn)生、抗氧化酶的增多等[13-18]。目前，有關(guān)二化螟高溫脅迫響應(yīng)的研究主要集中在單個(gè)基因的功能，尚未通過轉(zhuǎn)錄組學(xué)的方法對二化螟全基因表達(dá)的影響進(jìn)行分析。

隨著第二代測序技術(shù)的成本降低，轉(zhuǎn)錄組高通量測序方法已經(jīng)被應(yīng)用到多種生物學(xué)試驗(yàn)中[19]。通過高通量測序，我們可以同時(shí)測定幾萬個(gè)基因的表達(dá)水平，并在全轉(zhuǎn)錄組范圍內(nèi)找出不同試驗(yàn)處理下的差異表達(dá)基因[20]。前人已經(jīng)有許多關(guān)于昆蟲熱脅迫響應(yīng)的比較轉(zhuǎn)錄組學(xué)研究[21-23]。在這些研究中，亞洲柑橘木虱(Diaphorinacitri)在熱脅迫下，它的分子伴侶、蛋白生物合成和抗氧化相關(guān)基因被誘導(dǎo)表達(dá)[21]。高溫也可以誘導(dǎo)桑絹絲野螟(GlyphodespyloalisWalker)熱激蛋白、解毒酶和氧化還原酶的表達(dá)，但糖酵解、糖異生和脂肪酸的生物合成等代謝活動被抑制[23]。東方粘蟲(MythimnaseparataWalker)在30 ℃熱脅迫下，其熱激蛋白和泛素介導(dǎo)的蛋白酶體的表達(dá)會顯著上調(diào)[22]。本研究通過比較轉(zhuǎn)錄組學(xué)的方法，分析二化螟幼蟲應(yīng)對熱脅迫的差異表達(dá)基因及相關(guān)代謝通路，從而更深入地了解二化螟幼蟲對高溫脅迫的適應(yīng)機(jī)理。

1 材料與方法

1.1 供試?yán)ハx與樣品準(zhǔn)備

二化螟幼蟲從浙江省溫州市種子種苗基地水稻田采集，并在室內(nèi)用稻苗飼養(yǎng)法進(jìn)行飼養(yǎng)，溫度設(shè)為(27±1)℃，光照黑暗時(shí)間比16 h/8 h[24]。試驗(yàn)時(shí)選取大小基本一致的二化螟四齡幼蟲。將熱脅迫組(heat stress，HS)和對照組(Control)二化螟幼蟲分別置于38 ℃和26 ℃中2 h，隨后，從對照組和熱脅迫組中分別取2頭二化螟放入離心管中，每組處理做3次生物學(xué)重復(fù)，加入TRIzol以備提取RNA。

1.2 cDNA文庫構(gòu)建與測序

RNA的提取、RNA質(zhì)量檢測、cDNA文庫的構(gòu)建及上機(jī)測序由武漢未來組生物科技有限公司(武漢)完成，測序平臺為Illumina NovaSeq 6000。

1.3 轉(zhuǎn)錄組拼裝與注釋

由Illumina平臺測序得到的原始圖像數(shù)據(jù)經(jīng)過堿基識別轉(zhuǎn)化為原始數(shù)據(jù)(raw data)，使用Trimmomatic軟件對原始數(shù)據(jù)進(jìn)行過濾得到clean data[25]。然后使用FastQC軟件對clean data進(jìn)行質(zhì)控[26]。使用HISAT2和StringTie軟件，將獲得的clean data與二化螟參考基因組(下載地址ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/004/000/445/GCA_004000445.1_ASM400044v1)進(jìn)行比對拼裝得到轉(zhuǎn)錄本，并計(jì)算出每個(gè)基因在不同處理中的FPKM值[27-28]。

1.4 差異基因表達(dá)分析及功能富集分析

通過DESeq2軟件對熱脅迫組和對照組進(jìn)行差異表達(dá)分析[29]。對差異檢驗(yàn)的P值作多重假設(shè)檢驗(yàn)校正，通過控制FDR(false discovery rate)來決定P值的閾值[30]。在本次分析中，差異表達(dá)基因的閾值為FDR≤0.05。利用NCBI的BLASTN軟件、GO(gene ontology)和KEGG(Kyoto encyclopedia of genes and genomes)數(shù)據(jù)庫對篩選出的差異表達(dá)基因進(jìn)行功能注釋[31-33]。利用顯著上調(diào)的差異表達(dá)基因和KEGG富集的結(jié)果，應(yīng)用超幾何檢驗(yàn)，找出顯著富集的KEGG 代謝途徑[32]。

2 結(jié)果與分析

2.1 二化螟轉(zhuǎn)錄組測序與序列組裝

我們應(yīng)用RNA-seq測序技術(shù)對兩組不同溫度下二化螟幼蟲的基因表達(dá)進(jìn)行了測定。轉(zhuǎn)錄組測序的原始數(shù)據(jù)已經(jīng)提交到NCBI，BioProject登錄號為PRJNA578414。通過測序，我們共得到34.2 GB 原始數(shù)據(jù)，425 986 598個(gè)clean reads(表1)。6個(gè)樣品(HS_1、HS_2、HS_3、Control_1、Control_2、Control_3) clean data的Q30值均大于93%(表1)。使用HISAT2軟件將6個(gè)樣品產(chǎn)生的測序序列與二化螟參考基因組進(jìn)行比對，能夠比對到參考基因組上的序列占總體的百分比均超過89%，其中，在參考基因組中有唯一比對位置的測序序列占總體的百分比均高于82%，有多個(gè)比對位置的序列占總體的百分比均不超過9%(表1)。

2.2 二化螟在熱脅迫下的差異表達(dá)基因

通過比較兩組不同溫度處理后的二化螟轉(zhuǎn)錄組，得到387個(gè)差異表達(dá)基因，其中221個(gè)基因?yàn)樯险{(diào)表達(dá)，166個(gè)基因?yàn)橄抡{(diào)表達(dá)(圖1)。通過BLASTN對這些差異表達(dá)基因進(jìn)行注釋，其中，有11個(gè)基因的注釋結(jié)果為熱激蛋白(heat shock protein)，它們在熱脅迫下上調(diào)表達(dá)。這些熱激蛋白中，有1個(gè)屬于熱激蛋白90家族，有4個(gè)屬于熱激蛋白70家族，2個(gè)屬于熱激蛋白40家族，還有4個(gè)是小分子熱激蛋白(圖2)。此外，在上調(diào)表達(dá)的差異表達(dá)基因中，有3個(gè)基因?yàn)榧?xì)胞色素P450(CYP，cytochrome P450)(圖3)。除了這兩類蛋白外，在上調(diào)表達(dá)的差異表達(dá)基因中，未找到其他與抗氧化相關(guān)的酶。

2.3 差異表達(dá)基因的GO注釋和KEGG通路分析

利用GO數(shù)據(jù)庫對所獲得的差異表達(dá)基因進(jìn)行功能注釋，并對差異表達(dá)基因進(jìn)行功能分類統(tǒng)計(jì)(圖4)。在上調(diào)的差異表達(dá)基因中，有136個(gè)基因比對到生物過程(biological process)數(shù)據(jù)庫，其中，代謝過程(metabolic process)這一條目中的差異基因數(shù)量最多，其次是單組織過程(single-organism process)和細(xì)胞過程(cellular process)。有90個(gè)上調(diào)差異表達(dá)基因比對到細(xì)胞組分(cellular component)數(shù)據(jù)庫中，其中，細(xì)胞膜(membrane)中的差異基因數(shù)量最多，其次是細(xì)胞(cell)和細(xì)胞部分(cell part)。有105個(gè)上調(diào)差異表達(dá)基因比對到分子功能(molecular function)數(shù)據(jù)庫中，其中，結(jié)合(binding)這一條目中的差異基因數(shù)量最多，其次是催化活性(catalytic activity)。

利用KEGG數(shù)據(jù)庫對差異表達(dá)基因進(jìn)行注釋，上調(diào)的221個(gè)差異表達(dá)基因分布在69條代謝通路中(圖5)。應(yīng)用超幾何檢驗(yàn)，找到上調(diào)表達(dá)基因顯著富集的代謝通路4條，分別為內(nèi)質(zhì)網(wǎng)內(nèi)蛋白質(zhì)加工、壽命調(diào)節(jié)通路-多物種、植物病原體相互作用和抗原處理及呈現(xiàn)(表2)。在內(nèi)質(zhì)網(wǎng)內(nèi)蛋白質(zhì)加工代謝通路的9個(gè)基因中， 8個(gè)為熱激蛋白，1個(gè)為熱激蛋白結(jié)合蛋白。

表1 二化螟轉(zhuǎn)錄組組裝結(jié)果統(tǒng)計(jì)

Table 1 Summary statistics of assembly results ofChilosuppressalistranscriptome

Q30比率表示質(zhì)量值大于等于30的堿基所占百分比。

Q30 rate， The percentage of bases with a quality score of 30 or higher.

紅點(diǎn)表示上調(diào)表達(dá)的差異表達(dá)基因；綠點(diǎn)表示下調(diào)表達(dá)的差異表達(dá)基因；藍(lán)點(diǎn)表達(dá)差異不顯著的基因；Fold change，差異倍數(shù)；padj，經(jīng)過FDR校正的P值。Red dot indicated up-regulated DEGs (differentially expressed genes); Green dot indicated down-regulated DEGs; Blue dot indicated non-DEGs; Fold change， The ratio of heat stress group versus control group; padj， FDR-adjusted P value.圖1 差異基因火山圖Fig.1 Volcano map of DEGs

HS，熱脅迫組；Control，對照組；GeneID，二化螟參考基因組注釋文件中的基因ID；HSP，熱激蛋白；FPKM，每1百萬個(gè)匹配上的讀長中匹配到外顯子的每1 000個(gè)堿基上的片段個(gè)數(shù)。HS， Heat stress group; Control， control group; GeneID， Gene ID from reference annotation file of Chilo suppressalis; HSP， Heat shock protein; FPKM， Fragments per kilobase of exon model per million reads mapped.圖2 二化螟幼蟲在熱脅迫下差異表達(dá)的熱激蛋白Fig.2 Differentially expressed genes encoding for heat shock proteins in Chilo suppressalis larvae under high temperature stress

HS，熱脅迫組；Control，對照組；GeneID，二化螟參考基因組注釋文件中的基因ID；CYP，細(xì)胞色素P450；FPKM，每1百萬個(gè)匹配上的讀長中匹配到外顯子的每一千個(gè)堿基上的片段個(gè)數(shù)。HS， Heat stress group; Control， control group; GeneID， Gene ID from reference annotation file of Chilo suppressalis; CYP， Cytochrome P450; FPKM， Fragments per kilobase of exon model per million reads mapped.圖3 二化螟幼蟲在熱脅迫下差異表達(dá)的細(xì)胞色素P450Fig.3 Differentially expressed genes encoding for CYPs in Chilo suppressalis larvae under high temperature stress

Up，上調(diào)表達(dá)的差異表達(dá)基因；Down，下調(diào)表達(dá)的差異表達(dá)基因。Up， Up-regulated DEGs; Down， Down-regulated DEGs.圖4 差異表達(dá)基因GO分類Fig.4 GO classification of DEGs

3 討論

本研究通過高通量測序技術(shù)對分別置于38 ℃和26 ℃中2 h的二化螟幼蟲進(jìn)行轉(zhuǎn)錄組測序和拼裝，并進(jìn)行了差異表達(dá)分析。測序結(jié)果顯示，6個(gè)試驗(yàn)樣品的Q30值均大于93%，說明轉(zhuǎn)錄組測序的質(zhì)量較高。將轉(zhuǎn)錄組數(shù)據(jù)與二化螟參考基因組進(jìn)行比對拼接，比對到參考基因組上的序列占總體的百分比均超過89%，在參考基因組中有唯一比對位置的測序序列占總體的百分比均高于82%，有多個(gè)比對位置的序列占總體的百分比均不超過9%，說明轉(zhuǎn)錄組拼接的質(zhì)量較高。通過比較兩組不同溫度處理后的二化螟轉(zhuǎn)錄組，我們得到了387個(gè)差異表達(dá)基因(221個(gè)上調(diào)表達(dá)、166個(gè)下調(diào)表達(dá))，部分差異表達(dá)基因與二化螟幼蟲的解毒、抗氧化應(yīng)激、免疫以及代謝相關(guān)。通過對上調(diào)表達(dá)的差異表達(dá)基因進(jìn)行代謝通路富集分析，發(fā)現(xiàn)內(nèi)質(zhì)網(wǎng)內(nèi)蛋白質(zhì)加工、壽命調(diào)節(jié)通路-多物種、植物病原體相互作用和抗原處理及呈現(xiàn)這4條通路在二化螟熱脅迫時(shí)顯著富集，說明這些代謝通路在二化螟熱脅迫響應(yīng)的過程中發(fā)揮著重要的作用。

Rich factor表示富集上的基因占注釋到的基因的百分比；縱坐標(biāo)表示富集上的條目；Gene number表示富集上的基因個(gè)數(shù)；顏色表示的Q value值，越低越顯著。Rich factor indicated the percentage of the gene on the enrichment to the annotated gene; the ordinate indicated the term of the enrichment; Gene number indicated the number of genes on the enrichment， and the Q value represented by the color， the lower the more significant.圖5 上調(diào)差異表達(dá)基因前20條富集的KEGG代謝通路Fig.5 The top 20 enriched KEGG pathways of the up-regulated differentially expressed genes

表2 二化螟上調(diào)表達(dá)差異基因的顯著富集代謝通路

Table 2 Summary of significantly enriched pathway terms associated with up-regulated differentially expressed genes

通路Pathway基因數(shù)Genenumber參考基因數(shù)Reference genenumber校正后的P值Corrected P valueMap登錄號Map ID內(nèi)質(zhì)網(wǎng)內(nèi)蛋白質(zhì)加工Protein processing in endoplasmic reticulum91300.000416map04141壽命調(diào)節(jié)通路-多物種Longevity regulating pathway-multiple species5580.008634map04213植物病原體相互作用Plant-pathogen interaction3140.008634map04626抗原處理及呈現(xiàn)Antigen processing and presentation3220.025823map04612

在昆蟲體內(nèi)，熱激蛋白通常和分子伴侶或輔助蛋白一起協(xié)作，通過促進(jìn)變形蛋白修復(fù)或參與蛋白的合成、折疊和組裝等生命活動來保持昆蟲細(xì)胞內(nèi)環(huán)境的穩(wěn)定[34]。熱激蛋白可以在昆蟲的正常生長發(fā)育過程中表達(dá)，也會在高溫、寒冷、種群密度過高、缺氧以及農(nóng)藥等脅迫環(huán)境下誘導(dǎo)表達(dá)[35-39]。根據(jù)相對分子質(zhì)量大小，熱激蛋白被劃分為許多家族，包括小分子熱激蛋白、熱激蛋白40、熱激蛋白60、熱激蛋白70和熱激蛋白90[34]。熱激蛋白能在高溫脅迫下誘導(dǎo)表達(dá)，已經(jīng)被許多文獻(xiàn)報(bào)道。在熱脅迫下，蔥蠅(Deliaantiqua)非滯育蛹的熱激蛋白90表達(dá)量增高[40]。肉蠅(Sarcophagacrassipalpis)非滯育蛹中的熱激蛋白70在高溫脅迫下上調(diào)表達(dá)[41]。在麥紅吸漿蟲(Sitodiplosismosellana)中也有類似結(jié)果，麥紅吸漿蟲滯育幼蟲的熱激蛋白70和熱激蛋白90分別在35 ℃和30 ℃表達(dá)量最高[42]。本研究中，1個(gè)二化螟熱激蛋白90和4個(gè)熱激蛋白70在高溫脅迫下上調(diào)表達(dá)。小分子熱激蛋白是一種十分普遍的蛋白，它能夠阻止蛋白聚集，幫助細(xì)胞適應(yīng)外界不良環(huán)境[43]。本研究中，4個(gè)二化螟小分子熱激蛋白在高溫條件下大量表達(dá)。上述結(jié)果說明，二化螟熱激蛋白70、熱激蛋白90和小分子熱激蛋白能夠提高二化螟幼蟲對高溫環(huán)境的適應(yīng)。

細(xì)胞色素P450是一類包含血紅蛋白的家族蛋白[44]。昆蟲的細(xì)胞色素P450基因家族有4個(gè)進(jìn)化枝，分別為CYP2家族、CYP3家族、CYP4家族和線粒體P450家族[45]。昆蟲細(xì)胞色素P450基因能夠參與昆蟲的多種生理活動，包括昆蟲的激素代謝、植物化感物質(zhì)的解毒、昆蟲的抗藥性以及昆蟲的抗氧化應(yīng)激反應(yīng)等[46-48]。細(xì)胞色素P450參與昆蟲的抗氧化應(yīng)激作用，已經(jīng)被許多文獻(xiàn)報(bào)道。在埃及伊蚊(Aedesaegypti)體內(nèi)發(fā)現(xiàn)了4個(gè)細(xì)胞色素P450基因，參與抗氧化應(yīng)激作用[48]。在中華蜜蜂(Apisceranacerana)體內(nèi)，也發(fā)現(xiàn)1個(gè)細(xì)胞色素P450基因AccCYP336A1，具有抗氧化應(yīng)激作用[49]。昆蟲在高溫環(huán)境下，會產(chǎn)生活性氧和其他有毒物質(zhì)，對昆蟲細(xì)胞造成損傷[50]。在本研究中，有3個(gè)細(xì)胞色素P450基因在38 ℃處理2 h后顯著上調(diào)表達(dá)，說明細(xì)胞色素P450可能在二化螟熱脅迫響應(yīng)和抗氧化應(yīng)激過程中發(fā)揮作用。

本研究利用轉(zhuǎn)錄組測序技術(shù)對分別置于38 ℃和26 ℃中2 h的二化螟幼蟲進(jìn)行轉(zhuǎn)錄組測序和分析，從轉(zhuǎn)錄組水平揭示短時(shí)高溫脅迫對二化螟幼蟲的影響。結(jié)果表明，二化螟幼蟲在高溫脅迫過程中，熱激蛋白和細(xì)胞色素P450這兩類與熱脅迫響應(yīng)相關(guān)的基因表達(dá)量發(fā)生了顯著變化。這些結(jié)果可為二化螟熱脅迫響應(yīng)相關(guān)的試驗(yàn)提供參考。