PEG脅迫下玉米轉(zhuǎn)錄組差異性分析

124

(1.河南牧業(yè)經(jīng)濟(jì)學(xué)院食品與生物工程學(xué)院(酒業(yè)學(xué)院),鄭州 450046;2.河南牧業(yè)經(jīng)濟(jì)學(xué)院河南省白酒風(fēng)格工程技術(shù)研究中心,鄭州 450046;3.河南牧業(yè)經(jīng)濟(jì)學(xué)院鄭州市白酒釀造微生物技術(shù)重點(diǎn)實(shí)驗(yàn)室,鄭州 450046;4.張弓老酒酒業(yè)有限公司,河南 寧陵 476733)

玉米是世界上最重要的糧食作物之一，不僅可食用，而且還廣泛用于白酒釀造及酒精生產(chǎn)；但其又是旱地作物中對水需求量較多，并且對水分脅迫較為敏感的作物[1-3]。在中國，玉米有著悠久的種植和培育歷史，其種植面積約占耕地面積的25%，年產(chǎn)量可達(dá)1.2億t。而干旱是限制玉米生產(chǎn)發(fā)展的第一要素，導(dǎo)致70%以上的玉米遭受干旱威脅，每年造成的產(chǎn)量損失總計達(dá)到1 500萬t以上[4-7]，因此研究玉米的抗旱機(jī)理尤為重要。

劉成等研究發(fā)現(xiàn)，在干旱條件下，各種酶活性增加并促進(jìn)體內(nèi)合成代謝，以增加玉米抗旱性[8]。賈曉燕等基于生理生化指標(biāo)研究得出，細(xì)胞中SOD、POD、CAT、MDA、脯氨酸積累和丙二醛含量變化與作物抗旱性相關(guān)的結(jié)論[9]。Ali Noman等發(fā)現(xiàn)，干旱脅迫增強(qiáng)了玉米植株SOD和POD的活性以改善玉米植株的抗氧化防御系統(tǒng)，減輕干旱脅迫對玉米植株的危害[10]。前人對玉米抗旱性研究多涉及生理生化指標(biāo)，但較少涉及從RNA水平研究玉米抗旱基因的表達(dá)情況。

而隨著后基因組時代的到來，轉(zhuǎn)錄組測序技術(shù)(RNA-seq)廣泛應(yīng)用于醫(yī)學(xué)[11]、藥物研發(fā)[12]和農(nóng)業(yè)科學(xué)[13]等基礎(chǔ)研究[14]。對玉米轉(zhuǎn)錄組序列的結(jié)構(gòu)進(jìn)行解釋、對轉(zhuǎn)錄組的功能進(jìn)行分析，能從分子水平揭示玉米植株內(nèi)部的生物抗性[15]。目前RNA-seq應(yīng)用于玉米相關(guān)性研究較多；谷榮歡利用RNA-seq發(fā)現(xiàn)玉米K 305 ms花藥與發(fā)育相關(guān)的代謝網(wǎng)絡(luò)中部分基因表達(dá)量不足，是導(dǎo)致雄配子形態(tài)和發(fā)育過程異常的原因[16]。李川等報道，花粒期的玉米在高溫脅迫下，熱激蛋白、熱激轉(zhuǎn)錄因子、生長素基因和磷脂?；?4,5-二磷酸基因顯著表達(dá)[17]。王其等研究發(fā)現(xiàn)，將杜仲內(nèi)生拮抗細(xì)菌DZSY 21引入玉米，可調(diào)節(jié)玉米葉片中抗性相關(guān)基因的表達(dá)，適應(yīng)生存環(huán)境的變化，減少傷害[18]。

但是利用轉(zhuǎn)錄組測序分析對比PEG脅迫前后玉米基因差異表達(dá)的研究較少報道。因此本研究在RNA-seq的基礎(chǔ)上，通過對玉米轉(zhuǎn)錄組進(jìn)行一系列的數(shù)據(jù)比對、分析，初步解釋玉米抗旱應(yīng)答機(jī)制，以期為玉米的抗旱性研究及玉米抗旱基因的發(fā)掘奠定了基礎(chǔ)。

1 材料與方法

1.1 實(shí)驗(yàn)材料

試驗(yàn)材料為玉米自交系鄭58。

1.2 實(shí)驗(yàn)設(shè)計

采用人工控水盆栽試驗(yàn)法：將玉米移栽在花盆中(花盆直徑30 cm，高25 cm，每盆裝15 kg培養(yǎng)基質(zhì))生長，盆中含有培養(yǎng)基質(zhì)(3園土∶1珍珠巖∶1泥炭土)，一次性施入3 g復(fù)合肥(N∶P2O5∶K2O=15∶9∶10)，種前保持土壤含水量不低于85%以保證出苗。設(shè)置樣品S 1處于正常供水和S 3處于干旱脅迫2個狀態(tài)，3次重復(fù)，對照組玉米保持整個生長期正常澆水的濕潤狀態(tài)，而干旱處理從苗期(5葉期)開始控制灌溉20%的PEG-6000滲透調(diào)節(jié)劑模擬干旱脅迫[19]。

1.3 RNA提取

改進(jìn)的SDS酚法?。?/p>

1) 取0.1 g新鮮玉米葉片，液氮研磨成粉末狀；

2) 轉(zhuǎn)入1.5 mL離心管中，迅速加入600μL SDS提取液(2% SDS，25 mmol·L-1EDTA，100 mmol·L-1Tris-HCl，2%PVP)和40μL巰基乙醇，混勻后再加入150μL無水乙醇和150μL 3 mol·L-1的KAc，混勻；

3) 在轉(zhuǎn)速12 000 r·min-1、4 ℃下離心15 min，取上清，加入250μL水飽和酚和250μL氯仿，混勻后冰上放置10 min；

4) 在轉(zhuǎn)速12 000 r·min-1、4 ℃下離心15 min，取上清，加入等體積的氯仿抽提，混勻后冰上放置10 min；

5) 在轉(zhuǎn)速12 000 r·min-1、4 ℃下離心15 min；取上清，加入等體積異丙醇，混勻后冰上放置10 min；

6) 在轉(zhuǎn)速12 000 r·min-1、4 ℃下離心20 min，棄上清，加入75%的乙醇1 mL，溫和洗滌沉淀；

7) 在轉(zhuǎn)速12 000 r·min-1、4 ℃下離心5 min，棄上清，干燥后用50μL DEPC水溶解RNA[20]。

1.4 文庫構(gòu)建及RNA-Seq測序

樣品檢測合格后，先采用Oligo(dT)方法分離并純化出mRNA，并將mRNA打斷為片段，通過RT-PCR擴(kuò)增得到玉米轉(zhuǎn)錄組cDNA文庫，再采用Illumina PE 150測序平臺技術(shù)對構(gòu)建好的文庫進(jìn)行高通量測序。

表1 轉(zhuǎn)錄組測序數(shù)據(jù)及質(zhì)量情況統(tǒng)計表

1.5 數(shù)據(jù)的處理

使用軟件：fastx-toolkit-0.0.14對每一個樣本的原始測序數(shù)據(jù)進(jìn)行測序相關(guān)質(zhì)量評估；使用StringTie-1.3.3 b軟件對每個樣本進(jìn)行拼接、組裝；使用軟件DEGseq 2-1.24.0進(jìn)行組間差異表達(dá)分析，默認(rèn)參數(shù)：p-adjust<0.001 &|log 2 FC| ≥1；使用BLAST 2 GO-2.5與GO數(shù)據(jù)庫進(jìn)行序列比對；使用DIAMOND-0.9.24將轉(zhuǎn)錄本與NCBI-NR，Swiss-Prot，COG數(shù)據(jù)庫進(jìn)行序列比對；采用HMMER-3.2.1與Pfam數(shù)據(jù)庫比對；采用軟件Goatools-0.6.5進(jìn)行GO功能富集分析；使用KOBAS-2.1.1獲得到KEGG Orthology結(jié)果。

2 結(jié)果與分析

2.1 玉米轉(zhuǎn)錄組測序數(shù)據(jù)及質(zhì)量情況

運(yùn)用統(tǒng)計學(xué)的方法對所有測序reads的每個cycle進(jìn)行堿基分布和質(zhì)量波動的統(tǒng)計；從表1可知，對測序數(shù)據(jù)過濾后S 1、S 3分別得到45 426 908、48 371 948條clean reads，且每個樣本的clean bases均在6.78 Gb以上；堿基百分比Q 20均大于98.26%、堿基百分比Q 30均大于94.65%，其GC含量在55.9%～56.09%之間，含量接近；堿基測序錯誤率在0.024 1%～0.024 5%的較低水平；以上數(shù)據(jù)直觀地反映出測序樣本的測序質(zhì)量和文庫構(gòu)建質(zhì)量都比較好，能用于后續(xù)的分析。

2.2 轉(zhuǎn)錄本組裝及功能注釋情況

使用StringTie軟件對2個樣本進(jìn)行拼接，最終merge在一起，結(jié)果見圖1，在該轉(zhuǎn)錄本長度范圍內(nèi)共得到161 823個轉(zhuǎn)錄產(chǎn)物(transcipts)，轉(zhuǎn)錄本長度主要集中在600～3 000 bp之間；利用生物信息學(xué)方法將轉(zhuǎn)錄本組裝獲得的轉(zhuǎn)錄本與六大數(shù)據(jù)庫NCBI-NR、Swiss-Prot、Pfam、COG、GO和KEGG進(jìn)行比對，全面獲得轉(zhuǎn)錄本的功能信息并對各數(shù)據(jù)庫注釋情況進(jìn)行統(tǒng)計。從表2可知，共有138 270個轉(zhuǎn)錄本得到注釋，在e-value≤1 e-5條件下，Nr注釋成功的transcript數(shù)目有133 087條，占組裝得到總transcript的96.25%；GO注釋成功的transcript數(shù)目有108 300條，占78.33%；Swiss-port注釋成功的transcript數(shù)目有109 166條，占78.95%；COG注釋成功的transcript數(shù)目有130 072條，占94.07%；KEGG注釋成功的transcript數(shù)目有65987條，占47.72%；Pfaw注釋成功的transcript數(shù)目有99 546條，占71.99%。

圖1 轉(zhuǎn)錄本組裝長度分布統(tǒng)計

表2 轉(zhuǎn)錄本功能注釋情況統(tǒng)計

2.3 玉米轉(zhuǎn)錄組響應(yīng)干旱的差異表達(dá)分析

利用DEGseq軟件，參數(shù)設(shè)置：p-adjust<0.001 &|log 2 FC|≥1，以未進(jìn)行干旱脅迫的玉米轉(zhuǎn)錄組S 1為參考，與轉(zhuǎn)錄組S 3進(jìn)行數(shù)據(jù)比對；從圖2可知，玉米轉(zhuǎn)錄組共獲得12 358個差異表達(dá)基因，其中上調(diào)基因5 275個，下調(diào)基因7 083個。

2.4 差異表達(dá)基因的COG功能分析

對差異表達(dá)基因進(jìn)行COG功能分類，結(jié)果見表3，共得到23個不同的COG功能，2 127個轉(zhuǎn)錄本參與細(xì)胞過程和信號傳導(dǎo)的9個COG功能，其中翻譯后修飾，蛋白質(zhì)轉(zhuǎn)換，伴侶蛋白681個，占7.74%；信號傳導(dǎo)機(jī)制功能有轉(zhuǎn)錄本675個，占7.67%；非細(xì)胞運(yùn)輸、分泌和囊泡運(yùn)輸391個，占4.44%。1 337個轉(zhuǎn)錄本歸納到新陳代謝類的8個COG功能，其中碳水化合物運(yùn)輸和代謝315個，占3.58%；氨基酸運(yùn)輸和代謝265個，占3.01%。1 296個轉(zhuǎn)錄本參與信息存儲與處理的5個功能，轉(zhuǎn)錄506個，占5.75%；翻譯，核糖體結(jié)構(gòu)和生物發(fā)生356個，占4.05%。而分類到POORLY CHARACTERIZED未知功能的轉(zhuǎn)錄本有4 039個(45.91%)。

表3 差異表達(dá)基因COG功能分類

圖2 基因差異表達(dá)火山圖分析

2.5 干旱脅迫下玉米差異表達(dá)基因的分析

2.5.1玉米轉(zhuǎn)錄組差異表達(dá)基因GO功能富集分析

采用軟件Goatools對干旱脅迫處理后玉米的差異表達(dá)基因集中的轉(zhuǎn)錄本進(jìn)行GO功能富集分析。結(jié)果見圖3，可以看出，玉米轉(zhuǎn)錄組差異表達(dá)基因富集在生物過程、細(xì)胞成分、分子功能三大類；差異表達(dá)基因富集的前20亞類中生物過程大類中有6個亞類，分別是光合作用、mRNA加工、tRNA代謝過程、ncRNA代謝過程、細(xì)胞表面受體信號通路、代謝過程。還富集在細(xì)胞成分大類中的8個亞類，分別為葉綠體、質(zhì)體、細(xì)胞質(zhì)、質(zhì)體部分、葉綠體的部分、胞質(zhì)部分、細(xì)胞內(nèi)的部分、細(xì)胞部分。分子功能大類的6個亞類得到富集，分別是GTP活性、核糖核酸結(jié)合、三磷酸鳥苷結(jié)合和鳥苷核糖核苷酸結(jié)合、鳥苷核苷酸結(jié)合、嘌呤核苷綁定。

圖3 差異表達(dá)基因GO功能富集分析

2.5.2玉米轉(zhuǎn)錄組差異表達(dá)基因KEGG功能富集分析

Pvalue corrected小于0.05時認(rèn)為顯著富集，基于此，對玉米轉(zhuǎn)錄組參與KEGG代謝通路進(jìn)行富集，結(jié)果見圖3，玉米轉(zhuǎn)錄組參與KEGG代謝通路分為三大分支：代謝、遺傳信息處理和細(xì)胞過程，共獲得28個標(biāo)準(zhǔn)KEGG代謝通路。其中代謝分支中富集到的通路有光合作用-天線蛋白、光合生物的固碳作用、卟啉和葉綠素代謝、脂肪酸降解、精氨酸和脯氨酸代謝、磷酸肌醇信號、糖酵解和糖異生、丙氨酸代謝、氨基糖和核苷酸糖代謝、淀粉和蔗糖代謝、氰氨基酸代謝、光合作用、類固醇生物合成、?；撬岷痛闻；撬岽x、精氨酸生物合成、花生四烯酸代謝、纈氨酸、亮氨酸和異亮氨酸的降解、硫代謝、亞麻酸代謝、纈氨酸、亮氨酸和異亮氨酸生物合成、亞油酸的新陳代謝、煙酸和煙酰胺代謝、苯丙合成。遺傳信息處理富集到的通路有氨酰生物合成、蛋白酶體、蛋白輸出、剪接體。細(xì)胞過程富集到的通路有吞噬體。

圖4 差異表達(dá)基因KEGG富集分析

3 結(jié)果與討論

近年來，轉(zhuǎn)錄組測序技術(shù)的應(yīng)用很大程度上推動了植物分子水平方面的研究[21]。本研究對正常灌溉和PEG模擬干旱脅迫下的玉米進(jìn)行RNA-sep測序，玉米轉(zhuǎn)錄組共獲得12 358個差異表達(dá)基因，其中上調(diào)基因5 275個，下調(diào)基因7 083個；可見，抗旱基因?yàn)樯险{(diào)表達(dá)，干旱脅迫限制了玉米的某些代謝活動；下調(diào)表達(dá)高于上調(diào)表達(dá)；這與王偉東[22]、吉福桑[23]和楊文軍等[24]對植物干旱脅迫下的基因表達(dá)上下調(diào)結(jié)果一致。對差異表達(dá)基因進(jìn)行COG功能分類，共得到23個不同的COG功能，細(xì)胞過程和信號傳導(dǎo)中的9個COG功能、新陳代謝類的8個COG功能、信息存儲與處理的5個功能以及POORLY CHARACTERIZED的未知功能；除未知功能的轉(zhuǎn)錄本4 039個外；翻譯后修飾，蛋白質(zhì)轉(zhuǎn)換，伴侶蛋白681個，占7.74%；信號傳導(dǎo)機(jī)制功能有轉(zhuǎn)錄本675個，占7.67%；轉(zhuǎn)錄506個，占5.75%；非細(xì)胞運(yùn)輸、分泌和囊泡運(yùn)輸391個，占4.44%。碳水化合物運(yùn)輸和代謝315個，占3.58%；翻譯，核糖體結(jié)構(gòu)和生物發(fā)生356個，占4.05%；氨基酸運(yùn)輸和代謝265個，占3.01%。植物可通過調(diào)整碳水化合物的合成利用來滿足植物對能量的需求，以緩解干旱脅迫[25]，同時干旱脅迫會增強(qiáng)玉米植株SOD等酶的活性以改善玉米植株的抗氧化防御系統(tǒng)[10]，而SOD的合成表達(dá)又由其伴侶蛋白來實(shí)現(xiàn)[26]。由此可見，玉米在PEG脅迫下通過信號傳導(dǎo)環(huán)境變化，進(jìn)而調(diào)控轉(zhuǎn)錄翻譯以合成SOD等酶，提高玉米適應(yīng)干旱環(huán)境的抗性。

差異表達(dá)基因GO功能富集分析，說明玉米轉(zhuǎn)錄組特性與生物過程、細(xì)胞成分、分子功能相關(guān)；進(jìn)一步對KEGG代謝通路富集，在Pvalue corrected小于0.05時顯著富集到28條標(biāo)準(zhǔn)KEGG代謝通路，這些基因主要參與光合作用-天線蛋白、光合生物的固碳作用、卟啉和葉綠素代謝、脂肪酸降解、精氨酸和脯氨酸代謝、磷酸肌醇信號、糖酵解和糖異生等生命代謝活動。其中多條代謝通路涉及光合作用等相關(guān)途徑，說明干旱脅迫對光合作用產(chǎn)生了顯著影響，劉素軍等報道，干旱脅迫會加強(qiáng)光合作用-天線蛋白途徑向相連的反應(yīng)中心輸送能量，從而響應(yīng)水分脅迫環(huán)境[27]。孫琪研究發(fā)現(xiàn)，光合生物的固碳作用用于合成更多的有機(jī)物滿足植物對能量的需求，以此來提高抗鹽能力[28]。張國偉研究發(fā)現(xiàn)，脂肪酸含量是影響光合作用的重要因素之一[29]。Trovato M等發(fā)現(xiàn)，脯氨酸的積累是生物有機(jī)體對環(huán)境脅迫產(chǎn)生的一種應(yīng)答反應(yīng)，能夠增強(qiáng)植物的抗逆能力[30]。王建等報道，精氨酸代謝會間接促進(jìn)PA、NO的合成，從而提高番茄對Cu脅迫的抵抗能力[31]。Ishitani M等發(fā)現(xiàn)，磷酸肌醇信號途徑在高等植物的生長、發(fā)育及對環(huán)境脅迫反應(yīng)中發(fā)揮重要作用[32]。上述代謝通路都與植物的抗逆性相關(guān)，但個別代謝通路在玉米中的調(diào)控機(jī)理暫不清楚，這也為進(jìn)一步了解玉米抗旱機(jī)理指明了研究方向。

RNA-Seq是一種在轉(zhuǎn)錄水平上更為精確的測定分析方法，其通量高、成本低、分辨率高、靈敏度高且不受物種限制[33-35]。本實(shí)驗(yàn)采用RNA-Seq研究玉米轉(zhuǎn)錄組，并進(jìn)行數(shù)據(jù)分析工作，其結(jié)果極大豐富了玉米轉(zhuǎn)錄組數(shù)據(jù)資源，為探討玉米抗旱機(jī)理打下了堅實(shí)的基礎(chǔ)。