利用高通量測序分析白花泡桐鹽脅迫相關microRNAs

王園龍， 曹 林， 鄧敏捷， 馬一平，趙振利， 牛蘇燕， 王曉丹， 范國強

(1.河南農業(yè)大學泡桐研究所，河南 鄭州 450002； 2.北京師范大學地理與遙感科學學院，北京 100875)

?

利用高通量測序分析白花泡桐鹽脅迫相關microRNAs

王園龍1， 曹 林1， 鄧敏捷1， 馬一平2，趙振利1， 牛蘇燕1， 王曉丹1， 范國強1

(1.河南農業(yè)大學泡桐研究所，河南 鄭州 450002； 2.北京師范大學地理與遙感科學學院，北京 100875)

為了鑒定和研究白花泡桐中與鹽脅迫相關的microRNAs(miRNAs)，構建了對照組PF2和處理組PF2-S4 2個小RNA文庫。采用高通量測序技術總共鑒定出來了33個已知miRNA和80個新miRNA，在這些miRNAs中分別有27個(3個已知miRNA和24個新miRNA)差異顯著表達和有45個(13個已知miRNA和32個新miRNA)差異極顯著表達. 最后采用qRT-PCR技術驗證了5個miRNA。

白花泡桐；microRNAs；鹽脅迫；高通量測序

木本植物不僅在水土保持和維持自然界碳循環(huán)中起關鍵作用，而且為工業(yè)、造紙業(yè)等民生行業(yè)提供重要的原材料[1]。但是它們經常遭受到環(huán)境中包括鹽脅迫在內的多種生物和非生物脅迫而導致減產甚至死亡。當今土壤鹽漬化不單是干旱半干旱地區(qū)的問題，已演變成為世界范圍內影響農林業(yè)發(fā)展的重要因素。全球范圍內超過10%的陸地正在遭受到不同程度的鹽漬化侵害而且還呈現不斷加重的趨勢[2]。所以對植物的抗鹽脅迫研究及育種已經成為擺在植物新品種培育科研工作者面前的一道急需解決的難題。泡桐因其生長迅速，材性優(yōu)良，適合農桐間作而被廣泛種植在世界范圍內250萬hm2的農田中[3,4]。泡桐產業(yè)的發(fā)展可以為帶來良好的經濟效益和生態(tài)效益，但是土壤鹽漬化約束著泡桐產業(yè)的發(fā)展。所以在泡桐的抗鹽脅迫基礎研究中首先鑒定出泡桐抗鹽脅迫相關基因，將有助于弄清楚泡桐的抗鹽機制，并且有針對性地進行泡桐的新品種培育。近年來，有關泡桐抗叢枝病[5-7]，抗干旱[3, 6,8,9]的轉錄組研究取得了一定的進展，找到了與泡桐叢枝病發(fā)病及干旱相關的基因，但是至今為止還沒有關于泡桐抗鹽脅迫等相關分子領域的研究報道。microRNAs(miRNAs)是一類內源性非編碼sRNAs，在轉錄后水平上起到負調控基因表達的作用[10]。植物中，miRNAs參與包括生長、發(fā)育、新陳代謝和抗脅迫等多種生物學過程[11]。目前，大量響應植物生長發(fā)育以及與脅迫相關的miRNAs被鑒定出來[10,12-14]。近年來, 范國強等[4,15, 16]采用高通量測序技術已鑒定了二、四倍體泡桐種的差異表達miRNAs并且預測分析了它們的功能。然而，至今還未有關于鹽脅迫下泡桐的miRNAs的研究報道。本研究構建了處理組和對照組2個sRNA文庫來鑒定分析鹽脅迫下白花泡桐的miRNAs，并采用qRT-PCR方法來驗證鹽脅迫下白花泡桐差異表達miRNAs的準確性。以期為揭示miRNAs調控鹽脅迫條件下白花泡桐體內基因表達的分子機制，闡明鹽脅迫的分子機理奠定理論基礎。

1 材料與方法

1.1 試驗材料

植物材料為河南農業(yè)大學林木生物技術實驗室培養(yǎng)的二倍體白花泡桐組織培養(yǎng)苗。在溫度(25±2)℃，光強130 μmol·m-2·s-1，光照時間16 h·d-1的條件下進行愈傷組織培養(yǎng)30 d。而后移植到室外培養(yǎng)30 d后選擇生長一致植株單株栽植于裝有等量普通園土的營養(yǎng)缽中，正常生長50 d后進行鹽脅迫處理。

以土壤含鹽量(NaCl質量/干土質量×100%)分別為0，0.4%設置處理試驗且命名為PF2和PF2-S4。每個處理組設置3個生物學重復。具體試驗方法如下：(1)分別稱取每缽中所需NaCl的質量，平均分成3份；(2)取1份完全溶于水后施入園土中，缽底墊有托盤，將滲水重新倒入缽中，避免鹽分損失，每3 d進行1次，共計3次；(3)每次處理時，對照組PF2只加入等量水；(4)鹽分添加完畢后，每2 d澆1次水，維持75%相對土壤含水量下培養(yǎng)20 d。處理結束后，分別摘取幼苗新梢頂端第二對葉片，立即用液氮速凍后保存于-86 ℃冰箱中備用。

1.2 白花泡桐sRNA文庫的構建及測序

本試驗采用pBIOZOL(杭州博日)植物提取試劑提取白花泡桐總RNA，而后使用NanoDrop 2000超微量分光光度計和Agilent 2100生物分析儀對總RNA進行質量檢測。本試驗的PF2和PF2-S4文庫構建及測序交由深圳華大基因研究院采用高通量測序平臺Illumina HiSeqTM2000完成。

1.3 白花泡桐sRNA信息分析

高通量測序所得的50 nt序列，通過過濾得到可信的目標序列，對這些序列的質量、長度及樣品間公共序列進行統(tǒng)計。通過目標序列分類注釋，可以獲得樣品中包含的各組分及表達量信息．對未注釋的小RNA片段來進行新miRNA預測．隨后對已知miRNA和新miRNA進行差異分析。

1.4 已知miRNAs的鑒定

利用SOAP軟件將長度在18～30 nt的高質量片段(clean reads)匹配到白花泡桐轉錄組(http://trace.ncbi.nlm.nih.gov/Traces/sra/)。隨后完美匹配序列用Blast軟件(http://www.ncbi.nlm.nih.gov/staff/tao/URLAPI/blastall/)與NCBI GeneBank (http://www.ncbi.nlm.nih.gov/)和Rfam 11.0 (http://rfam.janelia.org/)數據庫比對，篩選和去除同rRNAs, scRNAs, snoRNAs, snRNAs, tRNA等相關的序列。將這2文庫中剩下的序列比對到miRBase 21.0(http://www.mirbase.org/ftp.shtml)中的已知植物miRNAs上鑒定白花泡桐已知miRNAs。

1.5 新miRNAs的鑒定

將未注釋的sRNA用于白花泡桐新miRNA的預測。使用mireap軟件 (http://sourceforge.net/projects/mireap/files/mireap/)預測2個文庫中新miRNA，繪制新miRNA的二級結構圖。如果序列存在Meyers等人描述的標準：序列的長度范圍在18～25 nt之間，最低自由能(minimum free energy，MFE)不高于-18 kcal·mol-1，miRNA與miRNA*之間的對稱度為7 nt、最長距離為100 nt、堿基配對數為10 nt，可以看作白花泡桐中新miRNA[17]。

1.6 2個文庫中miRNA的差異表達分析

根據每個miRNA分別在PF2和PF2-S4庫中的相對表達豐度確定miRNA的差異表達，如果miRNA的fold change值 ≥1或者≤-1且P-values≤0.05，定義其在2庫中差異顯著表達；如果miRNA的fold change值≥1或者≤-1且P-values≤0.01，則定義其在2庫中差異極顯著表達。

具體計算方法如下：

(1)miRNA在各自庫中表達量歸一化值=miRNA表達量/clean reads總數*1 000 000，如果結果是0，將其變更為0.01。

(2)某一miRNA的fold change值=log2(PF2-S4庫中miRNA表達量的歸一化值/PF2庫中miRNA表達量的歸一化值)。

(3)P-values計算公式。

其中：N1代表PF2文庫序列總數；N2代表PF2-S4文庫序列總數；x代表某一miRNA在PF2文庫的表達量；y代表某一miRNA在PF2-S4文庫的表達量。

1.7qRT-PCR驗證差異表達miRNAs

為了驗證鹽脅迫下白花泡桐中差異表達miRNAs及其靶基因的存在和表達模式，隨機選擇了5個差異表達miRNA進行qRT-PCR驗證。將生長狀況良好且生長一致的白花泡桐組培苗進行0.4%濃度的鹽脅迫處理，處理方法同1.1，對對照組和處理組葉片進行總RNA的提取，且總RNA提取與檢測方法同1.2。

PCR操作程序參照SuperScriptIIIplatinumSYBRGreenone-stepqRT-PCRkit(Invitrogen，USA)試劑盒說明進行。使用CFX96TMReal-TimePCRSystem(Bio-Rad,Hercules,CA)進行實驗。其反應體系總體積為20μL，其中，10μLSSIIIandplatinumMix、0.4μL5×buffer、0.4μLReverseprimer(10μmol·L-1)、0.4μLForwardprimer(10μmol·L-1)、0.4μLReverseuniverseprimer(10μmol·L-1)、500ngRNA，最后，用DEPC處理過的無菌水補至20μL。擴增程序為：50 ℃，3min；95 ℃，5min；40循環(huán)(95 ℃，15s；55 ℃，30s); 40 ℃, 10min。每個樣品進行3次技術重復。數據處理使用2-△△CT方法計算相對表達量，以內參U6rRNA進行標準化，將白花泡桐(對照樣品PF2)的表達量作為“1”，計算出處理樣品的表達量，并作圖分析。

2 結果與分析

2.1 白花泡桐sRNA文庫的構建及測序

根據鄧敏捷等人的不同質量濃度鹽處理下白花泡桐生理響應研究結果，選用0.4%土壤含鹽量作為測序處理組處理質量濃度[18]。利用Solexa/Illumina測序技術分別構建PF2和PF2-S4 2個小RNA的cDNA文庫。在PF2和PF2-S4庫中，分別有12 248 047條和11 942 362條rawreads。然后分別去除測序質量較低的序列、帶有5′接頭污染的序列、無插入片段的序列、沒有3′接頭序列、包含polyA的序列、小于18nt的片段，共得到12 060 908條和11 855 951條cleanreads。在2個樣品間共有序列種數1 176 068條，占庫中所有序列百分數為16.53%；共有序列總量17 053 689條，占庫中所有序列百分數為71.30%。測序結果表明2個樣品間序列種類上的差異比較大，但公共部分的序列表達比較集中，說明了2個樣品在測序整體上的一致性是比較好的。

對PF2和PF2-S4文庫中sRNA進行長度作圖分析(圖1)，2個文庫中24nt長度的sRNA峰值最高，分別占到sRNA總量的43.70%(PF2)和41.76%(PF2-S4)；其次是21nt長度序列，分別占到sRNA總量的22.06%(PF2)和24.84%(PF2-S4)。相似的結果也出現在毛泡桐[4]和南方泡桐[15]等其他泡桐種的測序結果中。同時長度處于20～24nt間的sRNA總數分別占其sRNA總量的92.96%(PF2)和91.88%(PF2-S4)，表明了2個樣品的sRNA文庫中均富含miRNA。

圖1 對照組(PF2)和處理組(PF2-S4)白花泡桐樣品中小RNAs的長度分布

2.2 sRNA分類注釋

將sRNA匹配到Genbank數據庫和Rfam數據庫中進行sRNA的分類注釋確保每一個unique sRNA都有唯一的注釋。按照Genbank數據庫的優(yōu)先級大于Rfam數據庫且rRNAetc>known miRNA> piRNA>repeat>exon>intron的優(yōu)先級順序將sRNA遍歷，沒有比上任何注釋信息的sRNA用unannote表示。其中分類注釋結果中的rRNA總量可以作為一個樣品的質控標準，植物樣品中的rRNA總量所占比例應低于60%(表1)。PF2和PF2-S4中rRNA總量所占比例分別為3.43%和4.22%，遠遠低于其質控標準。說明2個樣品的sRNA是測序所得的真實可信的數據。

表1 白花泡桐對照組(PF2)與處理組(PF2-S4)小RNA分類

2.3 白花泡桐中已知miRNAs的鑒定

將PF2和PF2-S4庫中完全匹配到白花泡桐轉錄組的序列與miRBase 21.0數據庫成熟的miRNA序列進行同源性比對，鑒定出白花泡桐已知miRNA序列33條，其中對照組(PF2)和處理組(PF2-S4)中均特有6條，2庫共有21條。鑒定出的33個已知miRNA中，有20個miRNAs分別屬于12個家族，13個miRNAs的家族未知。平均前體長度和平均最低自由能分別為155 nt和-58.2 kcal·mol-1。MIR156家族、MIR166家族、MIR167-1家族、MIR168家族和MIR169-1家族在2庫中均有較高的表達量。不同的家族表達量差異較大，而且即使是同一家族中的不同成員，其表達量也有較大差異。比如MIR167-1家族中pfo-miR167d和pfo-miR167e在2 個文庫中的表達量明顯低于其他3個家族成員。結果中有趣的是已知miRNA中只有20 nt，21 nt和22 nt這3個序列長度，長度為21 nt的miRNA最多，共有21個，其次有9個22 nt長的miRNA。這個結果與其他植物的研究結果一致，即長度為21 nt 的miRNAs 表達豐度最高[19]，這也說明了鑒定出的白花泡桐中和與鹽脅迫相關的已知miRNAs的準確性。

2.4 白花泡桐中新miRNAs的鑒定

如表1所示PF2庫中沒有注釋(unannote)的sRNA有4 063 862種，共計8 899 321條；在PF2-S4庫中沒有注釋的sRNA有4 008 628種，共計8 346 779條。利用Mireap軟件對2個庫中unannote序列進行miRNA二級結構等特征的篩選，共獲得新miRNA序列80條，其中PF2庫中特有序列34條，PF2-S4庫中特有序列29條，PF2庫和PF2-S4庫中共有序列17條。新miRNA的序列長度在20 nt到23 nt之間，其中長度為21 nt的新miRNA數量最多，有37條。平均前體長度和平均最低自由能分別為156 nt和-47.7 kcal·mol-1。與預測出來的2個庫中已知miRNAs的表達量相比，新miRNAs的表達量總體上呈現較低的水平。

2.5 PF2和PF2-S4文庫中miRNAs的差異表達分析

對白花泡桐miRNA在PF2庫和PF2-S4庫中表達量進行差異分析。與對照組PF2庫相比，在處理組PF2-S4庫中有27個miRNAs(3個已知miRNA和24個新miRNA)差異顯著表達；45個miRNAs(13個已知miRNA和32個新miRNA)差異極顯著表達(表2)。這些顯著差異和極顯著差異表達的miRNAs可能在白花泡桐的抗鹽脅迫中發(fā)揮重要的調控作用。

2.6 qRT-PCR驗證差異表達miRNA

為了進一步驗證測序結果的準確性，本研究隨機選取了5個差異表達的脅迫相關miRNAs進行實時熒光定量PCR驗證，分析鹽脅迫相關miRNAs的相對表達差異(圖2)。結果顯示，與測序結果相比，除pfo-mir890之外，pfo-miR159、pfo-miR167d、pfo-miR5021a、pfo-mir5301的表達量變化與測序的結果完全相同，說明測序結果的可靠性。

表2 PF2和PF2-S4文庫中miRNA的差異表達分析

續(xù)表 Continuing table

注： * ：差異顯著；**：差異極顯著。

Note：*：significant difference； **：more significant difference.

圖2 差異表達miRNAs的qRT-PCR驗證

3 結論與討論

本試驗采用高通量測序技術對PF2和PF2-S4 2個sRNAs文庫中的已知miRNA和新miRNA進行鑒定。通過對2個文庫的小RNA長度分布圖進行分析得出了2樣品中均富含miRNA，說明測序結果是值得信賴的；盡管有一些sRNA可能是由降解途徑產生，但PF2和PF2-S4的長度分布沒有明顯的差異，說明了sRNAs的長度分布可能獨立于鹽脅迫而被植物其他機制控制[20]。

根據PF2和PF2-S4中miRNAs的表達分析，找到了72個與鹽脅迫相關的miRNA。其中，一些miRNAs在白花泡桐正常的生長發(fā)育中不容易檢測到，但在受到鹽脅迫時會受到特異性表達而被鑒定出來；相反，一些miRNA在白花泡桐正常的生長發(fā)育中可以特異表達，但受到鹽脅迫后會被抑制從而在處理后不能被鑒定出來；還有一些miRNAs在受到鹽脅迫前后發(fā)生差異表達，這些差異顯著表達和差異極顯著表達的miRNA通過上調或者下調其表達量來響應鹽脅迫。對鑒定出來的miRNAs進行差異表達分析，發(fā)現分別在PF2和PF2-S4庫中共12個特有的已知miRNA均差異顯著表達或差異極顯著表達；2個庫中共63個特有的新miRNA中，有52個miRNAs差異顯著表達或差異極顯著表達，這些miRNAs很可能在泡桐抗鹽脅迫相關的過程中發(fā)揮重要的調控作用。

總之，本試驗是國內外首次基于泡桐轉錄組對白花泡桐抗鹽脅迫miRNA進行研究，采用高通量測序分析方法鑒定出了白花泡桐與鹽脅迫相關的miRNA。然后采用qRT-PCR方法驗證了白花泡桐miRNAs的存在。通過miRNA在2個小RNAs文庫的差異表達情況，找到與鹽脅迫相關的miRNA，有關白花泡桐抗鹽脅迫的miRNA研究將為研究木本植物的耐鹽機制提供直接和有效的參考。

[1] ATES S, NI Y, AKGUL M, et al. Characterization and evaluation ofPaulowniaelongotaas a raw material for paper production [J]. African Journal of Biotechnology, 2008, 7(22): 4153-4158.

[2] 孫文喜, 孫文鑫, 周大旺, 等. 植物抗鹽性生理生化研究進展及提高植物抗鹽性對策 [J]. 農業(yè)科技通訊, 2013 (6): 174-178.

[3] DONG Y, FAN G, DENG M, et al. Genome-wide expression profiling of the transcriptomes of fourPaulowniatomentosaaccessions in response to drought [J]. Genomics, 2014, 104(4): 295-305.

[4] FAN G, ZHAI X, NIU S, et al. Dynamic expression of novel and conserved microRNAs and their targets in diploid and tetraploid ofPaulowniatomentosa[J]. Biochimie, 2014, 102: 68-77.

[5] CAO X, FAN G, DENG M, et al. Identification of Genes Related toPaulowniaWitches’ Broom by AFLP and MSAP [J]. International Journal of Molecular Sciences, 2014, 15(8): 14669-14683.

[6] CAO X, FAN G, ZHAO Z, et al. Morphological changes ofPaulowniaseedlings infected phytoplasmas reveal the genes associated with witches' broom through AFLP and MSAP [J]. PloS one, 2014, 9(11): 533-542.

[7] LIU R, DONG Y, FAN G, et al. Discovery of Genes Related to Witches Broom Disease inPaulowniatomentosaxPaulowniafortuneiby a De Novo Assembled Transcriptome [J]. Plos One, 2013, 8(11):238-247.

[8] DONG Y, FAN G, ZHAO Z, et al. Transcriptome Expression Profiling in Response to Drought Stress inPaulowniaaustralis[J]. International Journal of Molecular Sciences, 2014, 15(3): 4583-4607.

[9] CUI L, SHAN J, SHIMIN G, et al. The miR156-SPL9-DFR pathway coordinates the relationship between development and abiotic stress tolerance in plants [J]. Plant Journal, 2014, 80(6): 1108-1117.

[10]LI B, DUAN H, LI J, et al. Global identification of miRNAs and targets inPopuluseuphraticaunder salt stress[J]. Plant Molecular Biology, 2013, 81(6): 525-539.

[11]SHUAI P, LIANG D, ZHANG Z, et al. Identification of drought-responsive and novelPopulustrichocarpamicroRNAs by high-throughput sequencing and their targets using degradome analysis [J]. Bmc Genomics, 2013, 14:233-241.

[12]BARRERA-FIGUEROA B E, GAO L, WU Z, et al. High throughput sequencing reveals novel and abiotic stress-regulated microRNAs in the inflorescences of rice [J]. Bmc Plant Biology, 2012, 12(4).132-145.

[13]HAN J, XIE H, SUN Q, et al. Bioinformatic identification and experimental validation of miRNAs from foxtail millet (Setariaitalica) [J]. Gene, 2014, 546(2): 367-377.

[14]JIA L, ZHANG D, QI X, et al. Identification of the Conserved and Novel miRNAs in Mulberry by-High-Throughput Sequencing [J]. Plos One, 2014, 9(8):409-417.

[15]NIU S, FAN G, XU E, et al. Transcriptome/Degradome-Wide Discovery of MicroRNAs and Transcript Targets in TwoPaulowniaaustralisGenotypes [J]. Plos One, 2014, 9(9):736-745.

[16]NIU S, FAN G, ZHAO Z, et al. High-throughput sequencing and degradome analysis reveal microRNA differential expression profiles and their targets inPaulowniafortunei[J]. Plant Cell Tissue and Organ Culture, 2014, 119(3): 457-468.

[17]GERMAN M A, PILLAY M, JEONG D H, et al. Global identification of microRNA-target RNA pairs by parallel analysis of RNA ends [J]. Nature Biotechnology, 2008, 26(8): 941-946.

[18]鄧敏捷, 張曉申, 范國強, 等. 四倍體南方泡桐和四倍體白花泡桐對鹽脅迫的生理響應研究 [J]. 河南農業(yè)大學學報, 2013, 47(6): 698-702.

[19]CHEN X. Small RNAs and Their Roles in Plant Development [J]. Annual Review of Cell and Developmental Biology, 2009,25: 21-44.

[20]HACKENBERG M, GUSTAFSON P, LANGRIDGE P, et al. Differential expression of microRNAs and other small RNAs in barley between water and drought conditions [J]. Plant Biotechnology Journal, 2015, 13(1): 2-13.

(責任編輯：蔣國良)

Analysis of salt stress-responsive microRNAs inPaulowniafortuneiby high-throughput sequencing

WANG Yuanlong1, CAO Lin1, DENG Minjie1, MA Yiping2,ZHAO Zhenli1,NIU Suyan1, WANG Xiaodan1, FAN Guoqiang1

(1.Institute of Paulownia, Henan Agricultural University, Zhengzhou 450002, China; 2.School of Geography,Beijing Normal University,Beijing 100875,China)

To identify and investigate salt-related microRNAs(miRNAs) inPaulowniafortunei(P.fortunei), we constructed and sequenced two small RNA libraries (PF2 and PF2-S4). A total of 33 conserved and 80 novel miRNAs were identified in the two libraries by high-throughput sequencing. Among these miRNAs, 27 miRNAs (3 conserved and 24 novel) were differentially expressed and 45 miRNAs (13 conserved and 32 novel) were significantly differentially expressed in response to salt stress inP.fortunei, respectively. Finally, 5 of the miRNAs were confirmed by qRT-PCR.

Paulowniafortunei; microRNAs; salt stress; high-throughput sequecning

2015-01-20

王園龍(1989-)，男，河南濮陽人，碩士研究生，主要從事泡桐生物技術方面的研究。

范國強(1964-)，男，河南禹州人，教授，博士生導師。

1000-2340(2015)04-0461-07

S792.43

A