白羊草葉片和根系干旱脅迫下microRNAs差異表達分析

李春艷， 董 潔， 鐘 華， 董寬虎*

(1.山西農業大學動物科技學院， 山西 太谷 030801；2. 山西省農業科學院畜牧獸醫研究所， 山西 太原 030032； 3.北京市科學技術情報研究所， 北京 100044)

眾所周知，干旱脅迫常影響植物生長發育，對農作物造成的損失在所有的非生物脅迫中占首位[1]。植物在漫長的進化過程中形成了適應各種脅迫的機制，許多脅迫響應正是通過誘導或抑制miRNAs的表達進而調控相關基因表達來實現的[2]。miRNAs是一類長度約為19～25nt的非編碼的單鏈小分子RNAs，廣泛存在于真核生物中。近年來與干旱脅迫相關的miRNA及其靶基因已經在不同植物中得到了驗證，如擬南芥[3]、棉花[4]、水稻[5]、苜蓿[6]和小麥[7]等，但這僅限于一些全基因組已知的物種，對于其他植物則研究的很少。

白羊草(B.ischaemum)為禾本科孔穎草屬多年生植物，分布廣泛，具有固土保水、生命力強，耐旱、高產耐牧、耐踐踏、適口性好等優點[8]。多年來對于白羊草的研究主要集中在生產性能[9-10]、營養價值[11]、生理生態[12]及遺傳多樣性上[13]，隨著分子生物技術的發展，對于干旱脅迫下白羊草基因的分離和克隆等研究取得一些成績[14]，但對白羊草耐旱的分子機理及相關miRNAs研究甚少，因此，本研究在結合干旱脅迫條件下白羊草葉片和根系轉錄組數據的基礎上，運用高通量Illumina測序技術對干旱脅迫和正常生長條件下白羊草葉片和根系的miRNAs進行測序分析，篩選出響應干旱脅迫的miRNAs，并對其預測的靶基因進行功能分析，為探索白羊草耐旱的分子機制和提供可利用的耐旱基因資源奠定理論基礎。

1 材料與方法

1.1 材料的制備

試驗材料為采集于太谷縣的白羊草種子，盆栽試驗在山西農業大學草業科學的日光溫室中進行，溫室中的平均溫度為15～30℃，平均光照強度為26 278 Lx，相對濕度為50%～60%。土壤取自山西農業大學動物科技試驗站牧草試驗田的耕層土(0～20cm)，pH為7.5，田間最大持水量為23.91%。每盆裝風干土7 kg，播種量為50粒，正常澆水，苗齊后間苗，每桶留20株長勢一致的幼苗。待白羊草幼苗平均生長至20 cm左右時開始干旱處理，處理前一次性澆水，使桶內土壤含水量保持在19.13±3.38%(相當于田間持水量的80%)。對照組持續每天澆水，使土壤含水量保持在田間含水量的80%。實驗組不澆水，當葉片開始出現萎蔫土壤的含水量為10.04±2.73%(相當于田間持水量的42%)時，取長勢一致的幼苗，分別剪取每株上全部的葉片和根系；對照組取同期沒有進行干旱處理的正常生長幼苗的葉片和根系，每個處理重復3次，液氮速凍后—70℃冰箱保存，用于RNA提取。

為了減少隨機取樣帶來的偏差，每次取樣從3株獨立的個體中取出，提取RNA后等量混合用于后續的建庫測序。在這次試驗中，我們構建了4個白羊草的小RNA文庫，干旱脅迫后的葉片和根系及其對照組的葉片和根系各1個(標記為TL-3、TR-3、CKL-3和CKR-3)。

1.2 樣品RNA提取與文庫的構建

采用Trizol法提取白羊草干旱脅迫和對照組的葉片和根系樣品總RNA，用1%的瓊脂糖電泳檢測RNA樣品是否有降解以及雜質；用2100 RNA Nano 6000 Assay Kit檢測RNA樣品的完整性和濃度，然后用凱奧K5500分光光度計檢測樣品純度。

總RNA樣本檢測合格后，首先對總RNA進行片段選擇，通過膠分離技術收集18～30nt RNA片段；在分離得到的RNA片段的兩端分別連接上接頭，然后反轉錄成cDNA，進行PCR擴增建立測序文庫；最后，對檢驗合格的測序文庫進行Illumina HiSeq高通量測序，采用SE50測序策略[15]。

1.3 白羊草miRNAs的鑒定、量化和差異表達分析

Illumina測序所得原始序列，通過去除低質量、接頭污染和含未知堿基比例大于10%的reads等過程完成數據處理，得到用于后續分析的目標序列。將過濾后的clean reads根據進一步的長度篩選(植物18～30nt)，因數據庫中沒有白羊草基因組和EST序列信息，選擇模式植物擬南芥的基因組參考序列，條件設置為允許1個錯配。通過基因組比對分析軟件SOAP(short oligonucleotide alignment program)將符合條件的clean reads與參考序列進行比對，得到定位信息[16]。匹配較好的reads用于后續分析。之后，將Rfam數據庫和GenBank的ncRNA序列比對到基因組[17]，得到該物種ncRNA在基因組中的定位信息，將匹配好的Clean reads和ncRNA的定位信息進行匹配，分別注釋為rRNA、tRNA、snRNA、snoRNA和其他ncRNA的小RNA。然后將miRbase數據庫中的miRNAs序列比對到基因組中，得到該物種目前已知的所有miRNAs在基因組中的定位信息，將剩余未注釋匹配好的Clean Reads與miRNAs的定位信息進行完全匹配[18]，就可以鑒定已知miRNAs。對于不能與已知注釋區域匹配的Clean Reads，采用軟件miRDeep2的方法進行新miRNA的預測[19]。針對每個樣本，統計與miRNAs匹配的Total Clean Reads的數量，并計算歸一化后的表達量RPM(Reads Per Million)值。然后使用DEGseq軟件包進行歸一化后的差異表達分析，即選取|log2Ratio|≥1且q<0.05的miRNAs作為差異表達miRNAs[20]。

1.4 miRNAs靶基因的預測與功能富集分析

通過psRobot[21]對差異表達顯著的miRNAs進行靶基因的預測，并對這些靶基因做GO(Gene Ontology)功能和KEGG(The Kyoto Encyclopedia of Gene and Genome)代謝通路富集分析。GO有三個類別，分別描述基因的分子功能、所處的細胞位置和參與的生物過程。對于差異表達顯著的miRNAs所預測的靶基因，可以采用Blast2GO得到每個基因對應的GO條目[22]。將上述靶基因比對到KEGG數據庫(http://www.geneome.jp/ kegg)進行通路富集分析，利用KAAS軟件找出顯著性富集的通路[23]。

1.5 用qRT-PCR法驗證miRNAs的表達

為了對miRNAs測序數據進行驗證，我們采用qRT-PCR法對隨機選取的8個miRNAs表達量進行測定。首先設計特異性的反轉錄引物RT primer(reverse transcription primer)(表1)，然后以建庫所剩的總RNA為模板，按照說明書的指導使用M-MuLV Reverse Transcription進行cDNA的合成。然后使用每個miRNA特異性的正向引物和通用的反向引物(表1)進行PCR擴增，qRT-PCR的反應體系為：總體積為20 μl，包括2 μl的cDNA，0.8 μl引物混合物，10μl的1×SYBR Mix。反應程序為兩步法：第一步95℃ 3 min，第二步95℃ 5 s，60℃ 20 s，40個循環。每個樣品3個重復，內參基因為18 s rRNA[14]，基因的相對表達量使用2-ΔΔCt法進行[24]。

表1 qRT-PCR驗證實驗中選取的miRNAs引物序列Table 1 Primer sequences of selected miRNAs used in qRT-PCR validation experiment

2 結果與分析

2.1 小RNA測序數據分析

本次測序共獲得61 609 058 raw reads和43 566 649 clean reads，堿基Q30(堿基錯誤率小于0.1%)都在95%以上(表2)。我們分析了18～30nt長度分布的小RNA，unique clean reads的長度分布24nt長度具有明顯的峰值(圖1)，這與其他植物的小RNA高通量測序研究結果相一致。

表2 白羊草葉片和根系測序數據統計Table 2 Summary of sequencing in B.ischaemum leaves and roots

圖1 白羊草葉片和根系小RNA過濾后reads種類長度分布Fig.1 Length distribution of unique small RNAs in B.ischaemum leaves and roots

2.2 小RNA分類注釋和miRNAs鑒定

Clean Reads中所有比對上小RNA的結果見表2，被注釋為miRNAs的平均數量為3 201 445個，達到了總數的7.35%。這表示本實驗所構建的小RNA文庫已經富集到了白羊草的miRNAs。然而unann的數量平均占到小RNA總數的89.50%，在所有的小RNA種類中最高，這說明本實驗所構建的小RNA文庫中含有未知的小RNA和潛在的新miRNAs。

總的來說，這次試驗鑒定出白羊草中有屬于20個miRNAs家族的79個已知miRNAs，在這20個家族中，最大的家族是含有13個成員的ata-miR169家族，其次是ata-miR396含有8個成員，再次是ata-miR156有和ata-miR167均含有7個成員。表達量超過10 000的miRNAs家族有6個，為一些在其它物種中也很保守miRNAs種類。此外，鑒定出92種新miRNA，其中表達量超過500的新miRNAs有10種。

表3 白羊草葉片和根系中小RNA的分類Table 3 Small RNA categorization inB.ischaemum leaves and roots

2.3 miRNAs在干旱脅迫下的差異表達分析

干旱脅迫后葉片和根系miRNAs差異表達數為分別為65和27個(表4)，葉片和根系中組織差異的miRNAs分別為55個和17個，共有的差異表達miRNAs為10個(圖2)。在這些差異表達的miRNAs中，除了ata-miR156c-3p是在干旱脅迫后的葉片中下調，根系中上調表達外，其余的9個miRNAs在干旱脅迫后的葉片和根系表達模式是一致的，表現為miR164家族和Novel 40均上調，剩余6個miRNAs均下調。

表4 干旱脅迫條件下白羊草葉片和根系差異表達miRNAs的個數Table 4 The number of differentially expressed miRNAs in leaves and roots of B.ischaemum under drought stress

圖2 白羊草葉片和根系中共同響應干旱脅迫的miRNAs的表達水平Fig.2 The expression level of drought-responsive miRNAs common in leaves and roots ofB.ischaemum

2.4 差異表達miRNA靶基因預測及GO功能富集分析

miRNAs調控的靶基因鑒定是了解miRNA功能的關鍵。本次試驗干旱脅迫后葉片和根系已知miRNAs預測的靶基因分別為430個和158個，新miRNAs所預測到的靶基因則多達幾千個，結合轉錄組注釋結果[25]，發現預測的這些靶基因多為轉錄因子、其他的轉錄調節物質和一些與脅迫相關的蛋白及與激素信號轉導相關的物質。

對這些預測到的靶基因進行GO富集分析(圖3)，發現無論是葉片還是根系，miRNAs預測到的靶基因在生物學過程中富集基因數最多的是cellular process(GO:0009987)，metabolic process(GO:0008152)，single-organism process(GO:0044699)，分子功能中富集最多的是binding(GO:0005488)，catalytic activity(GO:0003824)，在細胞組分中富集基因數最多的cell part(GO:0044464)，organelle(GO:0043226)，membrane(GO:0016020)。其中“cellular process”，“metabolic process”，“cell part”，“binding”和“catalytic”注釋到的靶基因數均占所有注釋靶基因數的50%以上。

圖3 干旱脅迫條件下白羊草差異表達miRNAs預測靶基因的GO功能富集分析Fig.3 Function analysis of target transcripts of differentially expressed miRNAs inB.ischaemum under drought stress.注：橫坐標為Ontology分類，縱坐標為注釋到該class中的靶基因占所有注釋的靶基因的比例Note：The X-axis is the Ontology classification,and the Y-axis is the ratio of the target gene annotated in the class to the target gene annotated in all annotations

2.5 KEGG富集分析

本實驗檢測出367條通路KEGG通路，其中與干旱脅迫相關且靶基因顯著富集的的通路在干旱脅迫后的葉片中有10條，在干旱脅迫后的根系中有1條(表5)。從表5中我們可以看出，Plant hormone signal transduction是干旱脅迫后葉片和根系共有的主要富集代謝通路，而Flavonoid biosynthesis則僅在干旱脅迫后葉片中顯著富集且具有最高的Rich factor(指差異表達的基因中位于該pathway條目的基因數目與所有注釋基因中位于該pathway條目的基因總數的比值)。

表5 干旱脅迫條件下白羊草葉片和根系(A：葉片；B：根系)代謝通路富集分析Table 5 The enriched pathway in leaves and roots (A:leaves;B:roots) of B.ischaemum under drought stress

注：q≤0.05作為KEGG通路顯著富集的標準，q代表調整后的P值，a代表富集到某一特定通路的差異表達基因的數目，b代表富集到該通路上的所有unigene數

Note：The significantly enriched KEGG pathways were determined using q≤0.05,q represents the corrected P-value,a Number of DEGs assigned to certain KEGG pathway,b Number of all reference unigenes assigned to certain KEGG pathway

2.6 miRNAs的qRT-PCR驗證

我們對隨機選取的8個具有顯著差異表達的miRNAs進行定量qRT-PCR驗證，從高通量測序miRNAs表達量(Fig.4A)和qRT-PCR(Fig.4B)的結果可以看出，這8個miRNAs在高通測序和qRT-PCR檢測中表現出相似的表達模式。

3 討論

MiRNAs在植物適應逆境脅迫的過程中發揮著重要作用。有一些miRNAs家族在不同植物物種之間是保守的，甚至在一些進化距離較遠的苔蘚與開花植物之間依然存在保守性[26]，但這些保守的miRNA在不同物種和組織中表達是不同的。如干旱脅迫條件下miR398在西紅柿[27]和棉花[4]中是下調的，而在小麥[28]和蒺藜苜蓿[6]中則是上調表達。這次研究發現miR398g-3p和miR398f-3p在干旱脅迫后的葉片和根系中均是下調。Kantar在對干旱脅迫后大麥葉片和根系miRNAs的研究中，發現四種miRNAs在脫水條件下呈現出組織特異的調控：miR166在葉中上調根中下調；而miR156a，miR171與miR408在葉中發生了誘導改變表達而在根中則沒有變化[29]。這與水稻干旱脅迫后的結果有一些差異，水稻不論是地上部分還是根系在干旱脅迫后miR156和miR171的表達均下調[30]。在本試驗中，miR156家族在干旱脅迫后葉片中下調根中上調，miR167，miR172，miR390，miR394，miR408在葉中發生了誘導改變表達在根中沒有發生表達變化，而miR166，miR171則是在根中發生了誘導改變表達在葉中沒有發生表達變化。這些具有沖突性的結果表明miRNAs的復雜調控并非僅與不同的物種和脅迫條件及組織差異性有關，也與物種在長期進化適應過程中所形成的同一家族中不同miRNAs成員的差異有關。

圖4 白羊草葉片和根系響應干旱脅迫的miRNA定量PCR驗證Fig.4 qRT-PCR validation of drought-responsive miRNAs in leaves and roots of B. ischaemum

研究表明，保守的miRNAs在大量物種中都有保守的靶向基因[31]，這表明miRNAs與其所調控靶基因之間的關系在植物長期進化的過程中趨于穩定[32]，如miR156的靶基因SPL (squamosa promoter-binding-like protein)在植物發育過程、合成代謝及非生物脅迫中起到了重要的作用，成為植物生長與發育的調控樞紐[33]。我們這次對白羊草miR156靶基因的預測表明，SPL是miR156所預測的主要靶基因。但除SPL外，我們所預測的miR156還有一些新的靶基因如脂質轉運蛋白、質膜ATPase等。預測的miR164的靶基因除了在大多數植物中已經驗證了的NAC轉錄因子，還有絲氨酸-蘇氨酸蛋白激酶和假想蛋白等。雖然這些新預測的靶基因還有待進一步的驗證，但靶基因預測的結果給我們一個暗示，一個miRNA家族中不同的成員可能能夠調控各自不同但是較為相似的靶基因，從而參與多個生理生化途徑的調控。

雖然一些新預測的miRNAs的成熟序列在miRBase中并未找到對應的序列，但通過對其預測的靶基因進行GO功能聚類分析和KEGG代謝途徑分析得知，這些新miRNAs主要在植物生長、發育、脅迫響應和激素信號轉導中發揮作用。如miRNA＿39被預測可調控植物激素信號轉導中ABA信號通路，miRNA＿36被預測與白羊草葉片響應干旱脅迫的黃酮類化合物代謝途徑密切相關。這說明這些新預測的miRNAs對白羊草適應干旱脅迫有著重要的作用，研究和探索這些新預測miRNAs的表達模式和功能對進一步闡明白羊草適應干旱的分子機制具有重要意義。

4 結論

這次通過對干旱脅迫條件下白羊草葉片和根系miRNAs表達特性的研究發現，白羊草響應干旱脅迫的miRNAs除了一些在其他植物中已經被證實與干旱脅迫相關的miRNAs外，還發現一些與干旱脅迫相關的新miRNAs，通過對miRNAs靶基因生物學功能和代謝途徑的富集分析發現，植物的次生代謝產物尤其是黃酮類化合物合成和植物信號轉導途徑是白羊草響應干旱脅迫富集的主要通路。