朱 琳，袁 夢，高紅秀，曲 悅，何沈雨，王 劍，王 珊， 王思宇，楊 玲，張忠臣，金正勛

(東北農業(yè)大學 農學院，黑龍江 哈爾濱 150030)

低溫是影響農作物生產的重要危害因素之一。在前期育苗過程中，水稻經常會受到低溫影響，導致水稻生長遲緩或停止生長[1-3]。早稻直播面積的逐漸擴大，水稻苗期受低溫脅迫的可能性日益增加[4]。因此，提高水稻抗低溫能力、深入研究水稻耐低溫機制成為許多研究人員關注的熱門話題[5]。水稻是世界上最重要的糧食作物之一，全世界超過1/3的人口以水稻為主食[6]。然而，水稻產量受到包括低溫等不同非生物脅迫的影響很大[7]。水稻苗期在遭受低溫脅迫后會打破細胞滲透壓穩(wěn)態(tài)、調節(jié)滲透壓的物質積累及代謝異常等，導致葉片失綠卷曲、生長緩慢、育性下降和產量降低等[8]。在農業(yè)生產中，低溫脅迫會致使幼苗枯黃、秧苗腐爛等現象，影響其穗發(fā)育和產量等[9]。低溫冷害是水稻種植所面臨的主要困難之一[10]。目前，關于作物耐冷性研究主要停留在生理階段[11]。水稻抗逆相關功能基因組學的研究一直是水稻研究核心內容。分析低溫脅迫下幼苗體內代謝物變化特征和研究水稻苗期耐低溫機理有助于解析水稻應對低溫脅迫的生理機制與分子機制，為以后水稻分子育種奠定基礎，這對農業(yè)經濟和社會都有著深遠意義。水稻可以通過抵御代謝反應、改變體內激素的含量及活性來調節(jié)生理過程，使其在低溫條件下能維持正常的生理活動[12-13]。轉錄水平的調控是生物體最重要的調控方式，與基因組學相比具有針對性強、研究范圍小的優(yōu)點[14]。因此，被認為是一種在轉錄水平上更精準的測定分析方法[15]。目前，利用高通量轉錄組測序技術已發(fā)現了細枝木麻黃[16]、玉米[17]、西葫蘆[18]、油菜[19]等植物與耐低溫相關的轉錄、病害防御、信號轉導、轉運、次級代謝等大量相關基因。

本研究采用高通量測序方法-轉錄組測序(RNA-Seq)技術，在轉錄組水平上研究水稻苗期應答低溫差異基因的表達特性，為增強對水稻耐低溫機理的認識，明確水稻應答低溫的分子機制，以及培育水稻耐低溫品種提供參考。

1 材料和方法

1.1 試驗材料

供試材料為粳稻品種中花-11，由浙江大學生命科學學院植物生理與生化國家重點實驗室提供。

1.2 試驗方法

1.2.1 培養(yǎng)條件 種子催芽方法和水培方法參考國際水稻所的水稻水培方法[20]。幼苗長至一葉一心時，移栽到規(guī)格為6 cm×6 cm×6 cm的培養(yǎng)缽中，三葉一心時小心將生育進程一致的幼苗取出避免傷到根系。將其中1/2的幼苗放置在光溫培養(yǎng)箱(HPG-280，哈爾濱市東聯電子技術開發(fā)有限公司)作為對照組(WTCK-2)，培養(yǎng)條件為：晝夜溫度為28 ℃/22 ℃，12 h/12 h，相對濕度為80%；另1/2的幼苗作為低溫處理組(WTCOLD-2)，移入另一光照培養(yǎng)箱(HPG-280，哈爾濱市東聯電子技術開發(fā)有限公司)，晝夜溫度保持在17 ℃，12 h/12 h，相對濕度80%，保持10 d。

1.2.2 RNA提取、文庫構建 水稻幼苗在低溫處理后，分別取對照和處理后的苗期葉片放入液氮速凍，存至-80 ℃冰箱保存?zhèn)溆?。委托北京諾禾致源生物信息科技有限公司完成RNA的提取、質控、建庫及Illumina HiSeqTM4000測序工作。

1.2.3 原始數據的組裝 由高通量測序(Illumina HiSeqTM)得到的原始數據文件，利用CASAVA堿基識別(Base calling)分析轉化成原始測序序列(Sequenced reads)，即Raw reads。將低質量和帶接頭的reads進行過濾得到Clean reads，從而提高信息分析質量，Clean reads將作為后續(xù)分析的研究基礎。本研究使用Q30(堿基準確率達到99.9%)的質量控制標準。

1.2.4 差異表達基因的篩選 通過使用DESeq自帶的標準化方法對原始數據進行標準化處理。在差異分析過程中，使用負二項分布法對Read count的分布進行估計，評估且計算完Pvalue后，對Pvalue進行多重假設檢驗校正，來減少假陽性。差異基因篩選標準為：|log2(Fold Change)|>1并且校正的P值(Padj)<0.05。

1.2.5 GO基因功能注釋分析 使用GOseq[21]將組裝后的Unigene序列按照Corrected_Pvalue<0.05篩選差異基因，為尋找出基因顯著富集的GO Term，對篩選的差異基因進行GO基因功能注釋分析，獲得所有差異表達基因的功能注釋、生物學功能以及相關的代謝途徑，從而推測這些基因對應的特定過程。

1.2.6 KEGG pathway分析 使用軟件KOBAS(2.0)進行特定基因的通路定位，找出顯著富集的通路，最后對KEGG注釋的基因功能信息進行生物學通路的注釋與預測，解釋低溫脅迫下植物對應的調控機制。

1.2.7 聚類分析 使用聚類軟件Pheatmap對差異基因的相對表達水平值log2(ratios) 進行聚類。以FPKM值的讀數來計算差異表達水平，顏色從紅到藍，表示log10(FPKM+1)從大到小。使用與之匹配的距離算法，計算出基因間的距離，反復迭代，算出基因之間的相對距離，最后依據基因的相對距離遠近來分成不同的Subcluster，從而實現聚類。

1.2.8 qRT-PCR驗證 利用TRIzol Reagent提取對照組和處理組葉片的總 RNA，cDNA合成利用PrimeScriptTMRT reagent Kit with gDNAEraser(Perfect Real Time)，熒光定量檢測使用abm?EvaGreen qPCR MasterMix-no dye試劑盒。利用PrimerQuest Tool(http：//sg.idtdna.com/Primerquest/Home/Index)設計引物，優(yōu)先選擇靠近基因3′端的引物并進行 Real-time PCR 分析。反應程序為95 ℃，10 min；反應循環(huán)數為45次，循環(huán)過程是94 ℃，15 s，60 ℃，1 min；72 ℃，20 s。用ABI StepOne Plus(ABI 公司，美國)進行熒光定量檢測。利用2-ΔΔCt法，以對照組中各基因表達量為1，Log2FC表示log2倍，誤差線表示標準偏差，以OsActin為內參基因對目的基因進行轉錄水平的相對定量分析。

1.2.9 新基因功能預測 將測序所得的全部reads數據和已知的基因模型進行比較，對比原有的基因注釋文件，可以發(fā)現新基因。根據轉錄組數據所提供的新基因堿基序列信息，在MSU(http：//rice.plantbiology.msu.edu/analyses_search_locus.shtml)網站上查找相關需要信息。

1.2.10 蛋白互作預測 應用STRING蛋白質互作數據庫(http：//string-db.org/)中的互作關系分析差異基因蛋白互作網絡。利用數據庫中提取目標基因集互作關系和目標基因集中的序列，應用BlastX比對到String數據庫中包含的參考物種的蛋白質序列上來構建蛋白質互作關系。

2 結果與分析

2.1 低溫影響水稻幼苗期的生長發(fā)育

17 ℃低溫處理10 d后，水稻在正常條件和低溫處理下，幼苗長勢差異明顯(圖1)。低溫脅迫后的幼苗葉片偏黃、葉尖變黃枯且株高較正常條件下偏低。表明低溫脅迫嚴重抑制了水稻幼苗的生長發(fā)育。

圖1 正常條件(左)與低溫處理(右)下的中花-11表型Fig.1 Zhonghua-11 phenotype under normal condition (left) and low temperature (right) treatment

2.2 測序結果的質量分析

分別收集處理組和對照組水稻苗期葉片并制備其RNA樣品。這些樣品通過illumina HiSeqTM測序平臺進行測序，一共生成約1.01億條長度約100 bp的原始數據，2種條件下的原始片段分別是47 576 488，53 790 724條。根據生物信息學分析，舍棄低質量的讀數，經過序列拼接，最后獲得約0.99億干凈的高質量序列，分別得到46 384 074，52 483 052條干凈的高質量序列，占原始數據的97.53%。其中，平均Q20、Q30都分別可達到96%，91%以上，GC含量分別為54.91%和55.65%，不確定堿基比例僅為0.02%。使用TopHat2將過濾后的數據定位到中花-11(對照組和處理組)的基因組。在2個樣本組的高質量序列讀數中，在參考序列上有唯一比對位置的分別占86.29%和84.95%、有多個比對位置的占1.57%和1.75%、總測序序列定位數占87.86%，86.70%。結合以上數據可以說明轉錄組測序質量較高，符合進一步的生物信息學分析(表1)。

表1 測序數據統(tǒng)計Tab.1 Sequencing date statistics

注：Q20.質量值≥20的堿基所占比例；Q30.質量值≥30的堿基所占比例；GC含量.計算堿基G和C的數量總和占總的堿基數量的百分比。

Note：Q20.The proportion of the bases whose quality value≥20； Q30.The proportion of the bases whose quality value≥30； GC content.Calculate the percentage of the total number of bases G and C as a percentage of the total number of bases.

2.3 差異表達基因的GO功能注釋分析

通過GO功能富集分析，對照組與處理組存在2 044個差異表達基因(Differentially expressed genes，DEGs)，其中，1 252個基因表達上調，792個基因表達下調。被富集到3 405個GO terms中。其中生物學過程(Biological process，BP)占63.29%、細胞組分(Cellular component，CC)占8.81%、分子功能(Molecular function，MF)占27.90%。通過數據所占比例可以初步判斷大多數差異表達基因與一些生物學功能顯著相關。單一生物體代謝過程(9.55E-16)、質體(3.94E-21)和陽離子結合(4 .02×10-3)分別是BP、CC、和MF中最顯著富集的功能類別。低溫脅迫下DEGs涉及的20個最顯著富集的功能如圖2所示。

圓點所在的位置代表著富集的條目，圓點的大小表示差異基因的數量。 The position of the dot represents the enriched item，the size of the dot indicates the number ofDifferential genes.

2.4 低溫脅迫下差異表達基因的鑒定與篩選

通過將表達模式相同或相近的基因聚集成類，用以推測已知基因的新功能或未知基因的功能[22]。為了鑒定具有功能富集的聚類，利用差異基因的相對表達水平值log2(ratios)進行層次聚類(圖3)，從圖3可以看出，該材料在低溫脅迫下基因表達量有變化。

基于以上發(fā)現，在對照組和處理組中鑒定了2 044個差異表達基因，其中有關于光合作用、苯丙氨酸代謝、WRKY轉錄因子等多種功能途徑。其中右側虛線的右側圓點代表差異表達上調的基因有1 252個，左側虛線的左側圓點代表差異表達下調的基因有792個，兩側虛線中間表示無顯著性差異表達的基因(圖4)。

2.5 差異表達基因的KEGG通路分析

為了精準分析低溫脅迫對水稻發(fā)育過程中參與的代謝途徑發(fā)生的變化，對得到的Unigene進行KEGG代謝通路分析(表2)。結果表明，在代謝通路中，主要是光合作用天線蛋白(ko00196)、光合作用(ko00195)、類胡蘿卜素生物合成(ko00906)、次生代謝產物的生物合成(ko01110)。利用KEGG數據庫對差異基因進行數據分析，得到差異基因KEGG注釋通路圖，其中水稻苗期對受低溫脅迫影響較大的代謝通路分析如下：

通路一：光合作用

低溫脅迫下，此通路中的光合系統(tǒng)Ⅰ存在11個差異表達基因，且全部表現為下調；光系統(tǒng)Ⅱ中有24個差異表達基因，其中3個基因表現為上調，21個基因表現為下調；細胞色素b6/f復合體中含有1個差異表達基因且表現為下調；光合電子運輸中有20個差異表達基因，有8個基因表現為上調，其余基因均表現為下調；F型ATP酶中有6個表現為下調基因(圖5)。上述結果得出，在光合作用通路里這些差異表達基因可能是受低溫處理影響較大的基因。

WTCK_2.對照條件下的中花-11；WTCOLD_2.低溫脅迫下的中花-11。 WTCK_2.Zhonghua-11 under control conditions； WTCOLD_2.Zhonghua-11 under low temperature stress.

空心箭頭.上調的顯著差異表達基因；實心箭頭.下調的顯著差異表達基因；虛線空心箭頭.無顯著性差異表達的基因；橫坐標.基因在不同樣本中表達倍數變化；縱坐標.基因表達量變化差異的統(tǒng)計學顯著性。

Hollow arrow.Significant Different expressed genes up-regulated； Solid arrow. Significantly Different expressed genes down-regulated； Dotted hollow arrow. Genes with no significant
Differential expression Abscissa. Gene expression fold change in Different samples； Ordinate.Statistically significant difference in gene expression changes.

圖4差異基因的火山
Fig.4Volcanicmapofdifferentialgenes

PhotosystemⅡ.光系統(tǒng)Ⅱ；PhotosystemⅠ.光系統(tǒng)Ⅰ；Cytochrome b6/f complex.細胞色素b6/f復合體；Photosynthethic electron transport.光合電子運輸；F-type ATPase. F型ATP酶；Carbon fixation in photosynthetic organisms.光合生物中的碳固定；Chloroplast stroma.葉綠體基質；Thylakoid membrane.類囊體膜；Thylakoid lumen.類囊體腔。

圖5差異表達基因參與的光合作用代謝途徑
Fig.5Photosynthesismetabolicpathwaysinvolvedindifferentiallyexpressedgenes

通路二：苯丙氨酸代謝

苯丙氨酸代謝途徑是植物最重要的次生代謝途徑之一，苯丙氨酸解氨酶(Phenylalanine ammonia-lyase，PAL)基因表達量高，與之對應的次生代謝產物產量就高，水稻幼苗的整體抗逆能力就強[23-24]。低溫脅迫下，在此通路中發(fā)現17個差異表達基因，8個DEGs表達上調，其余9個DEGs表達下調。分別在4.3.1.24(K01110)和4.3.1.25(K01110)處PAL相關基因OS02G0626400(1.654 5)、OS02G0626100(1.229 9)、OS04G0518100(1.544 4)上調表達，進而調節(jié)水稻次生代謝反應(圖6)。從這個通路圖可以看出，在低溫脅迫下PAL基因表達上調，說明低溫與苯丙氨酸代謝通路具有相關性，差異性表達的基因可能是低溫應答的關鍵基因。

圖6 差異表達基因參與苯丙氨酸代謝通路Fig.6 Differentially expressed genes involved in phenylalanine metabolic pathway

序號No.通路名稱Pathway nameko_ID差異基因數DEGsP-value1 光合作用-天線蛋白質ko00196121.26×10-62 光合作用ko00195212.06×10-63類胡蘿卜素的生物合成ko00906122.88×10-44次生代謝產物的生物合成ko011101256.33×10-35代謝途徑ko011002091.22×10-26乙醛酸和二羧酸代謝ko00630122.59×10-27植物激素信號轉導ko04075314.66×10-28晝夜節(jié)律-植物ko0471285.60×10-29氮代謝ko0091075.63×10-210醚脂質代謝ko0056566.23×10-211抗壞血酸和aldarate新陳代謝ko0005387.58×10-212碳固定在光合生物中ko00710138.23×10-213亞油酸代謝ko0059148.47×10-214半胱氨酸和蛋氨酸代謝ko00270159.07×10-215α-亞麻酸代謝ko0059289.92×10-216淀粉和蔗糖代謝ko00500241.21×10-117其他聚糖降解ko0051141.30×10-118卟啉和葉綠素代謝ko0086071.36×10-119苯丙素生物合成ko00940261.41×10-120葉酸生物合成ko0079041.64×10-1

表2(續(xù))

表2(續(xù))

2.6 qRT-PCR驗證基因表達

為了驗證轉錄組分析結果中差異基因表達的可靠性，上調基因中選取了一個與植物病原體相互作用的基因OS05G0502000、一個與氮代謝相關的基因OS02G0770800、2個參與苯丙氨酸代謝的基因OS01G0869800、OS07G0677200；下調基因中挑選了4個與光合作用有關的基因OS07G0558400、OS01G0501800、OS03G0592500、OS02G0581100，利用qRT-PCR方法進行驗證(表3)。使用RNA-Seq和qRT-PCR比較基因表達水平，以對照組中各基因表達量為1，統(tǒng)計低溫處理下各基因的相對表達量。變化結果表明，qRT-PCR驗證的基因表達量變化趨勢和RNA-Seq差異基因表達趨勢保持一致，如圖7所示。表明轉錄組分析所得的數據具有較高的可靠性。

表3 實時定量PCR引物Tab.3 Primers used for qRT-PCR analysis

圖7 RNA-Seq和qRT-PCR基因表達水平比較Fig.7 Comparison of RNA-Seq and qRT-PCR gene expression levels

2.7 新基因功能預測

綜合各種外顯子預測的算法和人們對基因結

構信息的認識，可預測出可能的完整基因[25]。新基因的預測有助于分析的完整性。轉錄組數據分析發(fā)現了25個新基因，它們的長度在779～6 551 bp不等，主要出現在1，2，3，4，5，6，7，8，10，12號染色體上，有6個基因沒有具體的信息描述，其他的主要表現在線粒體前體、鈣轉運ATP酶、氯通道蛋白和轉綠因子WRKY等方面(表4)。

水稻的WRKY、MYB、bHLH、bZIP、C2H2等家族調控水稻的生長發(fā)育和耐逆抗病過程[26-28]。本試驗發(fā)現，2個WRKY轉錄因子，即Novel00699和Novel00786。這2個WRKY轉錄因子可能參與了應答低溫環(huán)境的信號刺激，被誘導表達，提高作物抗逆能力。

表4 新基因功能預測Tab.4 Novel gene function prediction

2.8 預測的互作蛋白

本研究分析得到5個與光合作用有關的基因，分別為LOC_Os01g31690、LOC_Os03g56670、LOC_Os04g33830、LOC_Os07g04840、LOC_Os08g44680，預測有16個基因與其互作；一個與PAL有關的基因LOC_Os02g41630，且預測到有3個基因與其互作(表5)。這些結果有助于對水稻幼苗響應低溫脅迫的調節(jié)機制的進一步研究。

表5 預測的互作蛋白Tab.5 Predicted interaction proteins

3 結論與討論

3.1 結論

本研究對正常條件和低溫脅迫下水稻苗期葉片進行轉錄組測序。通過GO功能注釋、KEGG Pathway、聚類分析等方法在正常條件和低溫脅迫下進行差異基因表達水平上的比較，篩選出2 044個DEGs，其中，1 252個基因表達上調，792個基因表達下調；發(fā)現低溫脅迫后差異表達基因主要參與光合作用、苯丙氨酸代謝等途徑；在低溫脅迫下，水稻幼苗的光合作用受到影響，其中該作用中光系統(tǒng)Ⅱ與光合運輸電子受到的影響較大，基因變化趨勢主要體現為下調；并且低溫誘導了PAL活性升高，促進苯丙烷類次生代謝產物的合成，增強水稻幼苗在低溫脅迫下的耐受力。通過新基因預測發(fā)現，差異表達基因可能參與轉錄因子的調控。本研究提供了水稻苗期低溫研究方面一定的表達數據，同時也有助于對水稻苗期低溫響應的分子調控機理的進一步研究。

3.2 討論

國內外學者針對低溫脅迫對植物光合作用的影響有著大量研究。王春萍等[29]研究了低溫脅迫對水稻幼苗不同葉齡葉片葉綠素熒光特性，結果發(fā)現，低溫脅迫對水稻幼苗三葉期的光合利用能力造成了一定影響。PSⅡ是植物光合作用過程中進行光反應的重要結構，低溫脅迫會導致光能和光合系統(tǒng)Ⅱ(PSⅡ)電子傳遞轉化效率降低，造成水稻的凈光合速率、氣孔導度與蒸騰速率均明顯下降[30-32]。Kaniuga等[33-34]研究證明，低溫條件下降低了葉綠體的希爾反應活性和類PSⅡ的電子傳遞活性。本研究發(fā)現，低溫處理水稻幼苗后在光合作用通路上的差異表達基因，分別體現在光系統(tǒng)Ⅰ、光系統(tǒng)Ⅱ、細胞色素b6/f復合體、光合運輸電子和F型ATP酶5個方面，受溫度影響DEGs數量較多的為光系統(tǒng)Ⅱ和光合電子運輸2個部分。光系統(tǒng)Ⅱ中的DEGs有21個表現下調，3個表現上調；光合電子運輸中有20個差異表達基因，有8個表現上調，其余均為下調。通過這些數據可以補充說明：低溫脅迫下，水稻幼苗的光合作用受到影響，且該作用中光系統(tǒng)Ⅱ與光合運輸電子受到的影響較大，基因變化趨勢主要體現下調。

PAL在植物生長發(fā)育和抵御環(huán)境脅迫等方面具有重要的生理意義[35-36]。黃小貞等[24]發(fā)現在不同的組織和時間段PAL基因的表達具有特異性，并且受內部發(fā)育信號等的調節(jié)。董春娟等[37]研究得出，黃瓜幼苗中的PAL基因表達和酶活性可在10 ℃低溫條件下被誘導升高，證明了在低溫脅迫下，黃瓜中的PAL基因參與黃瓜幼苗對低溫脅迫的響應且存在基因功能的重疊。油菜在2 ℃處理后，葉片經PAL抑制劑氨基茚磷酸處理，葉片物質積累會發(fā)生遲緩，且Fv/Fm的值下降。證實PAL活性受到抑制會導致幼苗對低溫脅迫的抗性降低[38]。本試驗發(fā)現，水稻幼苗在低溫脅迫下，有17個DEGs參與苯丙氨酸代謝通路，其中有3個注釋到PAL功能，表達量均表現為上調。結合本試驗結果說明，水稻幼苗為在低溫條件下提高自身的抗逆性，保證其正常生長發(fā)育。所以，低溫誘導了PAL的活性升高，促進苯丙烷類次生代謝產物的合成，緩解低溫引起的不利因素。

眾多研究表明，擬南芥與水稻 WRKY 基因家族的生物學功能和復雜相關特征已被廣泛的挖掘并進行了功能鑒定[39-42]。目前，粳稻和秈稻分別鑒定出98，102個WRKY轉錄因子[41]，生物和非生物逆境因子可大量誘導WRKY基因表達[43-45]。鄂志國等[46]將水稻WRKY家族按WRKY結構域基本結構特征的差異為3類，且得出它們通過調控生長調節(jié)物質所介導的信號途徑，廣泛的影響水稻的生長發(fā)育、生物與非生物脅迫應答。彭喜旭等[47]研究表明，WRKY80可能通過茉莉酸/乙烯介導的信號途徑參與或調控防御反應。Shimono等[48-54]在水稻中找到WRKY45、WRKY13、WRKY03等多個WRKY家族基因被證明參與抗病防御調節(jié)。仇玉萍等[55]證實水稻中的7個WRKY因子(OsWRKY8、OsWRKY12、OsWRKY13、OsWRKY14、OsWRKY17、OsWRKY23 和OsWRKY45)能被4 ℃低溫誘導表達；Zhou等[56]研究發(fā)現GmWRKY21過量表達可提高擬南芥抵抗低溫的能力。本試驗研究發(fā)現，對照組和處理組的水稻幼苗共有10個WRKY轉錄因子發(fā)生差異性表達，有8個DEGs表現上調，分別為WRKY21(OS01G0821600)、WRKY24(OS01G0826400)、WRKY121(OS03G0741400)、WRKY68(OS04G0605100)、WRKY7(OS05G0537100)、WRKY8(OS05G0583000)、WRKY69(OS08G0386200)、WRKY97(OS12G0116400)；2個DEGs表現下調，分別是WRKY30(OS08G0499300)、WRKY90(OS09G0481700)。新基因預測中也發(fā)現了2個WRKY轉錄因子WRKY115(Novel00699)和WRKY117(Novel00786)，經過低溫處理后，水稻幼苗中部分WRKY轉錄因子的表達量發(fā)生變化，說明一些WRKY轉錄因子參與低溫脅迫。