徐 媛，陳錦玲，陳玉梅，李璐璐，李惠敏，秦新民

(廣西師范大學 生命科學學院，廣西 桂林 541004)

0 引言

【研究意義】花生(ArachishypogaeaL.)是我國重要的油料作物，也是植物油和蛋白質的來源之一。2014年我國花生種植面積達460.39萬hm2，總產量為1 648.2萬t，成為產量最大的油料作物[1]。我國花生產區(qū)主要分布在干旱、半干旱地區(qū)。干旱是花生生產上影響最大的限制因素，據統(tǒng)計，由于干旱引起的花生減產占全國總產的20%以上[2]。干旱除引起減產外，還使花生品質下降、黃曲霉污染及病蟲害加重等[3]。因此，開展花生響應干旱的分子機制研究，在提高花生抗逆性、增加產量和提高品質方面均極為重要?！厩叭搜芯窟M展】近年來，轉錄組測序技術在花生干旱研究中得到廣泛應用，為進一步認識花生抗旱分子機制提供了重要依據。趙小波等[4]利用轉錄組測序技術對抗旱花生品種J11在干旱脅迫下轉錄因子的表達變化分析，發(fā)現236個差異表達轉錄因子。孫愛清等[5]利用 Solexa 高通量測序技術對干旱脅迫下的花生葉片 cDNA文庫進行差異基因表達譜分析，發(fā)現9個類黃酮代謝相關基因在干旱脅迫下顯著上調表達，推測類黃酮在干旱脅迫下可能發(fā)揮重要作用。萬麗云等[6]通過對10個不同的花生材料苗期進行干旱-復水試驗，篩選抗感材料進行轉錄組分析，發(fā)現抗感材料的差異表達基因主要富集在氧化磷酸化、光合作用和植物代謝途徑。胡述浩[7]研究發(fā)現花生60條保守miRNA序列，其靶基因為392個，這些靶點包括編碼轉錄因子和其他功能蛋白，如NAC轉錄因子、乙酰輔酶A羧化酶、色氨酸合成酶α鏈、葉綠體核糖核酸蛋白等，同時發(fā)現隨著高滲溶液濃度升高，miR2111、miR482的表達量明顯降低，而分別對應的靶基因表達量則顯著升高，植株的抗旱能力也隨之增加?！狙芯壳腥朦c】花生在干旱脅迫過程中，差異基因和miRNA的表達變化以及差異表達的miRNA所調控基因種類和功能尚需深入研究?！緮M解決的關鍵問題】近年來，PEG6000常用于植物模擬干旱研究[8-11]，通過對干旱脅迫下花生根系進行轉錄組測序和miRNA測序，結合轉錄組測序數據中基因的表達水平和miRNA測序中靶基因的表達水平進行統(tǒng)計分析，旨在挖掘出與抗旱相關的關鍵基因和miRNA，為花生的抗旱分子機理研究與抗旱育種提供分子依據。

1 材料與方法

1.1 材料

花生品種為桂花37，廣西農業(yè)科學院經濟作物研究所選育的珍珠豆型高油酸花生新品種。

1.2 方法

1.2.1 花生種子干旱處理 花生種子經75%酒精消毒10 min，在28℃培養(yǎng)箱黑暗條件下蒸餾水浸泡10 h后種在滅過菌的蛭石中，澆灌蒸餾水，待其發(fā)芽后澆灌霍格蘭營養(yǎng)液；用營養(yǎng)液配置15%濃度的PEG6000，待花生長至3葉一心時，選擇生長健壯、長勢一致的花生苗小心沖洗根部，將其進行15% PEG6000脅迫處理0 d、1 d、2 d和 3 d。收集處理后的花生根系，貯存在-80℃的冰箱，作為轉錄組和miRNA測序的材料。0 d、1 d、2 d、3 d處理后的材料分別以HS15-0、HS15-1、HS15-2和HS15-3命名。

1.2.2 花生總RNA的提取、文庫構建及測序 樣品總RNA的提取按照trizol試劑盒說明書進行。通過Agilent 2100 Bioanalyzer分析儀和紫外分光光度計對總RNA的質量、濃度和完整性進行檢測。4個RNA文庫(HS15-0、HS15-1、HS15-2、HS15-3)構建及測序由華大基因科技有限公司完成。

1.3 數據處理和分析

轉錄組測序和miRNA測序數據的處理分別參照文獻[12-13]進行。轉錄組數據使用過濾軟件SOAPnuke對BGISEQ-500平臺測序獲得的原始數據進行篩選,去除包含接頭的reads、未知堿基N含量大于5%的reads,以及去除低質量的reads。得到clean reads之后,以花生基因組 (ftp://ftp.ncbi.nlm.nih.gov/genomes/Arachis_duranensis)作為參考基因組，使用HISAT (http://www.ccb.jhu.edu/software/hisat)將 clean reads比對到參考基因組序列,使用Cufflinks(http://cole-trapnell-lab.github.io/cufflinks)將比對結果組裝構建轉錄本。miRNA測序將獲得的sRNA經過blast或者bowtie比對到miRBase數據庫，鑒定得到已知的miRNA。新miRNA采用華大自主開發(fā)的Mireap進行預測。

2 結果與分析

2.1 轉錄組序列的拼接與組裝

轉錄組HS15-0、HS15-1、HS15-2和HS15-3 4個文庫分別獲得43 625 768、44 792 456、43 100 476和36 638 768條Clean reads，過濾后的Reads數比對到參考基因組(花生：NCBI_Aradu1.1 ftp://ftp.ncbi.nlm.nih. gov/genomes/ Arachis_duranensis/)，比對上參考基因組的占比分別是77.19%、74.76%、69.83%和51.78%。同時，測序數據表明具有蛋白編碼潛力的轉錄本有13 387條，其中屬于已知基因的具有編碼潛力的轉錄本為11 590條，占86.57%，而屬于新基因的具有編碼潛力的轉錄本為1 797條，占13.42%。

2.2 轉錄組差異表達基因

按照FPKM法計算差異基因的表達量，當倍數變化>0時為上調表達，反之，則為下調表達。同時以log(FPKM)≥2、FDR≤0.001為條件，篩選顯著差異表達基因。在HS15-0與HS15-1比較組中，差異基因總數為14 018個，其中上調和下調差異表達的基因分別為3 836個和10 182個；HS15-1與HS15-2比較組中，差異基因總數是8 811個，其中上調和下調差異表達的基因分別有2 068個和6 743個；HS15-2與HS15-3比較組中，基因總數是5 936個，其中上調和下調差異表達的基因分別為2 264個和3 672個。

2.2.1 差異基因的GO功能 將3個比較組的差異基因進行GO功能分類，在HS15-0與HS15-1組的差異表達基因中，共有44個功能小類：其中分子功能有12類，共28 479個差異表達基因；細胞組分有14類，共20 730個差異表達基因；生物過程有18類，共24 065個差異表達基因。HS15-1與HS15-2組的差異表達基因中，共有為44小類：11類分子功能、14類細胞組分及19類生物過程，差異表達基因分別為19 032個、13 898個和16 107個。HS15-2與HS15-3組的差異表達基因中，共有42小類：10類分子功能、14類細胞組分以及18類生物過程，差異表達基因分別為12 065個、8 841個和10 631個。

2.2.2 差異基因的GO功能富集 對各個比較組中的差異表達基因富集到的GO Term進行統(tǒng)計分析發(fā)現，HS15-0與HS15-1組的差異表達基因被富集到4 160個GO Term中，其中生物過程2 377個、細胞組分545個、分子功能1 238個；HS15-1與HS15-2組的差異表達基因被富集到3 867個GO Term中，其中生物過程2 210個、細胞組分500個、分子功能1 157個；HS15-2與HS15-3組中的差異表達基因被富集到3 398個GO Term中，其中生物過程1 946個、細胞組分472個、分子功能980個。從3個功能(生物過程、細胞組分和分子功能)被富集到的條目數發(fā)現，3個比較組中參與生物過程差異表達基因的數目最多，而細胞組分的差異基因數目最少。

2.2.3 差異基因的KEGG代謝通路 對HS15-0與HS15-1、HS15-1與HS15-2、HS15-2與HS15-3 3個比較組中的差異基因進行KEGG Pathway注釋分類，共分為6個分支：細胞過程(Cellular Processes)、環(huán)境信息處理(Environmental Information Processing)、人類疾病(Human Disease)(僅限動物)、遺傳信息處理(Genetic Information Processing)、有機系統(tǒng)(Organismal Systems)和代謝(Metabolism)。3個比較組差異基因在代謝分支中的注釋的數目最多。

根據 KEGG Pathway 注釋分類結果，采用函數進行富集分析得出，3個比較組中的差異表達基因共同富集到133個代謝通路中。在HS15-0與HS15-1比較組中，植物激素信號轉導(Plant hormone signal transduction)、MAPK信號通路-植物(MAPK signaling pathway - plant)、核糖體(Ribosome)、植物病原體互作(Plant-pathogen interaction)、碳代謝(Carbon metabolism)通路中差異基因數目比較多，分別為368個、357個、341個、339個和319個。在HS15-1與HS15-2比較組中，核糖體(Ribosome)氨基酸的生物合成(Biosynthesis of amino acids)、碳代謝(Carbon metabolism)植物激素信號轉導(Plant hormone signal transduction)和RNA轉運(RNA transport)通路中的差異基因數目比較多，分別為251個、241個、207個和202個；在HS15-2與HS15-3比較組中，核糖體(Ribosome)、氨基酸的生物合成(Biosynthesis of amino acids)、碳代謝(Carbon metabolism)、植物激素信號轉導(Plant hormone signal transduction)、RNA轉運遺傳信息處理(RNA transport Genetic Information Processing)通路中的差異基因數目比較多，分別為251個、241個、207個和202個；核糖體(Ribosome)、氨基酸的生物合成(Biosynthesis of amino acids)、碳代謝(Carbon metabolism)、RNA轉運、植物病原菌互作(Plant-pathogen interaction)通路中的基因數目最多，分別為251個、172個、172個、144個和143個。

2.2.4 持續(xù)上調和下調表達基因 在3個比較組中共篩選出13個持續(xù)上調的差異基因，其中7個為HSP分子伴侶家族基因、2個RING E3泛素連接酶、1個為谷胱甘肽-S-轉移酶、1個柱頭特異性STIG1樣蛋白1、1個免疫球蛋白結合蛋白和1個 DnaJ同源B亞家族成員13。同時篩選到13個持續(xù)下調的差異基因。其表達量及注釋見表1。

表1 持續(xù)上調和下調基因的表達量及注釋Table 1 Expression and annotation of continuously up- and down-regulated gene

2.3 樣品miRNA的測序

對HS15-0、HS15-1、HS15-2 和HS15-3 4個樣品進行的miRNA測序，分別獲得30 182 740條、29 044 095條、29 702 900條和28 200 685條 Raw reads，去除低質量、沒有3’接頭序列、5’接頭污染、小于18nt和包含poly A的序列，分別獲得28 005 836條、26 862 969條、26 668 939條和25 970 037 條clean reads。將4個樣本中所有的 clean reads 和已知的miRNA數據庫(miRBase、Rfam、siRNA、piRNA、snoRNA等)進行比對，能比對上的 reads 分別有21 804 361條、21 797 056條、21 192 973條和15 074 860條，分別占 clean reads 的 77.86%、81.14%、79.47%和 58.05%。

2.3.1 已知 miRNAs的鑒定 將HS15-0、HS15-1、HS15-2和 HS15-3 4個樣本文庫中比對上基因組的sRNA比對到miRbase數據庫中，分別鑒定到 28個、27個、27個和 26個已知miRNA，其中4個文庫中共有的已知 miRNA有25個。

2.3.2 新miRNA的鑒定 在樣本HS15-0、HS15-1、HS15-2和 HS15-3中分別鑒定出116個、113個、110個和106個新的miRNA，其中有105個新 miRNAs在4個樣本中為共表達。

2.4 差異表達miRNA

2.4.1 差異表達miRNA的篩選 在HS15-0與HS15-1比較組中，差異表達的miRNA有23個，其中有7個為上調差異miRNA，16個為下調差異miRNA；在HS15-1與HS15-2比較組中，差異表達的miRNA有36個，其中11個為上調差異miRNA，25個為下調差異miRNA；在HS15-2與HS15-3比較組中，差異表達的miRNA有57個，其中7個為上調差異miRNA，50個為下調差異miRNA。在3個比較組中，下調差異表達的miRNA數目多于上調差異表達的miRNA。

2.4.2 差異表達miRNA的靶基因 利用Target Finder等軟件對差異表達miRNA進行靶基因預測。HS15-0與HS15-1比較組中15個miRNAs為差異表達，共預測到 392個靶基因，其中5個已知 miRNA調控 330個靶基因，10個新 miRNA調控62個靶基因；HS15-1與HS15-2比較組中有25個差異表達 miRNAs，預測到的靶基因數為913個，其中6個已知miRNA調控440個靶基因，19個新miRNA調控473個靶基因；HS15-2與HS15-3比較組中的 44個miRNAs呈現出差異表達，預測的靶基因為1 718個，其中8個已知miRNA調控的的靶基因數為1 007個，而36個新 miRNA僅調控 711個靶基因。

2.4.3 差異表達miRNA靶基因的GO和KEGG富集 將3個比較組的差異表達 miRNA的靶基因進行GO功能分類。在HS15-0與HS15-1比較組差異表達miRNA靶基因中，生物過程、細胞組分和分子功能三大類所占的靶基因數目分別為495個、429個和246個；HS15-1與HS15-2比較組的靶基因數目分別是1 089個、860個和504個；HS15-2與HS15-3比較組的靶基因數目分別是1 479個、1 324個和766個。數據分析發(fā)現，HS15-2與HS15-3比較組中差異表達miRNA的靶基因數最多，HS15-0與HS15-1比較組最少。在GO分類中，3個比較組的功能分類模式相似。在生物過程分類中，3個比較組注釋到差異表達miRNA的靶基因數主要分布在細胞過程(cellular process)和代謝過程(metabolic process)類別中；在細胞組分分類中，3個比較組中的靶基因主要集中在細胞(cell)、膜(membrane)、膜部分(membrane part)和細胞器(organelle)4個小類別中；在分子功能類別中，3個比較組中的差異表達miRNA靶基因主要分布在結合(binding)和催化活性(catalytic activity)類型中。該分類結果與轉錄組測序得到的差異基因注釋分類結果一致。在HS15-0與HS15-1比較組中，共富集到76條代謝通路，其中代謝途徑(Metabolic pathways)、次生代謝產物的生物合成(Biosynthesis of secondary metabolites)、晝夜節(jié)律-植物(Circadian rhythm - plant)和類黃酮生物合成(Flavonoid biosynthesis)通路中差異miRNA靶基因的數目比較多，分別為98個、80個、27個和21個。在HS15-1與HS15-2比較組中，共富集到91條代謝通路，其中代謝途徑(Metabolic pathways)、次生代謝產物的生物合成(Biosynthesis of secondary metabolites)、內質網中的蛋白質加工(Protein processing in endoplasmic reticulum)和剪接體(Spliceosome)通路中差異miRNA靶基因數目比較多，分別為123個、53個、45個和44個；在HS15-2與HS15-3比較組中，共富集到105條代謝通路，其中代謝途徑(Metabolic pathways)、次生代謝產物的生物合成(Biosynthesis of secondary metabolites)、植物激素信號轉導(Plant hormone signal transduction)和胞吞作用(Endocytosis)通路中差異miRNA靶基因數目比較多，分別為205個、104個、69個和65個。數據分析發(fā)現，在花生干旱脅迫下，有較多的靶基因參與調控代謝和次生代謝產物的生物合成途徑。

2.5 與花生干旱相關的miRNA及靶基因預測

通過miRNA的靶基因功能注釋分析發(fā)現，有13條miRNA的靶基因注釋為與干旱相關的轉錄因子或基因，其中9條為持續(xù)下調miRNA，涉及35個基因，源異型域-亮氨酸拉鏈蛋白(homeobox-leucine zipper protein, HD-Zip)編碼基因1個、MYB類LHY基因2個，WRKY轉錄因子編碼基因3個、谷胱甘肽S-轉移酶(glutathione S-transferase，GST)基因2個、 SPL轉錄因子(Squamosa promoter binding protein-like)基因20個、 CCCH型鋅指蛋白(zinc finger CCCH domain protein)基因7個(表2)。ahy-miR156a分別參與WRKY轉錄因子、谷胱甘肽S-轉移酶和SPL轉錄因子基因的調控。

表2 花生干旱相關的靶基因功能注釋及miRNATable 2 Functional annotation and miRNA of target gene related to peanut under drought stress

續(xù)表2

3 討論

3.1 抗旱功能相關基因

植物熱激蛋白(HSPs)在植物減輕逆境脅迫引起的傷害方面發(fā)揮作用[14]。王利彬等[15]發(fā)現轉HSP70煙草的耐旱性和耐熱性得到顯著提高。WANG等[16]在研究大麥干旱脅迫差異基因時發(fā)現HSP90蛋白在不同處理上出現差異表達。西瓜在水分脅迫時熱激蛋白也出現上調表達[17]。此外，在小麥研究中也發(fā)現干旱脅迫會誘導抗旱小麥材料中HSP60和HSP90表達上調[18]。最近有研究表明，泛素連接酶E3在植物的抗旱功能中發(fā)揮至關重要的作用[19-20]，持續(xù)上調的RING finger E3泛素連接酶可通過調控脫落酸(ABA)信號通路參與干旱脅迫響應[21]。XERICO是一種RING-H2E3泛素連接酶，ZENG等[22]發(fā)現，過表達XERICO能增強ABA的生物合成，從而提高擬南芥的抗旱性。谷胱甘肽-S-轉移酶(GSTs)是一類重要的多功能超家族酶，GSTs可通過親電取代反應、降低毒害作用和過氧化酶清除功能，增加植物適應各種逆境脅迫的能力[23-24]。如過表達PjGSTU1基因可提高煙草的抗旱能力[25]。近年來，DnaJ蛋白在植物干旱脅迫的研究方面取得了較大進展[26]。SCHAFLITNER等[27]鑒定出6 個受干旱脅迫誘導上調表達的DnaJ蛋白基因。RODRIGUES等[28]從甘蔗中鑒定出 3個DnaJ蛋白基因，并發(fā)現其中有 1個在耐旱品種中下調表達，而在不耐旱品種中則上調表達。THUMMA等[29]鑒定出1個桉樹DnaJ蛋白，發(fā)現在干旱脅迫條件下，其3個等位基因的表達差異極其顯著。研究結果表明，在HS15-0與HS15-1、HS15-1與HS15-2和HS15-2與HS15-3 3個比較組中共篩選出13個持續(xù)上調表達基因，分別為7個HSP分子伴侶家族基因、2個RING E3泛素連接酶、1個GSTs，1個柱頭特異性STIG1樣蛋白1、1個免疫球結合蛋白和1個 DnaJ同源B亞家族成員13。其中HSPs分子伴侶家族基因、 RING E3泛素連接酶、 GSTs和DnaJ基因很可能與花生應激干旱脅迫相關。而柱頭特異性STIG1樣蛋白1和免疫球結合蛋白基因也為持續(xù)上調表達，但目前尚未見這2個基因與植物干旱脅迫相關的研究報道。此外，13個持續(xù)下調表達的基因其功能注釋與花生干旱脅迫的關系也不明確，有待深入研究。

3.2 花生根系抗旱功能相關miRNA

HD-Zip是植物特有的轉錄因子，其可通過與相關激素基因及其下游基因互作增強植物耐逆性[30]。研究中miRNA中novel_mir70的靶基因為同源異型域-亮氨酸拉鏈蛋白基因。隨著脅迫時間延長，novel_mir70表達持續(xù)減少，通過負調控機制使得靶基因表達增強，有利于植物抗旱能力的提高。MYB-related transcription factor LHY屬于生物鐘節(jié)律基因，其主要功能不僅調控生物晝夜節(jié)律，還可以通過識別脅迫信號通路基因上游的調控元件，調控植物對非生物脅迫作出反應[31]。PATADE等[32]在甘蔗的PEG脅迫中發(fā)現MYB和ahy-miR159的表達水平彼此相反。

ahy-miRNA159可以通過控制ABI3的表達來調控 MYB33轉錄因子的含量，從而參與ABA信號通路響應干旱[33]。研究中靶向MYB-related transcription factor LHY的2條miRNA ahy-miR156b-3p和ahy-miR159均呈持續(xù)下降的表達模式，可能通過負調控來增強花生的抗旱性。WKRY基因家族是高等植物中最大的轉錄因子家族之一，其蛋白結構基本含有 1～2個 WRKY結構域，參與非生物逆境應答如干旱、澇漬、高鹽、高溫、低溫、損傷等，如BhWRKY1 與BhGolS1(半乳糖醇合成酶)啟動子的 W-box 結合，轉錄激活BhGolS1的表達，可以提高擬南芥的耐旱性[34]。研究發(fā)現ahy-miR156a、novel_mir103和 novel_mir115靶向4個WKRY轉錄因子基因，其中ahy-miR156a和novel_mir103在脅迫的1～3 d過程中表達量也表現為逐步下降模式，使得靶基因WKRY轉錄因子基因高表達，從而提高花生的抗旱能力。谷胱甘肽S-轉移酶(GSTs)為細胞內多功能蛋白，在植物體內參與初生代謝、次生代謝、解毒以及非生物脅迫過程[35-36]。研究發(fā)現2個miRNA(novel_mir83和ahy-miR156a)的靶基因為谷胱甘肽S-轉移酶基因，其表達量在脅迫過程中呈持續(xù)下降趨勢。 該結果與轉錄組中GSTs基因表達持續(xù)上調一致，推測novel_mir83和ahy-miR156a通過負調控作用，增加GSTs基因的表達，增強花生的抗旱能力。

鱗狀啟動子結合蛋白樣(Squamosa promoter binding protein-like，SPL)是植物中特有的轉錄因子。已有研究證明SPL可以增強植物在干旱脅迫下的耐受性[37]。在植物生長發(fā)育過程中miR156與SPL的表達水平通常呈負相關，miR156可以通過mRNA剪切或翻譯抑制來調控SPL的表達[38]。研究共篩選到2個下調表達的miR156家族蛋白(ahy-miR156a和ahy-miR156b-5p)，所靶向的20個靶基因功能均注釋為SPL轉錄因子。由此推斷，干旱脅迫下ahy-miR156a、ahy-miR156b-5p通過下調，使其靶基因表達水平增強來提高自身的抗逆作用。CCCH型鋅指蛋白(zinc finger CCCH domain protein)是一類重要的調控因子，在動植物生長發(fā)育以及對外界逆境的反應過程中發(fā)揮著重要的作用。徐倩倩[39]研究發(fā)現，PEG可以誘導玉米中的ZmC3H54的表達，干旱處理下過表達的 Zm C3H54 轉基因水稻植株抗旱性增強。本研究中發(fā)現1個下調的ahy-miR167-5p，其7個靶基因功能注釋為CCCH型鋅指蛋白，由此推測干旱脅迫下ahy-miR167-5p可能負調控CCCH型鋅指蛋白，從而減少逆境對植物的傷害。

以上幾個家族的miRNA可能通過作用于各自的靶基因參與到花生的干旱脅迫中，但這些miRNAs參與花生干旱脅迫應答機制尚待進一步研究。

4 結論

研究采用15% PEG6000對花生幼苗進行了干旱脅迫處理，通過對材料的轉錄組和miNRA測序，在轉錄組測序數據中發(fā)現了13個顯著持續(xù)上調的差異基因和13個顯著下調的差異基因；在miRNA測序數據中篩選到8個持續(xù)下調表達的miRNA，分別為ahy-mir156-3p、 ahy-miR159、ahy-miR167-5p、ahy-miR156a、ahy-miR156b-5p、novel_mir103、novel_mir83和novel_mir70。這些差異表達基因和miRNA參與了花生的干旱脅迫過程。