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

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

早在2013年，本研究團(tuán)隊(duì)就對(duì)山欄稻展開了研究，劉欣欣等［9］克隆了海南山欄稻的HKT2基因片段，研究了海南山欄稻運(yùn)輸Na+和K+的分子機(jī)制；劉維俊等［4］研究了干旱脅迫下2個(gè)山欄稻與栽培水稻品種苗期表型性狀及生理差異；吳丹等［10］通過田間試驗(yàn)設(shè)計(jì)覆蓋秸稈旱作、水作和傳統(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ù)。我們?cè)趯?duì)山欄稻前期研究的基礎(chǔ)上，通過RNA-seq（RNA Sequencing）Illumina 二代高通量測(cè)序技術(shù)，對(duì)山欄稻抗旱相關(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)條件如下：相對(duì)濕度50%-70%；光照16 h，30℃/黑暗8 h，26℃培養(yǎng)。水培營(yíng)養(yǎng)液用Yoshida水稻營(yíng)養(yǎng)液配制，2 d更換一次。

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

1.2.3 數(shù)據(jù)過濾、重新序列組裝以及unigenes的功能注釋 將測(cè)序得到的原始序列Raw reads進(jìn)行過濾得到Clean reads。取大于200 bp且序列最長(zhǎng)的轉(zhuǎn)錄本作為unigene，將6個(gè)庫混合拼接成一個(gè)轉(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)行比對(duì)。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ù)測(cè)基因和表達(dá)水平分析 CDS（Coding sequence）預(yù)測(cè)是按照NR、Swissprot、KEGG GENES的優(yōu)先級(jí)順序，將unigene比對(duì)到蛋白數(shù)據(jù)庫。采用RSEM（RNA-Seq by Expectation-Maximization）軟件，對(duì)bowtie的比對(duì)結(jié)果進(jìn)行統(tǒng)計(jì)，得到每個(gè)樣品比對(duì)到每個(gè)基因上的readcount數(shù)目，并對(duì)其進(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á)水平計(jì)算的基礎(chǔ)上，我們采用DESeq進(jìn)行分析。為了篩選到更多的與干旱脅迫相關(guān)的候選基因，對(duì)干旱脅迫處理前后的閾值設(shè)置為P-value＜0.1，|log2FoldChange|＞1。

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

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

2 結(jié)果

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

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

表2 Reads與參考序列比對(duì)情況一覽表

2.2 SNP與InDel

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

2.3 基因表達(dá)對(duì)比

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 討論

本研究對(duì)山欄稻白沙糯進(jìn)行了干旱脅迫應(yīng)答的轉(zhuǎn)錄組分析發(fā)現(xiàn)，所有差異表達(dá)的基因都為單基因，有2 791個(gè)差異基因受干旱脅迫誘導(dǎo)，其中1 586個(gè)基因在脅迫后表達(dá)上調(diào)，1 065個(gè)基因在脅迫后表達(dá)下調(diào)。轉(zhuǎn)錄因子可以控制一系列下游脅迫反應(yīng)，從而啟動(dòng)信號(hào)傳導(dǎo)中的級(jí)聯(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個(gè)差異表達(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］對(duì)海南野生稻和黃立鈺［19］對(duì)水稻H471干旱處理后轉(zhuǎn)錄組測(cè)序的差異表達(dá)轉(zhuǎn)錄因子的類型相似，不同的是山欄稻干旱脅迫處理后差異表達(dá)的轉(zhuǎn)錄因子數(shù)量比前兩者的數(shù)量更多，這也從側(cè)面說明了山欄稻抗旱性更強(qiáng)。

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

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

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

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

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

4 結(jié)論

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