孕穗期低溫脅迫下寒地水稻轉錄組分析

摘要：為了明確寒地水稻孕穗期低溫應答機制，以孕穗期耐冷性差異較大的龍粳25、龍粳11水稻品種為材料，對幼穗穎花開展轉錄組測序研究。結果表明，在龍粳11中共檢測到9 468個差異表達基因（DEG），在龍粳25中共檢測到6 784個DEG。其中3 508個DEG為二者共同表達，3 276個DEG在龍粳25中特異表達。GO富集分析表明，龍粳25、龍粳11的共有DEG主要富集在生物進程調控、氧化還原反應、非生物刺激應答等類別中，KEGG通路主要富集在植物激素信號轉導、碳代謝、淀粉蔗糖代謝等途徑。而龍粳25的特有DEG多富集在生物進程調控、非生物刺激應答、氧化還原反應進程等類別，KEGG通路主要富集在植物激素信號轉導、內質網蛋白加工、苯丙烷合成等途徑。綜上表明，寒地水稻對孕穗期的低溫應答是一個復雜過程，包含多個代謝調控途徑及多個相應基因。

關鍵詞：寒地；水稻；轉錄組；孕穗期低溫

中圖分類號：S511.01" 文獻標志碼：A

文章編號：1002-1302（2024）19-0034-07

收稿日期：2023-10-31

基金項目：黑龍江省博士后面上資助項目（編號：LBH-Z21027）；黑龍江省省屬科研業務費項目（編號：CZKYF2023-1-C007、CZKYF2021-2-C002）；黑龍江省農業科技創新跨越工程重大需求科技創新攻關項目（編號：CX23ZD01）；黑龍江省自然科學基金優秀青年項目（編號：YQ2021C033）。

作者簡介：郭震華（1985—），男，內蒙古集寧人，博士，副研究員，主要從事水稻分子育種研究。E-mail：hljsdsgzh@163.com。

通信作者：劉傳雪，碩士，研究員，主要從事水稻育種研究，E-mail：liuchuanxue2007@163.com；馮延江，博士，研究員，主要從事水稻育種栽培研究，E-mail：zixuanfeng2008@163.com。

水稻是全世界最重要的糧食作物之一^［1］，全球一半以上的人口以水稻為主食。水稻起源地多為熱帶或亞熱帶地區，因此對低溫十分敏感。通常認為，秈稻低于18 ℃，粳稻低于15 ℃，就會對其正常生長造成不利影響^［2］。目前，水稻冷害已然成為世界性問題，其中尤以日本、韓國及我國東北稻區最為突出^［3］。黑龍江省地處我國最北部，屬于寒地稻作區，是世界水稻栽培的北限，年平均氣溫為全國最低；低溫冷害是限制黑龍江水稻安全生產的重要因素。

水稻在全生育期內都可能受到低溫冷害的威脅。孕穗期是花粉形成的關鍵時期，若遭受低溫會直接影響花粉發育，造成花粉不育或花粉管無法萌發，導致結實率下降，最終造成水稻減產^［4］。具體來說，以減數分裂期對低溫最為敏感；若此時遭遇低溫，水稻花粉壁上的絨氈層細胞不能正常降解，影響淀粉在花粉中的正常轉運，花粉發育所需的能量供給不足，最終造成花粉敗育^［5-6］。

低溫脅迫除對結實率等直觀表型有影響外，對作物各方面生理功能都會造成一系列的影響，包括光合作用、能量代謝、滲透調節、過氧化物酶系統保護等^［7-8］。其中涉及到的不同代謝途徑中多個調控基因，均被報道在低溫脅迫應答中起到不同的作用。Ctb1則是第1個通過圖位克隆得到的水稻孕穗期耐冷調控基因，該基因編碼1個F-box蛋白，通過與E3泛素連接酶亞基Skp1間的互作，來啟動耐冷應答反應^［9］。Osmyb4被稱作水稻耐冷應答的總開關，低溫脅迫下正向調控下游多個基因的表達^［10］。過氧化物酶基因OsPOX1、抗壞血酸過氧化物酶基因OsAPXa、過氧化氫酶基因OsCATC等均受低溫誘導^［11-13］。

轉錄組測序技術是近年來在基因芯片、基因表達譜等技術上發展起來的新一代測序技術。與之前的測序方法相比，轉錄組測序通量高，能夠快速全面地對某一物種的特定組織或器官在某一特定狀態下進行全轉錄本測序分析。目前此項技術已經成功應用在多種作物中^［14-16］。在水稻低溫脅迫應答方面，亦有眾多相關報道^［17-18］。本研究針對黑龍江的生態類型和氣候條件，以當地孕穗期耐冷性差異較大的粳稻品種龍粳25（耐冷）、龍粳11（不耐冷）作為研究對象，通過在孕穗期低溫脅迫下的轉錄組測序分析，進一步闡明黑龍江水稻孕穗期耐冷性的遺傳調控分子機制，探索發掘寒地水稻孕穗期耐冷調控基因，以期為后續耐冷分子育種提供理論參考。

1 材料與方法

1.1 試驗材料

試驗材料為黑龍江省農業科學院水稻研究所育成并提供的龍粳25、龍粳11，二者均為黑龍江省第三積溫區適宜種植品種，孕穗期耐冷性差異明顯。

1.2 材料種植

本研究采用盆栽試驗。試驗地點為黑龍江省農業科學院水稻研究所（下簡稱水稻所）田間試驗場。盆栽所需土壤均采自水稻所田間。土壤進行干燥粉碎并過篩混勻后，裝入直徑34 cm、高35 cm的塑料桶中，每桶裝干土量約為12 kg。試驗材料于2018年4月20日播種，待生長到3葉1心時（5月15日），挑選長勢一致的幼苗移栽至塑料桶中，每桶均勻栽植15株。每品種分別種植12桶，各取6桶作為對照。生長過程中，為保證各單株生育進程一致，所有材料均不留分蘗，只留主莖。孕穗期判定依據Sato等的方法^［12］，略作改動。當劍葉與倒2葉的葉枕距達-4～0 cm時，即認為其進入孕穗期，并對該單株掛牌標記。將同時進入孕穗期的試驗材料搬至光照培養箱（HPG-280，哈爾濱市東聯電子技術開發有限公司生產）中，12 ℃低溫持續處理 2 d，相對濕度為75%，光—暗周期設定為8 h—16 h。低溫處理結束后，將試驗材料移至正常環境生長至成熟。收獲時計算結實率。

1.3 材料取樣及RNA抽提

取水稻穎花作為試驗材料進行轉錄組測序分析。分別在未處理和低溫處理2 d的龍粳25、龍粳11中選取生育進程一致的水稻主莖，剝出幼穗，取幼穗上1/3的穎花，從中挑取長度為3～5 mm 的穎花，稱重約0.5 g，錫箔紙包裹，過液氮迅速冷凍，于-80 ℃保存，每處理重復3次。RNA抽提采用Trizol試劑盒（Thermo Fisher Scientific公司生產）的，方法參照試劑盒提取說明，RNA的純凈和完整程度通過瓊脂糖凝膠電泳檢測。

1.4 文庫構建及測序

文庫構建及轉錄組測序委托北京百邁克生物技術有限公司（www.biomarker.com.cn）進行。文庫構建完成后，對其質量進行檢測，檢測合格后的文庫按照目標下機數據量進行pooling，用Illumina HiSeq 2500平臺根據操作標準進行測序，產生 150 bp 的雙端reads。

1.5 測序數據過濾及組裝

測序后下機產生的原始數據（Raw Data）通常以FASTQ格式提供。過濾掉包含接頭的低質量數據后即為可用數據（clean reads），通過Q20、Q30、GC 含量以及序列重復水平等參數來評估所得可用數據能否用于進一步分析。以Oryza sativa L. ssp. japonica （Oryza sativa IRGSP-1.0）作為本試驗水稻的參考基因組，通過TopHat2將可用數據比對到參考基因組上。以FPKM 值判定轉錄本或基因表達水平。

1.6 差異表達基因的篩選及生物信息學分析

通過DEGseq2軟件分析基因表達水平。參數設置具體為|log₂ Fold Change （FC）|≥1，1 discovery rate （FDR）≤0.01。GO功能分類中基因以及基因數目的篩選是基于GO數據庫（http：//www.geneontology.org/）的映射完成。GO功能富集是通過GOSeq R軟件包，以Wallenius 非中心超幾何分布檢驗進行分析。KEGG Pathway 的富集分析是通過KOBAS軟件（http：//kobas.cbi.pku.edu.cn/home.do）進行。GO及KEGG的富集標準皆為P≤0.01。

1.7 差異表達基因的qRT-PCR分析

選用Thermo Fisher Scientific公司的Trizol試劑盒，分別對龍粳25、龍粳11在處理前（0 d）和處理后（2 d）的穎花進行抽提RNA，穎花選取部位及方法同“1.3”節。純化后的RNA通過 HiScript Q RT SuperMix for qPCR試劑盒反轉錄合成 cDNA。選取表達差異大的基因，通過primer3軟件設計qRT-PCR引物（表1），以actin為內參。qRT-PCR反應在ABI 7500 qPCR thermal cycler （Applied Biosystems Inc.，Carlsbad，CA，USA）儀器上運行，重復3次，通過2^-ΔΔC_T法計算基因的相對表達量變化。

1.8 數據處理

用WPS Office分析數據，用SPSS 22.0進行顯著性檢驗（α=0.05），用GraphPad Prism 9.0作圖。

2 結果與分析

2.1 孕穗期低溫下結實率變化

孕穗期低溫直接影響水稻的結實率。通過對龍粳25、龍粳11在孕穗期進行12 ℃低溫處理2 d，由圖1可以看出，未經過低溫處理時，龍粳25、龍粳11的結實率分別為94.5%、93.05%，二者間未有明顯差異。當低溫處理2 d后，龍粳25的結實率為89.27% 較對照未有明顯降低。而龍粳11在低溫處理2 d后，結實率降為37.85%，較對照呈極顯著降低。可以看出，龍粳25、龍粳11在正常環境中結實率沒有明顯差異，低溫是造成二者結實率差異的原因，并且龍粳25的孕穗期耐冷性要明顯強于龍粳11。

2.2 測序結果組裝分析

本研究通過對龍粳25、龍粳11低溫處理2 d的穎花進行轉錄組測序分析，共構建12個文庫，去掉raw reads，共得到689.23 Mb的clean reads，GC含量為54.52%～56.58%，Q30均高于91.10%，平均81.70%的clean reads比對到參考基因組上（表2）。以上結果表明，龍粳25、龍粳11的轉錄組測序質量較高，可以繼續進行后續分析。

2.3 差異表達基因識別與分析

低溫處理下，基因的表達量會根據自身對低溫的敏感程度發生相應變化，其差異表達程度可反映對低溫的敏感度。因此，在數據質控分析之后，進一步對差異表達基因（differentially expressed gene，DEG）進行統計分析。由圖2-a可知，以|log₂ FC|≥1，FDR≤0.01為篩選閾值，在龍粳11中，低溫處理2 d后共有9 468個DEG，其中包括上調表達基因 4 915 個、下調表達基因4 553個。而在龍粳25中，低溫處理2 d后共有6 784個基因的表達量發生變化，其中包括上調表達基因3 499個、下調表達基因 3 285 個（圖2-b）。進一步通過韋恩圖（圖2-c）分析可以看出，龍粳25、龍粳11共有的DEG數目為3 508個，而龍粳25、龍粳11中特有的DEG分別為3 276、5 960個。結果表明，低溫處理2 d后，龍粳11的DEG數目明顯多于龍粳25，表明龍粳11對低溫更為敏感，低溫反應更加明顯。

2.4 差異表達基因的GO富集分析

通過分析DEG的韋恩圖，發現有3 508個DEG是龍粳25、龍粳11在低溫處理后共有的，而在龍粳25中有3 276個DEG是其低溫下特有的。為近一步明確其功能分類，對這些基因分別展開GO富集分析。通過對二者共有的DEG進行GO分析，以KS≤0.01為篩選標準，生物進程、細胞組分、分子功能的GO功能類別明顯富集到的類別分為95、20、41條。圖3-a分別列出每個類別中含有基因數最多的前20條terms。生物進程方面，主要富集在生物進程調控（GO：0050789，631個DEG）、氧化還原反應（GO：0055114，301個DEG）、非生物刺激應答（GO：0009628，302個DEG）、蛋白磷酸化（GO：0006468，279個DEG）等類別中。細胞組分方面，主要富集在膜（GO：0016020，967個DEG）、細胞質膜包圍的囊泡（GO：0016023，773個DEG）、原生質膜（GO：0005886，427個DEG）、細胞核（GO：0005634，382個DEG）等類別中。在分子功能方面，主要富集在ATP結合功能（GO：0005524，439個DEG）、絲氨酸/蘇氨酸蛋白激酶活性（GO：0004674，211個DEG）、轉運活性（GO：0005215，189個DEG）、序列特異性DNA結合轉錄（GO：0003700，156個DEG）等類別中。

由于龍粳25具有孕穗期強耐冷性，為了發現其低溫下特有的功能類別，對其特有的DEG進行GO富集分析。同樣以KS≤0.01為篩選標準，明顯富集到的生物進程、細胞組分、分子功能的GO功能類別分別為67、24、32條。由圖3-b可見，在生物進程方面，龍粳25特有基因多富集在生物進程調控（GO：0019222，226個DEG）、非生物刺激應答（GO：0009628，222個DEG）、氧化還原反應（GO：0055114，215個DEG）、碳水化合物代謝過程（GO：0005975，176個DEG）等類別中。在細胞組分方面，龍粳25特有基因富集最多的前4條類別與龍粳25、龍粳11共有基因富集的類別相同，而所包含的基因數則分別為723、613、324、289個。在分子功能方面，龍粳25特有基因主要富集在絲氨酸/蘇氨酸蛋白激酶（GO：0004674，158個DEG）、序列特異性DNA結合轉錄（GO：0003700，90個DEG）、轉移酶活性/轉移糖基（GO：0016757，81個DEG）、血紅素結合部位（GO：0020037，69個DEG）等類別中。

2.5 差異表達基因的KEGG富集分析

為近一步探明低溫對寒地水稻代謝途徑的影響，更加明確地闡釋DEG的生物學功能，對龍粳25、龍粳11低溫處理下共有的DEG以及龍粳25特有的DEG分別進行KEGG富集分析。結果表明，二者共有的3 508個DEG明顯富集到112個KEGG通路中。通過列出富集程度最高的前20個通路（圖4-a）可以看出，DEG主要富集在植物激素信號轉導（ko04075，46個DEG，富集因子1.81）、碳代謝（ko01200，44個DEG，富集因子1.61）、淀粉蔗糖代謝（ko00500，43個DEG，富集因子1.86）、苯丙烷合成（ko00940，35個DEG，富集因子1.61）、糖酵解/糖原異生通路（ko00010，33個DEG，富集因子2.67）等途徑。而在龍粳25特有的DEG中，共富集到110個KEGG通路中。同樣列出富集程度最高的前20個通路，DEG多富集在植物激素與信號轉導（ko04075，41個DEG，富集因子1.94）、內質網蛋白加工（ko04141，34個DEG，富集因子1.61）、苯丙烷合成（ko00940，32個DEG，富集因子1.76）、氨基糖和核酸糖代謝（ko00520，21個DEG，富集因子1.87）、谷胱甘肽代謝（ko00480，17個DEG，富集因子2.10）等通路中（圖4-b）。綜合比較2個部分KEGG通路，發現植物激素信號轉導（ko04075）在龍粳25、龍粳11共有DEG及龍粳25特有DEG中都被明顯富集到含DEG均為最多的通路，表明這條通路在寒地水稻低溫脅迫應答中起著重要的作用。

2.6 差異表達基因的qRT-PCR驗證

為驗證轉錄組測序結果的準確性，從測序結果中隨機選取6個表達差異大的基因進行qRT-PCR分析。分別以龍粳11未做低溫處理時的各個基因表達量作為對照（相對表達量為1），qRT-PCR結果表明，經過低溫處理后2個品種的6個基因表達均產生差異變化；與轉錄組測序結果相比，雖然基因表達量的變幅與轉錄組結果不盡一致，但變化趨勢與轉錄測序結果均一致，證實了轉錄組測序結果的準確性及可靠性（圖5）。

3 結論與討論

本研究通過RNA-seq技術對寒地水稻低溫下展開轉錄組測序，比較2個品種低溫脅迫后的差異表達基因。研究結果表明，低溫脅迫后，龍粳11的DEG為9 468個，龍粳25為6 784個；低溫處理后，龍粳11的DEG明顯多于龍粳25，表明龍粳11對低溫更加敏感，而龍粳25對低溫則表現出更強的耐受性，這也可以從二者在低溫下的結實率得到驗證。GO富集分析發現，低溫脅迫下，生物進程調控、氧化還原反應、非生物刺激應答、蛋白磷酸化等GO類別的反應比較大。KEGG富集分析得出，檸檬烯和蒎烯的降解、糖酵解/糖異生、半乳糖代謝、（角質、木質、蠟質）生物合成、植物激素信號轉導等途徑在低溫脅迫后富集的程度最高。以上分析表明，低溫對寒地水稻的影響涉及到生長發育的各個方面，包括各類復雜的新陳代謝變化。

本研究所用的材料為孕穗期耐冷的龍粳25和不耐冷的龍粳11，在探究寒地水稻特有的低溫應答防御機制及相關基因時，龍粳25特有的DEG及其GO類別值得進一步深入研究。本研究中，龍粳25富集到67條生物進程、24條細胞組分、32條分子功能類別。由于生物進程能夠直接反映生物體生長發育情況，對此類別展開分析。氧化還原反應進程不但在二者共有的DEG中得到明顯富集，并且在龍粳25特有的DEG中也被明顯富集，表明此GO類別在寒地水稻低溫應答中起重要作用。

植物體內源激素雖然含量很低，但是在生物脅迫、非生物脅迫等眾多脅迫因子方面起到積極的防衛作用^［19］。目前研究較多的激素包括脫落酸（ABA）、赤霉素（GA）、生長素（IAA）、乙烯（ETH）、水楊酸（SA）、茉莉酸（JA）、細胞分裂素（CK）、油菜素內酯（BR）等^［20］。而IAA、ABA、ETH等激素在水稻低溫應答中起積極作用^［21-22］。本研究中，植物激素信號轉導通路在龍粳11、龍粳25共有的DEG以及龍粳25特有的DEG都被明顯富集，且含有基因數也是最多的，表明此通路在寒地水稻低溫應答中同樣起關鍵調控作用，這一結果也與前人研究結果^［23］一致。

低溫脅迫會對水稻植株的正常生長產生各種生理及分子調控機制方面的影響；水稻通過自身調整，盡可能適應溫度脅迫，以保持其正常生長^［24-25］。黑龍江地處寒稻作區，由于其特殊的地理環境及生態類型，低溫是寒地水稻正常生長的重要限制因子^［26］。本研究以寒地水稻孕穗期耐冷差異明顯的品種龍粳25、龍粳11為研究材料，在前人對其結實率研究的基礎上，對低溫脅迫下穎花進行轉錄組測序分析。通過差異表達基因分析，發現龍粳11的DEG數量明顯多于龍粳25，表明龍粳11對低溫反應更為敏感。GO、KEGG分析表明，氧化還原反應進程、低溫應答、植物激素信號轉導等眾多GO類別及KEGG通路均參與到低溫脅迫應答中。此外，通過qRT-PCR驗證，驗證了RNA-seq結果的準確性，期待本研究結果對寒地水稻孕穗期耐冷的分子遺傳調控機理的解析有一定的推動作用。

參考文獻：

［1］程式華，李 建. 現代中國水稻［M］. 北京：金盾出版社，2007.

［2］Howarth C J，Ougham H J. Gene expression under temperature stress［J］. The New Phytologist，1993，125（1）：1-26.

［3］Guo Z H，Cai L J，Chen Z Q，et al. Identification of candidate genes controlling chilling tolerance of rice in the cold region at the booting stage by BSA-Seq and RNA-Seq［J］. Royal Society Open Science，2020，7（11）：201081.

［4］Guo Z H，Liu C X，Xiao W M，et al. Comparative transcriptome profile analysis of anther development in reproductive stage of rice in cold region under cold stress［J］. Plant Molecular Biology Reporter，2019，37（3）：129-145.

［5］Mamun E A，Alfred S，Cantrill L C，et al. Effects of chilling on male gametophyte development in rice［J］. Cell Biology International，2006，30（7）：583-591.

［6］Oda S，Kaneko F，Yano K，et al. Morphological and gene expression analysis under cool temperature conditions in rice anther development［J］. Genes amp; Genetic Systems，2010，85（2）：107-120.

［7］Qiu C P，Ethier G，Pepin S，et al. Persistent negative temperature response of mesophyll conductance in red raspberry （Rubus idaeus L.） leaves under both high and low vapour pressure deficits：a role for abscisic acid？［J］. Plant，Cell amp; Environment，2018，41（4）：876.

［8］李健陵，霍治國，吳麗姬，等. 孕穗期低溫對水稻產量的影響及其生理機制［J］. 中國水稻科學，2014，28（3）：277-288.

［9］Saito K，Hayano-Saito Y，Kuroki M，et al. Map-based cloning of the rice cold tolerance gene Ctb1［J］. Plant Science，2010，179（1/2）：97-102.

［10］Vannini C，Campa M，Iriti M，et al. Evaluation of transgenic tomato plants ectopically expressing the rice Osmyb4 gene［J］. Plant Science，2007，173（2）：231-239.

［11］Kim S H，Choi H S，Cho Y C，et al. Cold-responsive regulation of a flower-preferential class Ⅲ peroxidase gene，OsPOX1，in rice （Oryza sativa L.）［J］. Journal of Plant Biology，2012，55（2）：123-131.

［12］Sato Y，Masuta Y，Saito K，et al. Enhanced chilling tolerance at the booting stage in rice by transgenic overexpression of the ascorbate peroxidase gene，OsAPXa［J］. Plant Cell Reports，2011，30（3）：399-406.

［13］Guo Z H，Cai L J，Liu C X，et al. Global analysis of differentially expressed genes between two Japonica rice varieties induced by low temperature during the booting stage by RNA-Seq［J］. Royal Society Open Science，2020，7（6）：192243.

［14］韋秀葉，趙信林，郭 媛，等. 麻類作物轉錄組測序分析研究進展［J］. 中國麻業科學，2020，42（3）：128-134.

［15］吳 倩，張 磊，黃志平，等. 轉錄組測序及其在野生大豆基因資源發掘中的應用［J］. 大豆科學，2013，32（6）：845-851.

［16］許 波，張偉強，馮曉曦，等. 轉錄組測序技術在玉米中應用研究進展［J］. 玉米科學，2014，22（1）：67-72，78.

［17］Bai B，Wu J，Sheng W T，et al. Comparative analysis of anther transcriptome profiles of two different rice male sterile lines genotypes under cold stress［J］. International Journal of Molecular Sciences，2015，16（5）：11398-11416.

［18］Zhang L，Li X Y，Li D K，et al. CARK1 mediates ABA signaling by phosphorylation of ABA receptors［J］. Cell Discovery，2018，4：30.

［19］叢 慶，張 琪，宋麗莉，等. 激素在植物冷脅迫應答中的角色［J］. 核農學報，2016，30（3）：614-619.

［20］張明菊，朱 莉，夏啟中. 植物激素對脅迫反應調控的研究進展［J］. 湖北大學學報（自然科學版），2021，43（3）：242-253，263.

［21］劉 芬. 低溫脅迫對細枝木麻黃無性系生理指標和轉錄組的影響［D］. 長沙：中南林業科技大學，2015：5-7.

［22］Jain M，Khurana J P. Transcript profiling reveals diverse roles of auxin-responsive genes during reproductive development and abiotic stress in rice［J］. The FEBS Journal，2009，276（11）：3148-3162.

［23］高紅秀，朱 琳，劉天奇，等. 水稻植物激素響應低溫脅迫反應的轉錄組分析［J］. 分子植物育種，2021，19（13）：4188-4197.

［24］胡瀟婕，毛東海. 基于RNA-Seq技術分析植物激素信號途徑在水稻幼苗中對低溫脅迫的應答規律［J］. 農業現代化研究，2019，40（5）：878-890.

［25］劉琳帥，卞景陽，孫興榮，等. 水稻低溫冷害的研究進展［J］. 江蘇農業科學，2022，50（24）：9-15.

［26］賀 梅，宋冬明. 寒地稻作區冷害防御措施［J］. 北方水稻，2015，45（1）：43-44.