海南山欄稻抗旱基因轉(zhuǎn)錄組分析

翟楠鑫 遲會 夏玥琳 劉彩月 裴新梧 袁潛華

（1. 海南大學(xué)熱帶作物學(xué)院，海口 570228；2. 中國農(nóng)業(yè)科學(xué)院生物技術(shù)研究所，北京 100081）

山欄亦作“砍山”或“砍山欄”，山欄稻是海南特有的旱稻，主要分布于保亭、白沙和樂東等山坡地，盛產(chǎn)于五指山地區(qū)。不同于水稻的育秧插秧，山欄稻是點(diǎn)播直播，主要種植在山坡地，產(chǎn)量低，加之種植規(guī)模小，品種老化，導(dǎo)致山欄稻種資源流失嚴(yán)重［1］。山欄稻具有抗旱、耐熱、抗病、抗蟲等優(yōu)良特性，是珍貴的旱稻種質(zhì)資源［2］。

山欄稻的農(nóng)藝性狀一般表現(xiàn)為植株偏高，穗形披散，結(jié)實(shí)率較高，粒型多為大長粒，葉片寬大而披垂，冠層密封分蘗力弱，穗大而彎垂等特點(diǎn)［2］。劉維俊等［4］研究了干旱脅迫條件下山欄稻與栽培稻品種苗期表型性狀及生理差異，結(jié)果表明山欄稻苗期具有較強(qiáng)的抗旱性。唐力瓊等［5］采用PEG人工模擬干旱的方法，研究20個山欄稻資源在芽期和苗期的抗旱性，發(fā)現(xiàn)山欄稻品種總體抗旱性較強(qiáng)。徐建欣等［3］對22個海南山欄稻品種進(jìn)行了全生育期抗旱性鑒定，發(fā)現(xiàn)山欄稻品種不同生育期間所表現(xiàn)出的抗旱能力具有極顯著的正相關(guān)關(guān)系，這表明山欄稻可以在早期（芽期或苗期）快速篩選抗旱品種。在山欄稻品種選育方面，研究人員也成功培育出了“熱大99 W”［6］和“山欄陸1號”［7］等抗旱豐產(chǎn)品系。劉躍林等［8］通過基因芯片技術(shù)，研究海南山欄稻在干旱逆境條件下的基因表達(dá)譜，將其中的OsMS17在水稻中超表達(dá)，轉(zhuǎn)基因植株的苗期和成株期的抗旱性明顯增強(qiáng)。

早在2013年，本研究團(tuán)隊(duì)就對山欄稻展開了研究，劉欣欣等［9］克隆了海南山欄稻的HKT2基因片段，研究了海南山欄稻運(yùn)輸Na+和K+的分子機(jī)制；劉維俊等［4］研究了干旱脅迫下2個山欄稻與栽培水稻品種苗期表型性狀及生理差異；吳丹等［10］通過田間試驗(yàn)設(shè)計覆蓋秸稈旱作、水作和傳統(tǒng)旱作3種栽培方式，發(fā)現(xiàn)山欄稻在水田環(huán)境下具有減少干旱脅迫、優(yōu)化農(nóng)藝性狀和提高產(chǎn)量等諸多優(yōu)點(diǎn)并建立了山欄稻水田灌溉種植的高產(chǎn)栽培方法［11］，為山欄稻從傳統(tǒng)山區(qū)旱地種植引向普通旱作和水田種植提供科學(xué)的依據(jù)。我們在對山欄稻前期研究的基礎(chǔ)上，通過RNA-seq（RNA Sequencing）Illumina 二代高通量測序技術(shù)，對山欄稻抗旱相關(guān)基因轉(zhuǎn)錄組進(jìn)行表達(dá)分析及鑒定，以進(jìn)一步明確海南山欄稻在干旱逆境條件下的基因表達(dá)譜，為解析山欄稻的抗旱機(jī)理奠定基礎(chǔ)。

1 材料與方法

1.1 材料

山欄稻品種白沙糯，由海南大學(xué)熱帶作物學(xué)院生物多樣性與水稻種質(zhì)創(chuàng)新實(shí)驗(yàn)室收集保存，該品種遺傳穩(wěn)定，根系發(fā)達(dá)，耐旱性強(qiáng)。

1.2 方法

1.2.1 植物材料培育 山欄稻品種白沙糯在光照培養(yǎng)箱中水培，培養(yǎng)條件如下：相對濕度50%-70%；光照16 h，30℃/黑暗8 h，26℃培養(yǎng)。水培營養(yǎng)液用Yoshida水稻營養(yǎng)液配制，2 d更換一次。

1.2.2 建庫和轉(zhuǎn)錄組測序 山欄稻白沙糯培養(yǎng)至三葉一心期后，對照組繼續(xù)營養(yǎng)液培養(yǎng)，處理組加20%的PEG-6000（聚乙二醇6000）培養(yǎng)液模擬干旱脅迫處理。取0 h、24 h、48 h、72 h的水稻樣本分別作為不同脅迫時段材料并提取總RNA。RNA樣品檢測合格后，用3個生物學(xué)重復(fù)的對照組和處理組（脅迫72 h）構(gòu)建6個cDNA文庫，庫檢合格后用Illumina HiSeq 2500進(jìn)行測序。

1.2.3 數(shù)據(jù)過濾、重新序列組裝以及unigenes的功能注釋 將測序得到的原始序列Raw reads進(jìn)行過濾得到Clean reads。取大于200 bp且序列最長的轉(zhuǎn)錄本作為unigene，將6個庫混合拼接成一個轉(zhuǎn)錄本。序列拼接后與NR（Non-redundant）、PFAM（Protein Families）、Swiss-Prot（http：//www.expasy.ch/sprot/）和KEGG（Kyoto Encyclopedia of Genes and Genomes）等數(shù)據(jù)庫進(jìn)行比對。GO（Gene Ontology）功能注釋基于NR和PFAM兩部分的蛋白注釋結(jié)果，軟件為Blast2GO v2.5和WEGO。并根據(jù)它們參與的KEGG代謝通路進(jìn)行分析，KEGG相關(guān)注釋通過KEGG Automatic Annotation Server完成。

1.2.4 CDS預(yù)測基因和表達(dá)水平分析 CDS（Coding sequence）預(yù)測是按照NR、Swissprot、KEGG GENES的優(yōu)先級順序，將unigene比對到蛋白數(shù)據(jù)庫。采用RSEM（RNA-Seq by Expectation-Maximization）軟件，對bowtie的比對結(jié)果進(jìn)行統(tǒng)計，得到每個樣品比對到每個基因上的readcount數(shù)目，并對其進(jìn)行FPKM（Fragments Per Kilobase of exon model per Million mapped reads）轉(zhuǎn)換，進(jìn)而分析基因的表達(dá)水平。

1.2.5 差異表達(dá)unigenes的驗(yàn)證和表達(dá)分析 在基因表達(dá)水平計算的基礎(chǔ)上，我們采用DESeq進(jìn)行分析。為了篩選到更多的與干旱脅迫相關(guān)的候選基因，對干旱脅迫處理前后的閾值設(shè)置為P-value＜0.1，|log2FoldChange|＞1。

為保證轉(zhuǎn)錄組測序的差異表達(dá)基因表達(dá)趨勢可靠性，我們從測序獲得基因中隨機(jī)選取差異表達(dá)基因進(jìn)行qRT-PCR驗(yàn)證其是否和轉(zhuǎn)錄組庫趨勢一致；另外選取了已報道干旱脅迫相關(guān)基因進(jìn)行驗(yàn)證，鑒定它們與干旱脅迫的相關(guān)性。qRT-PCR采用SYBR Green結(jié)合相對定量法在ABI 7500 Real Time System上進(jìn)行。選用OsActin作為內(nèi)參，所有樣本作3個重復(fù)。序列見表1。

表1 qRT-PCR分析特異性引物

2 結(jié)果

2.1 轉(zhuǎn)錄組拼接質(zhì)量分析及參考序列比對

將Illumina測序結(jié)果去除銜接子和質(zhì)量差的讀數(shù)后，把Clean reads映射到參考基因組數(shù)據(jù)中，對照組的3個重復(fù)基因序列個數(shù)分別為：38 504 957、35 184 648和26 794 172個；脅迫組的3個重復(fù)基因序列個數(shù)分別為：32 620 538、29 253 860和30 195 571個。這6個庫的測序錯誤率均小于0.05%，Q20均大于95.00%，GC含量也在50%左右，表明這6個庫的測序結(jié)果可靠，能夠用于后續(xù)分析。將Clean reads與參考序列進(jìn)行比對（表2），比對到參考序列上的Reads數(shù)量占總數(shù)的92%左右，其中比對到多個位置的Reads占總數(shù)的2.6%以下。

表2 Reads與參考序列比對情況一覽表

2.2 SNP與InDel

SNP（Single Nucleotide Polymorphisms），指 基因組水平上由單個核苷酸變異引起的DNA序列多態(tài)性。SNP主要為轉(zhuǎn)換和顛換兩種類型，一般二者的比例為2∶1。因?yàn)槠鋸V泛存在于基因組上，現(xiàn)在成為了遺傳研究的主要標(biāo)記。InDel（Insertion and Deletion）指基因組上小片段的插入和缺失，編碼區(qū)的InDel往往會引起基因功能的重大改變。我們首先通過samtools和picard-tools等工具對比對結(jié)果進(jìn)行染色體坐標(biāo)排序，然后通過變異檢測軟件GATK分別進(jìn)行SNP Calling和InDel Calling，最后進(jìn)行過濾，得到可靠的SNP和InDel檢測結(jié)果（圖1）。

2.3 基因表達(dá)對比

使用DESeq對得到的數(shù)據(jù)進(jìn)行基因差異表達(dá)分析，結(jié)果（表3）表明：山欄稻白沙糯干旱脅迫前后有2 791個差異表達(dá)的基因，其中1 669個表達(dá)上調(diào)，1 122個表達(dá)下調(diào)。與目前已知的水稻基因組數(shù)據(jù)對比，其中已知基因共2 650個，未知基因共141個。

圖1 SNP和InDels的統(tǒng)計數(shù)據(jù)

表3 干旱脅迫前后差異表達(dá)基因列表

2.4 差異基因GO富集分析和KEGG富集分析

為研究山欄稻基因調(diào)控與干旱脅迫的關(guān)系，在NR功能注釋基礎(chǔ)上，我們對脅迫前后差異表達(dá)基因進(jìn)行了GO分析。干旱脅迫前后差異表達(dá)基因共2 791個，有GO注釋的差異基因有1 599個，被分配到了1 524個GO terms功能組，所占比例最大的為生物過程和分子功能。我們統(tǒng)計出這1 524個組中有25個GO terms與脅迫響應(yīng)相關(guān)，其中對應(yīng)激的反應(yīng)所占比例最大（基因個數(shù)與所占比例分別為：164，10.26%，下同），其次從大到小分別為：對壓力的響應(yīng)（84，5.25%）、細(xì)胞對應(yīng)激的響應(yīng)（70，4.38%）、對化學(xué)應(yīng)激響應(yīng)（57，3.56%）、對非生物應(yīng)激的反應(yīng)（30，1.88%）、氧化脅迫響應(yīng)（23，1.44%）和對水的反應(yīng)（20，1.25%）。

如圖2-A所示，縱坐標(biāo)為富集的GO term，橫坐標(biāo)為該term中差異表達(dá)基因個數(shù)。帶“*”為顯著富集的GO term。在1 524個GO terms中達(dá)到顯著富集有24個，可大致分為以下幾類：蛋白質(zhì)磷酸化、蛋白激酶活性、磷酸化與磷代謝過程、修飾過程、催化活性、對非生物刺激的反應(yīng)、轉(zhuǎn)移酶活性、離子結(jié)合、陽離子結(jié)合、幾丁質(zhì)分解代謝過程、水解酶活性、氨基聚糖分解代謝過程和對水的反應(yīng)。

KEGG是系統(tǒng)分析基因功能、基因組信息的數(shù)據(jù)庫，它有助于研究者把基因及表達(dá)信息作為一個整體網(wǎng)絡(luò)進(jìn)行研究。在生物體內(nèi)，不同基因相互協(xié)調(diào)行使其生物學(xué)功能，通過Pathway顯著性富集分析，找出與整個轉(zhuǎn)錄組背景相比顯著性富集的Pathway。如（圖2-B）縱軸表示Pathway名稱，橫軸表示Rich factor，點(diǎn)的大小表示此Pathway中差異表達(dá)基因個數(shù)多少，而點(diǎn)的顏色對應(yīng)于不同的q-value范圍，q-value的值越接近于零，說明富集程度越顯著。可以看出干旱脅迫前后差異表達(dá)基因KEGG顯著富集在植物激素信號轉(zhuǎn)導(dǎo)，次生代謝物的生物合成，角質(zhì)、木栓質(zhì)和蠟生物合成中。

2.5 差異基因轉(zhuǎn)錄因子對比

轉(zhuǎn)錄因子（Transcription factor，簡稱TF）可以通過結(jié)合靶基因啟動子的順式作用元件，從而調(diào)控基因的表達(dá)。轉(zhuǎn)錄因子可以控制一系列下游脅迫反應(yīng)，從而啟動信號傳導(dǎo)中的級聯(lián)反應(yīng)，使細(xì)胞產(chǎn)生相應(yīng)的抗逆性，分析這些轉(zhuǎn)錄因子可以為水稻干旱脅迫下復(fù)雜的調(diào)控網(wǎng)絡(luò)提供重要的信息。對差異基因進(jìn)行篩選，發(fā)現(xiàn)可能與抗旱有關(guān)的轉(zhuǎn)錄因子，如表4所示為差異表達(dá)轉(zhuǎn)錄因子基因信息（部分），將庫中與抗旱有關(guān)的轉(zhuǎn)錄因子進(jìn)行分類統(tǒng)計得到如下結(jié)果：Znf（Zinc finger）家族轉(zhuǎn)錄因子在干旱脅迫后上調(diào)的數(shù)量最多，達(dá)41個，占相關(guān)上調(diào)轉(zhuǎn)錄因子總數(shù)的35.65%（表5）。

2.6 差異表達(dá)基因驗(yàn)證

為了驗(yàn)證測序結(jié)果的可靠性，我們從干旱脅迫相關(guān)的已知基因中隨機(jī)選取6個基因并利用qRTPCR進(jìn)行驗(yàn)證，這6個基因的相關(guān)注釋信息以及在庫中的表達(dá)量如表6所示；另外選取已報道的與干旱脅迫相關(guān)的基因RAB21和OsLEA3-2進(jìn)行試驗(yàn)，驗(yàn)證干旱脅迫處理是否適當(dāng)。圖3是8個基因的qRT-PCR結(jié)果：RAB21和OsLEA3-2在脅迫后表達(dá)量顯著升高，表明模擬干旱脅迫處理適當(dāng)。基因Os02g0693700、Os07g0658600和Os11g0115100為已注釋的上調(diào)基因；基因Os09g0351700和Os07g025-1900為已注釋下調(diào)基因；基因Os01g0910800為未注釋表達(dá)上調(diào)基因，且上調(diào)倍數(shù)明顯。這6個基因與測序的表達(dá)趨勢基本一致，說明本次研究的RNAseq高通量測序結(jié)果可靠。

圖2 差異基因GO富集柱狀圖（A）及差異基因KEGG富集散點(diǎn)圖（B）

2.7 干旱脅迫相關(guān)差異基因不同脅迫時段的表達(dá)檢測

驗(yàn)證測序可靠后，運(yùn)用qRT-PCR技術(shù)，我們從干旱脅迫相關(guān)的候選基因中選取了4個基因來檢測這些基因在山欄稻白沙糯不同脅迫時期（0 h、24 h、48 h和72 h）的表達(dá)情況，這4個 基因的相關(guān)注釋信息以及在庫中的表達(dá)量如表6所示。其中基因Os01g0910800為未注釋上調(diào)基因；基因Os02g0693700和基因Os07g0658600為上調(diào)已注釋基因；基因Os07g0251900為下調(diào)已注釋基因。從圖4的結(jié)果可以看出：基因Os07g0658600、Os02g0693700和Os01g0910800隨著脅迫時間的延長，表達(dá)逐漸上升與干旱脅迫正相關(guān)；而基因Os07g0251900隨著脅迫時間的延長，表達(dá)逐漸下降與干旱脅迫負(fù)相關(guān)。

表4 差異表達(dá)轉(zhuǎn)錄因子基因信息（部分）

表5 不同庫中差異表達(dá)轉(zhuǎn)錄因子分布

3 討論

本研究對山欄稻白沙糯進(jìn)行了干旱脅迫應(yīng)答的轉(zhuǎn)錄組分析發(fā)現(xiàn)，所有差異表達(dá)的基因都為單基因，有2 791個差異基因受干旱脅迫誘導(dǎo)，其中1 586個基因在脅迫后表達(dá)上調(diào)，1 065個基因在脅迫后表達(dá)下調(diào)。轉(zhuǎn)錄因子可以控制一系列下游脅迫反應(yīng)，從而啟動信號傳導(dǎo)中的級聯(lián)反應(yīng)并使細(xì)胞產(chǎn)生相應(yīng)的抗逆性。目前與水稻干旱脅迫相關(guān)的轉(zhuǎn)錄因子主要有：bZIP、NAC、DREB/CBFs、Zinc finger蛋白、HSP和MYB等轉(zhuǎn)錄因子［12-17］。在本研究中獲得的山欄稻抗旱轉(zhuǎn)錄組數(shù)據(jù)中，有184個差異表達(dá)轉(zhuǎn)錄因子，在干旱脅迫后上調(diào)的轉(zhuǎn)錄因子主要有：Znf、AP2/ERF、MYB、NAC、bZIP等家族轉(zhuǎn)錄的因子，其中Znf家族轉(zhuǎn)錄因子在脅迫后上調(diào)的數(shù)量最多，占相關(guān)上調(diào)轉(zhuǎn)錄因子總數(shù)的35.65%。這些上調(diào)的轉(zhuǎn)錄因子與田新杰［18］對海南野生稻和黃立鈺［19］對水稻H471干旱處理后轉(zhuǎn)錄組測序的差異表達(dá)轉(zhuǎn)錄因子的類型相似，不同的是山欄稻干旱脅迫處理后差異表達(dá)的轉(zhuǎn)錄因子數(shù)量比前兩者的數(shù)量更多，這也從側(cè)面說明了山欄稻抗旱性更強(qiáng)。

山欄稻白沙糯差異基因的GO分類分析顯示，1 599個單基因在功能上用GO術(shù)語注釋，分屬1 524個GO terms功能組。山欄稻白沙糯顯著富集在25個terms，包括蛋白質(zhì)磷酸化、非生物應(yīng)激反應(yīng)、幾丁質(zhì)分解代謝、蛋白激酶活性、金屬離子結(jié)合和幾丁質(zhì)酶活性等terms。蛋白激酶主要起著催化蛋白質(zhì)磷酸化的作用，而蛋白質(zhì)的可逆磷酸化是細(xì)胞信號識別與轉(zhuǎn)導(dǎo)的重要環(huán)節(jié)，與植物干旱、高鹽、低溫、激素等反應(yīng)的信號傳遞有關(guān)［20］；大多數(shù)幾丁質(zhì)酶也可受某些非生物因素的誘導(dǎo)，其誘導(dǎo)表達(dá)是植物對逆境的一般反應(yīng)［21］；金屬離子如K+、Ca2+和Mg2+與作物的耐旱性極為密切，它們是維持作物細(xì)胞滲透壓的最主要離子［22］。這些研究結(jié)果表明上述通路的確與抗旱有關(guān)，與本研究結(jié)果一致。

表6 干旱脅迫前后差異表達(dá)基因信息

圖4 干旱脅迫相關(guān)基因不同脅迫時段的表達(dá)檢測

KEGG通路分析結(jié)果顯示，山欄稻白沙糯差異基因的富集在109個通路中。在植物激素信號轉(zhuǎn)導(dǎo)，次生代謝物的生物合成，角質(zhì)、木栓質(zhì)和蠟生物合成這3個途徑中顯著富集。植物激素是在植物體內(nèi)合成，對生長發(fā)育產(chǎn)生顯著作用的微量有機(jī)物質(zhì)。前人研究表明，干旱條件下各種植物激素調(diào)節(jié)通過相互影響協(xié)同作用來抵御干旱［23］。角質(zhì)層覆蓋在植物最表層，是所有陸生植物最重要的保護(hù)層之一。角質(zhì)層包括角質(zhì)和蠟質(zhì)，在植物生長發(fā)育過程中可以保護(hù)植物免受環(huán)境壓力，使其地上組織適應(yīng)水分脅迫的環(huán)境及其它非生物和生物逆境，我們研究的結(jié)果也很好的印證了這一點(diǎn)。在水分脅迫下，植物組織中次生代謝產(chǎn)物的濃度常常上升，包括氰甙、其他硫化物、萜類化合物、生物堿、單寧和有機(jī)酸等［24］。對比上述研究及本試驗(yàn)結(jié)果，我們認(rèn)為山欄稻白沙糯主要通過滲透調(diào)節(jié)物質(zhì)（金屬離子）、角質(zhì)蠟質(zhì)保護(hù)和激素調(diào)節(jié)來抵御干旱脅迫。

本研究首次利用高通量測序技術(shù)對海南旱稻山欄稻進(jìn)行轉(zhuǎn)錄組分析及表達(dá)鑒定，我們發(fā)現(xiàn)有1 599個單基因在功能上有注釋，有1 192個差異基因是未注釋基因。其中以基因Os01g0910800為代表的若干未注釋基因在干旱脅迫誘導(dǎo)后上調(diào)，且倍數(shù)非常高，說明這些基因可能受干旱脅迫誘導(dǎo)表達(dá)，與水分脅迫相關(guān)，這些未注釋基因的功能有待我們?nèi)パ芯亢吞剿鳌Ｎ覀儷@得了海南山欄稻在干旱脅迫條件下的基因表達(dá)譜，期望為旱稻的進(jìn)一步深入研究及利用提供理論依據(jù)。

4 結(jié)論

山欄稻白沙糯干旱脅迫后有2 791個基因差異表達(dá)，未注釋基因1 192個。其中1 669個基因表達(dá)上調(diào)，1 122個基因表達(dá)下調(diào)；對差異基因中的轉(zhuǎn)錄因子進(jìn)行分析，發(fā)現(xiàn)山欄稻白沙糯中有184個差異表達(dá)的轉(zhuǎn)錄因子，其中Znf家族轉(zhuǎn)錄因子是在干旱脅迫后上調(diào)的數(shù)量最多的，共41個。本研究發(fā)現(xiàn)山欄稻白沙糯通過植物激素、角質(zhì)阻力增大和滲透壓調(diào)節(jié)來提高抗旱性。一系列未注釋基因（如基因Os01g0910800）隨著脅迫時間的延長，表達(dá)逐漸上升且增長倍數(shù)明顯。