水稻干旱脅迫的small RNA轉錄組分析

劉源, 張秀妍, 徐妙云, 鄭紅艷, 鄒俊杰, 張蘭, 王磊

(中國農業科學院生物技術研究所， 北京 100081)

干旱脅迫影響植物生長發育、生理代謝、地理分布、品質產量，是威脅糧食安全的主要環境因素之一[1]。在干旱條件下，植物會產生一系列結構和生理改變，如葉片和根系形態改變、氣孔運動、滲透損傷、物質吸收失衡、光合能力降低等[2-5]。植物在干旱響應過程中，大量抗旱相關調控基因的表達水平會發生改變，包括轉錄因子基因、滲透調節基因、抗氧化代謝基因和逆境誘導基因等，其中，許多基因在轉錄和轉錄后水平受small RNA調控[2, 5]。因此，small RNA及其靶基因組成調控模塊是植物干旱脅迫響應過程中的關鍵因子。

microRNA(miRNA)是廣泛分布于植物基因組中的內源性單鏈非編碼small RNA，長度18～24個核苷酸(nt)。miRNA通過結合靶基因上的互補序列誘導轉錄水平基因沉默(transcriptional gene silencing，TGS)和轉錄后水平基因沉默(post-transcriptional gene silencing，PTGS)而發揮功能[6]。miRNA不僅是植物生長發育和生理代謝活動的重要調節因子，同時也參與了各種生物和環境脅迫(干旱、高鹽、高溫、冷凍等)的響應過程[6-7]。目前，許多植物物種的干旱脅迫緊密相關或響應miRNAs已被鑒定和研究[2-4]。然而，干旱逆境適應和響應是一個非常復雜的過程，更多干旱調控關鍵miRNAs及其功能有待鑒定和解析，而且miRNA-靶基因調控干旱脅迫的分子機制仍不完全清楚。

本研究擬通過高通量RNA-seq技術對水稻根、莖和葉片中small RNA在干旱處理后全基因組水平上的表達變化情況進行研究，鑒定已知miRNAs和新miRNAs、篩選和挖掘干旱響應差異表達miRNAs并對其預測的調控靶基因進行功能注釋和富集分析，為進一步探明水稻干旱脅迫應答和適應過程中的miRNAs調控機制及其分子調控通路和網絡提供理論和數據參考。

1 材料與方法

1.1 水稻材料培養

水稻材料OsPIR1-RNAi(背景為OryzasativaL. ssp.japonica)最初來自本課題組構建的RNAi突變體庫[8]，為OsPIR1基因(Os03g0845000)的RNA干擾(RNA interference)植株。首先，選取該材料T2代的種子經75%乙醇和2.5% NaClO(體積分數)消毒后在濾紙上室溫萌發約48 h，然后轉移至土壤，在16 h光照/8 h黑暗，光照強度約8 000 lx，28 ℃光照培養箱生長至幼苗二葉期。

1.2 干旱處理與樣品收集

選取二葉期生長相對一致的水稻植株轉移至中國農業科學院生物技術所溫室內繼續培養，從該時期的水稻幼苗開始持續斷水處理。處理0、24、96 h采集根(root，R)、莖(stem，S)及葉(leaf，L)，分別命名為R-0 h、S-0 h、L-0 h、R-24 h、S-24 h、L-24 h、R-96 h、S-96 h、L-96 h，共9個樣品，立即置于液氮中保存備用。

1.3 Small RNA文庫構建及測序

采用EASYspin Plus多糖多酚復雜植物RNA快速提取試劑盒(RN5301，Aidlab)提取9個樣品的總 RNAs，通過反轉錄合成cDNA第一鏈。PCR擴增后經過純化和大小選擇構建small RNA文庫并進行質量評估[9-10]。將構建合格的文庫使用北京貝瑞和康生物技術有限公司Illumia Solexa測序平臺進行高通量single-end測序。

1.4 測序數據質量控制及注釋

首先將Illumia測序得到的圖像經堿基識別[11]轉化為FASTQ格式的原始測序數據(raw reads)。去除原始數據中低質量讀數、接頭序列、冗余片段及短于18或長于30個核苷酸的序列，獲得高質量序列(clean reads)，統計測序堿基質量值Q30等[11-12]，盡可能地保證測序結果的準確性。將clean reads利用Bowtie軟件(v1.0.0)[13]比對Silva(http://www.arb-silva.de/)、GtRNAdb(http://lowelab.ucsc.edu/GtRNAdb/)、Rfam(http://rfam.xfam.org/)和Repbase(http://www.girinst.org/repbase/)等數據庫進行注釋，包括核糖體RNA(rRNA)、核內小RNA(snRNA)、核仁小RNA(snoRNA)、轉運RNA(tRNA)等非編碼RNA以及重復序列(repbase)，未比對上的clean reads則被認定為未注釋序列(unannotated reads)。

1.5 已知miRNAs和新miRNAs鑒定

使用Bowtie軟件[13]將未注釋序列與水稻參考基因組(Oryza_sativa.IRGSP-1.0)進行序列比對，獲取其在基因組上的位置信息。將比對到參考基因組的序列(mapped reads)與miRNA數據庫miRBase(release 22；http://www.mirbase.org/)中的已知成熟miRNA序列及其上游2 nt與下游5 nt范圍內進行比對，最多允許一個堿基錯配，能比對到的序列即被鑒定為已知miRNAs。利用miRDeep2軟件(v2.0.5)[14]，通過mapped reads在參考基因組上的位置信息得到潛在前體序列，基于miRNA前體特征二級發夾結構及前體結構能量信息(RNAfold randfold)采用貝葉斯模型經打分最終實現新miRNAs的預測[15]。

1.6 差異表達miRNAs分析

首先將所有樣本中miRNAs表達量采用TPM(transcript per million)值進行標準化或歸一化處理[16]。將干旱處理后24 和96 h的根莖葉(R/S/L-24 h、R/S/L-96 h)分別與0 h對應的組織(R/S/L-0 h)進行兩兩比較，使用edgeR軟件(v3.8.6)[17]用來確定差異表達，以|log2(FC)|≥1.0(fold change, FC)和FDR≤0.1(false discovery rate, FDR)為標準篩選差異表達的miRNAs。利用在線平臺(http://bioinformatics.psb.ugent.be/webtools/Venn/)構建韋恩圖。使用R packages(http://www.rproject.org/)對差異表達的miRNAs以log2(TPM+1)值進行層次聚類分析，TPM (transcripts per million)為標準化的miRNA表達值。

1.7 靶基因預測及功能注釋

利用TargetFinder平臺(v1.6)[18-19]對獲得的差異表達的已知miRNAs和新miRNAs靶基因進行預測。然后通過BLAST軟件將靶基因序列與NR(ftp://ftp.ncbi.nih.gov/blast/db/)、Swiss-Prot(http://www.uniprot.org/)、GO (http://www.geneontology.org/)、COG(http://www.ncbi.nlm.nih.gov/COG/)、KEGG(http://www.genome.jp/kegg/)、KOG (ftp://ftp.ncbi.nih.gov/pub/COG/KOG/)、Pfam(http://pfam.xfam.org/)等數據庫進行檢索比對，進而對靶基因進行功能注釋。

1.8 靶基因GO與KEGG富集分析

使用GOseq[20]和KOBAS[21]對靶基因進行GO富集分析(Wallenius non-central hyper-geometric distribution)和KEGG富集分析。

2 結果與分析

2.1 水稻干旱處理后small RNA的數量變化

干旱處理0、24和96 h的水稻根組織分別獲得4 120 407、2 538 241和2 130 755條高質量序列(clean reads)；莖組織分別產生2 506 166、1 371 020和1 407 041條高質量序列；葉組織分別產生2 995 912、1 741 575和1 617 673條高質量序列。通過質量控制，每個樣本堿基質量值Q30均大于94%，表明small RNA-seq建庫測序過程中堿基識別精度很高，數據質量可靠。

干旱處理后0、24和96 h的水稻根組織產生的未注釋序列(unannotated reads)占所有small RNA比例分別為87.15%、72.19%和76.57%；莖的占比分別為90.4%、71.78%和77.48%；葉的占比分別為92.06%、69.98%和79.34%(表1)。將以上得到的unannotated reads映射到水稻參考基因組，能比對到基因組正鏈(+)或負鏈(-)上的即為mapped reads。干旱處理0、24和96 h的水稻根組織mapped reads占所有unannotated reads比例分別為55.75%、13.42%和21.04%；莖的占比分別為60.05%、50.77%和43.78%；葉的占比分別為53.18%、39.04%和46.05%(表1)。這些結果表明，干旱處理后水稻根、莖和葉組織中的clean reads、unannotated reads、mapped reads均存在不同程度的下降。

表1 水稻small RNA分類注釋和比對到水稻參考基因組信息統計Table 1 Summary of classification of rice small RNA and unannotated reads mapped on rice genome

2.2 水稻干旱處理后small RNA的長度分布變化

small RNA的長度一般為18～30 nt，其中21～24 nt的small RNA在植物體內豐度較高且發揮了重要的調控作用。對干旱處理后水稻根、莖和葉組織中small RNA長度分布的變化趨勢進行分析，結果(圖1)表明：①R/S/L-0 h、S/L-24 h和R/S/L-96 h中24 nt small RNA豐度最高，而R-24 h中21 nt small RNA豐度最高；②隨著干旱處理時間的延長，水稻根和莖中21 nt small RNA豐度先升高后下降，而葉片中的21 nt small RNA則不斷下降；③根、莖和葉中22 nt small RNA豐度在處理后均呈現持續下降趨勢；④隨著干旱處理時間的增加，根中23 nt small RNA豐度先迅速下降后緩慢上升，而莖和葉中23 nt small RNA豐度則不斷下降；⑤根、莖和葉中24 nt small RNA的豐度在干旱處理后均呈先下降后上升趨勢。

圖1 干旱處理后水稻根、莖和葉中small RNA長度分布Fig.1 Length distribution of small RNAs in rice root, stem and leaf under drought

2.3 水稻干旱處理后miRNAs變化

將mapped reads與miRBase數據庫進行序列比對和莖環結構預測，結果(表2)發現，所有樣品共鑒定到403個miRNAs，其中已知miRNA 352個，新預測miRNA 51個，而且水稻干旱處理后不同時間的根、莖和葉組織中已知miRNAs和新miRNAs表達數量不同。

表2 miRNAs數量統計Table 2 Summary of miRNAs number

2.4 水稻干旱處理后差異表達的miRNAs

根據表達量計算各miRNA在同一組織不同時間點樣本間的表達水平變化情況，結果(表3)顯示，水稻葉片干旱處理96 h后差異表達的miRNAs最多，而根部處理96 h后和莖部處理24 h后差異表達的miRNAs最少。其中，miRNAs表達上調數量最多的是干旱處理后96 h的葉，數量最少的是處理后24 h的根和莖；miRNAs表達下調數量最多的是干旱處理后96 h的莖部，數量最少的是處理后96 h的根部。為深入挖掘參與調控水稻干旱脅迫響應過程中的miRNAs，進一步對差異表達miRNAs在干旱處理后不同組織的變化趨勢進行詳細分析。

表3 干旱處理后水稻根、莖和葉差異表達的miRNAs數量Table 3 Number of differentially expressed miRNAs in drought-treated root, stem and leaf

2.4.1根差異表達miRNAs分析 從圖2可以看出，21個miRNAs在處理后兩個時間點的水稻根組織中都差異表達。根據表達量變化趨勢，將這21個差異表達miRNAs分為三類。第Ⅰ類包含Osa-miR11337-5p和Osa-miR11340-5p，它們在干旱處理24 h內表達量迅速降低，在24 至96 h內逐漸升高；第Ⅱ類包含12個miRNAs，它們的表達量在干旱處理96 h內持續升高；第Ⅲ類包含7個miRNAs，它們的表達量在干旱處理24 h內迅速升高，在24 至96 h內逐漸下降。

圖2 干旱處理后水稻根中差異表達miRNAs情況Fig.2 Changes in miRNAs expression of root under drought

對這21個差異表達miRNAs進行分析發現，Osa-miR166家族和Osa-miR820家族中表達發生變化的成員最多，其中Osa-miR166家族的4個成員表達量呈現先升高后下降趨勢，而Osa-miR820家族的3個成員表達則持續升高，表明這兩個miRNA家族可能在水稻根對干旱脅迫應答反應過程中發揮重要功能。

2.4.2莖差異表達miRNAs分析 從圖3可以看出，28個miRNAs在處理后兩個時間點的水稻莖組織中都呈現差異表達。根據表達水平變化趨勢，將這28個差異表達miRNAs分為四類。第Ⅰ類包含4個miRNAs：Osa-miR167i-3p、Osa-miR408-5p、Osa-miR2878-5p和Osa-miR11340-5p，它們的表達量在干旱處理96 h內持續下調；第Ⅱ類包含3個miRNAs：Osa-miR1425-5p、Osa-miR5788和Osa-miR11337-5p，它們的表達量在24 h內迅速降低，在24 至96 h內逐漸升高；第Ⅲ類包含10個miRNAs，它們在干旱處理96 h內表達量持續升高；第Ⅳ類包含了11個miRNAs，它們的表達量在96 h內呈現先升高后降低的趨勢。

差異表達的miRNA家族中，Osa-miR166家族的5個成員和Osa-miR159家族的2個成員在干旱處理后表達量呈現先升高后下降趨勢，而3個新miRNAsnovel_miR_13、novel_miR_24和novel_miR_29則受干旱脅迫誘導持續上調表達，推斷它們可能也參與了水稻莖部干旱響應分子調控網絡。

2.4.3葉差異表達miRNAs的變化 從圖4可以看出， 38個miRNAs在處理后兩個時間點的水稻葉組織中都呈現差異表達。根據表達量變化趨勢，將這38個差異表達miRNAs分為五類。第Ⅰ類包含8個miRNAs，它們的表達水平在干旱處理96 h內持續下降；第Ⅱ類包含4個miRNAs，它們在24 h內表達量迅速降低，在24 至96 h內逐漸升高；第Ⅲ類包含19個miRNAs，它們在干旱處理96 h內表達量持續上調；第Ⅳ類包含5個miRNAs，它們的表達量在96 h內先升高后降低；第Ⅴ類包含Osa-miR1425-3p和Osa-miR11342-5p，它們只在0 h的葉中表達，而在處理后24 和96 h均未檢測到表達。

對這38個差異表達miRNAs進行進一步家族分析發現，Osa-miR160家族的3個成員、Osa-miR166家族的5個成員和Osa-miR820家族的3個成員表達量受干旱誘導持續上升，而Osa-miR159家族的2個成員表達水平隨著處理時間的增加呈現先升高后降低趨勢，表明這些miRNA家族在水稻葉片中響應干旱脅迫。

2.4.4根、莖和葉中共同差異表達miRNAs的變化 進一步對21、28和38個分別在水稻根、莖和葉差異表達的miRNAs進行分析，結果(圖5)發現，共有15個miRNAs在干旱處理后24和96 h的各組織中均差異表達。其中Osa-miR159f、Osa-miR164f、Osa-miR5082和Osa-miR5493的表達水平隨處理時間的增加而持續升高，表明它們受干旱脅迫誘導。因此，這些miRNAs是水稻營養組織干旱逆境適應與響應機制中關鍵調控因子。進一步對比水稻不同組織的差異表達miRNAs發現：Novel_miR_37在0 h的根、莖和葉中表達水平很低，而干旱處理后迅速上升；Osa-miR11337-5p和Osa-miR11340-5p在0 h的根、莖和葉中表達量較高，而處理后迅速降至很低水平；Osa-miR166c-3p在根、莖和葉中的表達變化趨勢相同，均表現為先升后降；Osa-miR1876在根和葉中表達量隨處理時間持續上升，而在莖中則先顯著上調后緩慢下降。這些結果說明，不同處理時間的根、莖和葉組織具有明顯不同的miRNAs表達譜，而且同一個miRNA在不同組織部位的表達水平也并不相同。

圖5 干旱處理后水稻根、莖、葉共同差異表達的miRNAsFig.5 Differentially expressed miRNAs all in root, stem and leaf under drought treatment

2.5 預測的miRNAs靶基因及GO/KEGG富集

miRNA是通過調節靶基因的表達水平從而發揮生理功能，靶基因被進一步利用TargetFinder平臺進行預測，最終發現403個miRNAs共對應4 537個靶基因。其中，352個已知miRNAs共對應3 883個靶基因，51個新miRNAs共對應785個靶基因。進一步使用數據庫比對方法獲得靶基因的詳細功能注釋信息，4 537個靶基因中4 484個獲得注釋信息。在干旱處理24和96 h后,根中差異表達的51個miRNAs共對應923個靶基因，莖中差異表達的69個miRNAs共對應641個靶基因，葉中差異表達的83個miRNAs共對應1 247個靶基因。

為了進一步了解干旱處理后水稻根、莖和葉組織中差異表達miRNAs對應的靶基因主要富集到哪些代謝途徑和生物學過程，進行了GO和KEGG富集分析。結果(圖6)表明：對于根，在GO生物學過程中靶基因顯著富集于氧化還原、木質素分解代謝等過程，在GO細胞組分中顯著富集于細胞膜相關囊泡和質外體等；在GO分子功能中顯著富集于電子載體活性和血紅素結合等，而角質、木栓質和蠟質生物合成及堿基切除修復等是顯著富集的KEGG通路；對于莖，在GO分子功能中顯著富集于ADP結合和染色質結合等，在GO生物學過程中靶基因顯著富集于防御應答和細胞凋亡過程等，在GO細胞組分中顯著富集于核仁組織區和葉綠體基質類囊體等，而甘氨酸、絲氨酸和蘇氨酸代謝及脂肪酸生物合成等是顯著富集的KEGG通路；對于葉，在GO生物學過程中靶基因顯著富集于氧化還原和細胞凋亡過程等，在GO分子功能中顯著富集于染色質結合和質子反向轉運活性等，在GO細胞組分中顯著富集于泛素連接酶復合體和細胞膜相關囊泡等，而角質、木栓質和蠟質生物合成以及油菜素內酯生物合成等是顯著富集的KEGG通路。總之，這些結果為揭示水稻根、莖和葉組織響應干旱脅迫應答的分子細胞機制提供了線索和參考。

A:根；B：莖；C：葉；左圖：GO分析；右圖：KEGG富集分析。A: Root; B: Stem: C: Leaf; Left: GO analysis; Right: KEGG analysis.圖6 干旱處理條件下差異表達miRNAs對應靶基因GO和KEGG分析Fig.6 GO and KEGG analysis of genes targeted by differentially expressed miRNAs after drought stress

3 討論

miRNA不僅受干旱脅迫誘導而改變表達水平，也通過在轉錄和轉錄后水平調節干旱響應靶基因表達而參與植物干旱環境適應及干旱抗性應答過程[2, 5]。目前，已有研究使用芯片技術對斷水干旱處理條件下分蘗期和花序形成期的水稻植株[22]，以及利用small RNA測序對20% PEG模擬干旱脅迫條件下幼苗期的水稻根和苗組織[23]和斷水干旱處理條件下開花期的水稻劍葉、穗[24]與成熟花序[25]進行了表達譜分析。本研究通過分析全基因組范圍內miRNA轉錄組的變化深入挖掘干旱逆境相關miRNAs，并結合先前研究結果解析水稻不同組織響應干旱脅迫的調控網絡。

miR159受干旱脅迫誘導而上調表達，而且它通過抑制MYB類靶基因的表達來參與ABA信號轉導途徑進而在干旱響應過程中發揮重要功能[3, 26]。本研究結果顯示，Osa-miR159f的表達水平在營養組織中受干旱脅迫誘導持續上調，而Osa-miR159f的預測靶基因Os01g0812000、Os03g0578900、Os05g0490600和Os06g0605600編碼MYB類轉錄因子，其中Os03g0578900和Os06g0605600已被報道在干旱條件下的水稻葉和花序組織中表達下調[27-28]，表明干旱脅迫下Osa-miR159f與Os03g0578900和Os06g0605600可能存在負調控關系。miR164通過調節NAC類靶基因的表達水平來參與生長素信號并負調控干旱抗性[29-30]。本研究中Osa-miR164f在干旱處理后持續上調表達，而其預測靶基因Os02g0579000、Os04g0460600、Os06g0675600、Os08g0200600和Os12g0610600編碼NAC類轉錄因子，其中Os02g0579000和Os08g0200600的表達在水稻葉或幼苗中受干旱誘導下調[27, 31]，表明Osa-miR164f在干旱響應過程中可能負調控Os02g0579000和Os08g0200600。已有研究表明，miR166通過轉錄后調控Homeobox-leucinezipper(HD-ZIP)類靶基因的表達參與干旱脅迫下側根發育[2, 4]。本研究中五個miR166家族成員在干旱處理后的營養組織中都差異表達，而它們的預測靶基因編碼HD-ZIP類轉錄因子，其中Os12g0612700在干旱處理后的水稻葉中表達被下調[27]，表明在干旱響應過程中它們可能也存在負調控關系。因此上述調控模塊在水稻營養組織響應和抵御干旱脅迫過程中可能發揮關鍵作用。miR160、miR167、miR408和miR528也被證實涉及響應、緩解和抵御植物干旱逆境脅迫[2-4]，而本研究中這些miRNAs也在干旱處理后根、莖或葉組織中也差異表達。Osa-miR5082和Osa-miR5493在營養組織中受干旱誘導顯著上調表達，但它們在干旱脅迫調控網絡中的具體功能目前尚未被報道。而且本研究中發現的許多差異表達的新miRNAs可能與干旱脅迫調控緊密相關，它們的具體分子機制也有待于進一步挖掘和解析。

植物在遭遇干旱脅迫時會通過干旱響應miRNAs-靶基因在分子水平啟動干旱適應機制，進而改變生理代謝活動和形態結構如抗氧化物和滲透調節物質積累、細胞失水保護物質合成等來緩解或抵御干旱逆境[2,5,32-34]。本研究中差異miRNAs對應的靶基因顯著富集在氧化還原途徑、脯氨酸和α-亞麻酸(alpha-linolenic acid)代謝以及角質、木栓質和蠟質生物合成等，證實這些生物過程和代謝途徑可能參與干旱逆境響應，但它們的具體作用機制仍待進一步探究。

總之，miRNAs在水稻根、莖和葉組織對干旱脅迫的適應和響應過程中扮演了重要的角色，而且本文構建的miRNAs-靶基因的調控模塊也為進一步闡明和完善植物在干旱脅迫條件下結構、生理和分子改變的詳細機理提供了理論基礎，并將為分子遺傳育種改良水稻抗旱性提供了分子數據資源。