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

(1.山西農(nóng)業(yè)大學(xué)動物科技學(xué)院， 山西 太谷 030801；2. 山西省農(nóng)業(yè)科學(xué)院畜牧獸醫(yī)研究所， 山西 太原 030032； 3.北京市科學(xué)技術(shù)情報(bào)研究所， 北京 100044)

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

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

1 材料與方法

1.1 材料的制備

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

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

1.2 樣品RNA提取與文庫的構(gòu)建

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

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

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

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

1.4 miRNAs靶基因的預(yù)測與功能富集分析

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

1.5 用qRT-PCR法驗(yàn)證miRNAs的表達(dá)

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

表1 qRT-PCR驗(yàn)證實(shí)驗(yàn)中選取的miRNAs引物序列Table 1 Primer sequences of selected miRNAs used in qRT-PCR validation experiment

2 結(jié)果與分析

2.1 小RNA測序數(shù)據(jù)分析

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

表2 白羊草葉片和根系測序數(shù)據(jù)統(tǒng)計(jì)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的結(jié)果見表2，被注釋為miRNAs的平均數(shù)量為3 201 445個(gè)，達(dá)到了總數(shù)的7.35%。這表示本實(shí)驗(yàn)所構(gòu)建的小RNA文庫已經(jīng)富集到了白羊草的miRNAs。然而unann的數(shù)量平均占到小RNA總數(shù)的89.50%，在所有的小RNA種類中最高，這說明本實(shí)驗(yàn)所構(gòu)建的小RNA文庫中含有未知的小RNA和潛在的新miRNAs。

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

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

2.3 miRNAs在干旱脅迫下的差異表達(dá)分析

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

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

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

2.4 差異表達(dá)miRNA靶基因預(yù)測及GO功能富集分析

miRNAs調(diào)控的靶基因鑒定是了解miRNA功能的關(guān)鍵。本次試驗(yàn)干旱脅迫后葉片和根系已知miRNAs預(yù)測的靶基因分別為430個(gè)和158個(gè)，新miRNAs所預(yù)測到的靶基因則多達(dá)幾千個(gè)，結(jié)合轉(zhuǎn)錄組注釋結(jié)果[25]，發(fā)現(xiàn)預(yù)測的這些靶基因多為轉(zhuǎn)錄因子、其他的轉(zhuǎn)錄調(diào)節(jié)物質(zhì)和一些與脅迫相關(guān)的蛋白及與激素信號轉(zhuǎn)導(dǎo)相關(guān)的物質(zhì)。

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

圖3 干旱脅迫條件下白羊草差異表達(dá)miRNAs預(yù)測靶基因的GO功能富集分析Fig.3 Function analysis of target transcripts of differentially expressed miRNAs inB.ischaemum under drought stress.注：橫坐標(biāo)為Ontology分類，縱坐標(biāo)為注釋到該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富集分析

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

表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通路顯著富集的標(biāo)準(zhǔn)，q代表調(diào)整后的P值，a代表富集到某一特定通路的差異表達(dá)基因的數(shù)目，b代表富集到該通路上的所有unigene數(shù)

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驗(yàn)證

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

3 討論

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

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

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

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

4 結(jié)論

這次通過對干旱脅迫條件下白羊草葉片和根系miRNAs表達(dá)特性的研究發(fā)現(xiàn)，白羊草響應(yīng)干旱脅迫的miRNAs除了一些在其他植物中已經(jīng)被證實(shí)與干旱脅迫相關(guān)的miRNAs外，還發(fā)現(xiàn)一些與干旱脅迫相關(guān)的新miRNAs，通過對miRNAs靶基因生物學(xué)功能和代謝途徑的富集分析發(fā)現(xiàn)，植物的次生代謝產(chǎn)物尤其是黃酮類化合物合成和植物信號轉(zhuǎn)導(dǎo)途徑是白羊草響應(yīng)干旱脅迫富集的主要通路。