旱柳轉(zhuǎn)錄組測序及生物學(xué)分析

李敏　郭聰　李玉娟　馮新民　王瑩　張健　談峰

摘要：對鹽敏感旱柳沿江柳和耐鹽旱柳9901的雜交F代根部進(jìn)行了RNA-Seq測序和分析，共獲取了107950條Unigene，平均長度為1076.96bp。通過與COG、GO等8個(gè)數(shù)據(jù)庫比對，60848個(gè)基因獲得注釋信息，其中38182條Unigene在GO數(shù)據(jù)庫中獲得注釋，24101條Unigene在KEGG數(shù)據(jù)庫中獲得注釋。GO和KEGG富集分析結(jié)果表明，差異表達(dá)基因主要調(diào)節(jié)核糖體代謝、植物激素信號轉(zhuǎn)導(dǎo)等生物學(xué)功能。

關(guān)鍵詞：旱柳;轉(zhuǎn)錄組分析;鹽脅迫;差異表達(dá)基因

中圖分類號：S718.43

文獻(xiàn)標(biāo)識碼：A

文章編號：1000-4440（2019）02-0271-11

土壤鹽漬化是阻礙植物生長并嚴(yán)重降低植物生產(chǎn)力的主要環(huán)境脅迫因素之一。根據(jù)糧農(nóng)組織和教科文組織的不完全統(tǒng)計(jì)資料，全世界約有4.0x108hm鹽堿地，而中國約有3.6x107hm。由于自然環(huán)境變化和人為因素，導(dǎo)致土壤退化，土壤鹽堿化面積逐漸增加。鹽脅迫對光合作用、能量生產(chǎn)和營養(yǎng)獲取等植物基本過程具有很大的負(fù)面影響，它破環(huán)細(xì)胞膜的完整性和各種酶的活性，這可能進(jìn)一步導(dǎo)致許多有害作用，例如缺水、滲透脅迫、二次氧化脅迫以及動(dòng)態(tài)平衡和離子毒性的破壞。許多研究者借助分子和基因組學(xué)的方法闡明耐鹽植物的耐鹽機(jī)制。通過轉(zhuǎn)錄組分析，許多耐鹽基因已經(jīng)被發(fā)現(xiàn)，涉及鹽脅迫信號轉(zhuǎn)導(dǎo)的幾種重要途徑，例如鹽敏感性蛋白（SOS）、鈣依賴性蛋白激酶（CDPK）和有絲分裂原激活蛋白（MAP）。

過去幾十年，不同植物耐鹽脅迫的轉(zhuǎn)錄組分析已被報(bào)道，比如棉花、大豆、冰葉日中花、草地早熟禾、歐洲油菜、沙漠楊樹和小麥。這些研究結(jié)果表明，與氧化還原過程和脯氨酸代謝有關(guān)的差異表達(dá)基因以及許多轉(zhuǎn)錄因子都在植物耐鹽中發(fā)揮重要作用。目前，一些模式植物，比如水稻、玉米和擬南芥等已經(jīng)進(jìn)行了轉(zhuǎn)錄組測序。這些模式植物的DNA序列信息和轉(zhuǎn)錄組序列信息都已經(jīng)很完善，可用于測序、繪圖和組裝基因等研究。然而，非模式木本鹽生植物的轉(zhuǎn)錄組研究很少。非模式木本鹽生植物基因組和轉(zhuǎn)錄組數(shù)據(jù)庫的缺乏很大程度上限制了其相關(guān)的研究。

柳樹（Salix L.）屬于楊柳科物種，其種類繁多，在中國就存在256種以上，其生長潛力大，頂端優(yōu)勢比較明顯，腋芽萌發(fā)力比較強(qiáng)，具有抗性強(qiáng)、適應(yīng)性廣等特性。柳樹根部是應(yīng)對鹽脅迫的首要部位，根部對Na*的選擇性吸收外排和區(qū)域化對柳樹響應(yīng)鹽脅迫具有重要作用。本研究利用llumina HiSeq高通量測序平臺對旱柳沿江柳和耐鹽旱柳9901的雜交F，代根部進(jìn)行轉(zhuǎn)錄組測序分析，比較耐鹽柳樹和鹽敏感柳樹中基因表達(dá)的差異，以期通過RNA-Seq技術(shù)獲得與耐鹽性相關(guān)的基因序列，為揭示柳樹耐鹽分子機(jī)制和培育柳樹耐鹽新品系提供參考。

1 材料與方法

1.1 材料及處理

試驗(yàn)材料為鹽敏感旱柳沿江柳（采自江蘇沿江地區(qū)農(nóng)業(yè)科學(xué)研究所）和耐鹽旱柳9901（采自山東省林業(yè)科學(xué)研究院）的雜交F代。2017年6月，將F，代進(jìn)行耐鹽分離群體試驗(yàn)，選擇15個(gè)耐鹽材料和15個(gè)鹽敏感材料，每5個(gè)材料混合為1個(gè)樣品，作為1個(gè)重復(fù)。經(jīng)0mmol/L和150mmol/LNaCl溶液處理4h后，取柳樹的新鮮根，立即放入液氮中保存?zhèn)溆谩?/p>

1.2 總RNA的提取及檢測

使用TIANCEN公司的多糖多酚植物試劑盒提取柳樹根的總RNA，樣品的純度、濃度和完整度分別用Nanodrop、Qubit2.0、Aglient2100分析。柳樹總RNA的提取、質(zhì)量檢測及無參考基因組序列的轉(zhuǎn)錄組高通量測序分析均由北京百邁客生物科技有限公司進(jìn)行。

1.3 cDNA文庫構(gòu)建和轉(zhuǎn)錄組測序

用檢測合格的柳樹總RNA樣品構(gòu)建eDNA文庫。以mRNA為模板，利用AMPure XP beads試劑盒依次進(jìn)行純化、末端修復(fù)和片段大小選擇，通過PCR富集得到柳樹eDNA文庫。eDNA文庫質(zhì)量用Agilent2100檢測，再通過llumina Hi-Seq平臺完成RNA-Seq測序分析。

1.4 De novo組裝及質(zhì)量控制

獲得的原始數(shù)據(jù)（Rawreads）經(jīng)過濾得到Cleanreads。將其用Trinity拼接，得到contigs。利用DeBruijin圖方法和測序Read信息，在片段集合中讀取轉(zhuǎn)錄本序列，最長contig即為最終轉(zhuǎn)錄本的Unigene。每個(gè)樣品測序產(chǎn)出不少于6Gb Clean data，Q30堿基百分比達(dá)到85%。

1.5 Unigene功能注釋

使用BLAST軟件將差異表達(dá)基因（DEG）序列與NR、Swiss-Prot、GOCOG、KOG、eggNOG4.5、KEGG數(shù)據(jù)庫分別比對，獲得相關(guān)注釋信息。

2 結(jié)果與分析

2.1 旱柳轉(zhuǎn)錄組測序與de novo組裝

利用lluminaHi-Seq對12個(gè)旱柳根樣品進(jìn)行轉(zhuǎn)錄組測序，共獲得83.81Gb Clean data，各樣品Clean data均達(dá)到6.31Gb以上（表1）。Q30堿基百分比均不小于89.21%，GC含量百分比在44%以上，可以用于下一步的de novo組裝。利用Trinity軟件將所有短序列拼接組裝得到107950條Unigene，總長度116257763bp，平均長度為1076.96bp。其中長度在500～1000bp的Unigenes所占比例最高，為29.39%，Unigene的N50為1690bp，組裝完整性較高（圖1），獲得的Mapped reads可用于后續(xù)分析。

2.2 旱柳差異基因表達(dá)量分析

在差異表達(dá)分析過程中采用了公認(rèn)有效的Benjamini-Hochberg方法對原有假設(shè)檢驗(yàn)得到的顯著性p值進(jìn)行校正，并最終采用校正后的p值，即FDR（False discovery rate）作為差異表達(dá)基因篩選的關(guān)鍵指標(biāo)，以降低對大量基因的表達(dá)值進(jìn)行獨(dú)立統(tǒng)計(jì)假設(shè)檢驗(yàn)帶來的假陽性。

在篩選過程中，將FDR小于0.01且差異倍數(shù)FC[兩樣品（組）間表達(dá)量的比值]≥2作為篩選標(biāo)準(zhǔn)。在鹽敏感柳樹中，以0mmol/LNaCl溶液處理為對照，在150mmol/L鹽脅迫下，共有1235個(gè)DEG表達(dá)量發(fā)生了改變，其中1130個(gè)DEG上調(diào)，105個(gè)DEG下調(diào);在耐鹽柳樹中，2629個(gè)DEG表達(dá)量發(fā)生了改變，其中2059個(gè)DEG上調(diào)，570個(gè)DEG下調(diào)。在150mmol/L鹽脅迫下，鹽敏感柳樹與耐鹽柳樹相對比，共有193個(gè)DEG表達(dá)量發(fā)生了改變，其中69個(gè)DEG上調(diào)，124個(gè)DEG下調(diào)（表2）。將不同組的差異基因進(jìn)行對比并繪制維恩圖，展示對比后各組特有的差異基因數(shù)量，以及共有的差異基因數(shù)量。其中鹽敏感柳樹在不同鹽濃度處理下與耐鹽柳樹在不同鹽濃度處理下一共有853個(gè)重復(fù)差異基因;鹽敏感柳樹和耐鹽柳樹在150mmol/L鹽脅迫下與鹽敏感柳樹在不同鹽處理下有59個(gè)差異基因，而與耐鹽柳樹在不同鹽處理下有99個(gè)差異基因，三者有5個(gè)共同的差異基因（圖2）。

2.3 旱柳差異表達(dá)基因功能注釋

通過選擇BLAST參數(shù)（E值≤x10～'）和HMMER參數(shù)（E值≤x10），將全部Unigene序列與COG、GO等8個(gè)數(shù)據(jù)庫進(jìn)行比對（表3、表4），最后得到60848個(gè)有注釋信息的基因。其中長度在1000bp以上的基因注釋效率為52.39%，而未得到注釋的Unigene為47102條（43.63%）。

2.4 旱柳差異表達(dá)基因GO分類和功能富集

將旱柳差異表達(dá)基因在GO數(shù)據(jù)庫中進(jìn)行比對，有38182條Unigene在GO數(shù)據(jù)庫中得到相關(guān)注釋，其中每2個(gè)樣品相比得到的差異基因在GO數(shù)據(jù)庫中注釋率最多達(dá)69.28%，注釋的基因參與生物學(xué)過程、細(xì)胞組分和分子功能調(diào)控（表5），功能組涉及到54個(gè)（圖3、圖4、圖5）。生物學(xué)過程的注釋顯示，Unigene表達(dá)差異顯著順序依次為代謝過程、細(xì)胞過程和單組織過程;在細(xì)胞成分分類信息中，細(xì)胞部分和細(xì)胞顯著富集;在分子功能分類中，富集程度最高的是生物合成和催化活性。

2.5 旱柳差異表達(dá)基因KEGG注釋和富集分析

有24101條（39.61%）Unigene在KEGG數(shù)據(jù)庫得到注釋，涉及209條KEGG代謝通路。其中鹽敏感柳樹組，以0mmol/LNaCl溶液處理為對照，在150mmol/L鹽脅迫下，獲得注釋的189條差異表達(dá)基因（DEG）可以被分配到77個(gè)KEGG通路（圖6）;150mmol/L鹽脅迫下，鹽敏感柳樹組與耐鹽柳樹組對比，獲得注釋的72條差異表達(dá)基因可以被分配到35個(gè)KEGG通路（圖7）。耐鹽柳樹組，以0mmol/LNaCl溶液處理為對照，在150mmol/L鹽脅迫下，獲得注釋的487條差異表達(dá)基因可以被分配到97個(gè)KEGG通路（圖8），其中在同一個(gè)KEGG通路中10個(gè)以上差異表達(dá)基因主要涉及細(xì)胞途徑、遺傳信號響應(yīng)、環(huán)境信號響應(yīng)、新陳代謝和生物系統(tǒng)5個(gè)方面。這5個(gè)方面中新陳代謝類最大，包括碳水化合物、氨基酸、多糖生物合成、能量、次生物質(zhì)等;生物系統(tǒng)類最小。同時(shí)，KEGG注釋結(jié)果也顯示，差異基因主要參與核糖體、植物激素信號轉(zhuǎn)導(dǎo)、淀粉與蔗糖代謝、碳代謝以及植物病原菌的相互作用等通路（表6）。

2.6 旱柳差異表達(dá)基因注釋結(jié)果分析

在鹽脅迫下，差異表達(dá)基因，上調(diào)2倍以上才有意義，現(xiàn)篩選高水平差異表達(dá)（log2FC達(dá)7倍以上）的基因（表7）。高水平差異表達(dá)的注釋基因中，在150mmol/L鹽脅迫下，鹽敏感柳樹相對于耐鹽柳樹類似葉綠體的基因BMK_Unigene_064460表達(dá)下調(diào)水平顯著。與對照相比，鹽脅迫下鹽敏感柳樹中高水平差異基因表達(dá)都是上調(diào)的，主要包括穩(wěn)定蛋白、細(xì)胞色素p450家族蛋白、細(xì)胞周期蛋白等。與對照相比，鹽脅迫下耐鹽柳樹中高水平差異基因表達(dá)都是上調(diào)的，主要包括水通道蛋白、TPR結(jié)構(gòu)蛋白、磷酸烯醇式丙酮酸羧激酶、谷胱甘肽s轉(zhuǎn)移酶等。

3 討論

當(dāng)植物受到干旱、高鹽、高溫、低溫等非生物脅迫時(shí)，植物細(xì)胞表達(dá)相關(guān)抗逆基因來適應(yīng)外界不利刺激。這些抗逆基因與各種代謝途徑、抗氧化過程等密切相關(guān)。本試驗(yàn)中，將鹽脅迫下的F，代和正常生長的F，代旱柳植株轉(zhuǎn)錄組進(jìn)行比較，發(fā)現(xiàn)兩者之間的基因表達(dá)具有較大差異，包括差異表達(dá)基因數(shù)量、表達(dá)量、發(fā)生上下調(diào)的基因數(shù)量。植株對鹽脅迫的抵御能力往往需要這些差異表達(dá)的基因來調(diào)節(jié)。我們對正常生長和鹽脅迫下的F，代旱柳植株.的根進(jìn)行了轉(zhuǎn)錄組測序，并以FDR小于0.01且差異倍數(shù)大于等于2的標(biāo)準(zhǔn)進(jìn)行篩選，分別在0-鹽敏（清水處理的鹽敏感柳樹）F代柳樹相對于150-鹽敏（150mmol/L處理的鹽敏感柳樹）F，代柳樹，150-鹽敏F，代柳樹相對于150-耐鹽（150mmol/L處理的耐鹽柳樹））F，代柳樹以及0-耐鹽（清水處理的耐鹽柳樹）F，代柳樹相對于150-耐鹽F，代柳樹中得到符合該標(biāo)準(zhǔn)的差異表達(dá)基因。0-鹽敏F，代柳樹相對于150-鹽敏F，代柳樹和0-耐鹽F，代柳樹相對于150-耐鹽F代柳樹的差異表達(dá)基因數(shù)量分別為1235和2629，說明耐鹽F代柳樹的耐鹽機(jī)制中參與調(diào)控的基因較多，可能有更好的抵御能力。150-鹽敏F代柳樹相對于150-耐鹽F代柳樹的差異表達(dá)基因數(shù)目為193，進(jìn)一步說明兩種F，代柳樹對鹽脅迫的抵御能力是有差異的。

不司組織器官、發(fā)育階段以及不同處理下樣本間基因表達(dá)的比較能夠有效地揭示分子機(jī)制。本研究中，我們對正常生長和鹽脅迫下的柳樹F，植株進(jìn)行轉(zhuǎn)錄組測序分析，結(jié)果顯示60848個(gè)Unigene通過COG、GO等8個(gè)數(shù)據(jù)庫比對得到注釋信息，其中有38182個(gè)Unigene在GO數(shù)據(jù)庫中得到相關(guān)注釋，24101個(gè)Unigene在KEGG數(shù)據(jù)庫中得到注釋。通過進(jìn)一步GO富集分析發(fā)現(xiàn)，大量差異表達(dá)基因在代謝過程、細(xì)胞過程、單組織過程、細(xì)胞部分生物合成以及催化活性等GO條目中顯著富集。通過KEGG顯著性富集分析發(fā)現(xiàn)差異表達(dá)基因主要在核糖體、植物激素信號傳導(dǎo)、淀粉與蔗糖代謝、碳代謝以及植物病原菌的相互作用等通路較為活躍，說明它們在柳樹響應(yīng)鹽脅迫中發(fā)揮重要作用。

另外，我們對log2FC達(dá)到7倍以上的差異表達(dá)基因進(jìn)行了注釋分析。結(jié)果表明，鹽敏感柳樹中高水平差異基因都是上調(diào)表達(dá)，主要包括蛋白質(zhì)、細(xì)胞色素p450家族蛋白以及細(xì)胞周期蛋白等。耐鹽柳樹中高水平差異基因也都是.上調(diào)表達(dá)的，主要包括水通道蛋白、TPR蛋白、磷酸烯醇式丙酮酸羧激酶以及谷胱甘肽s轉(zhuǎn)移酶等。因此，這些差異表達(dá)基因可能是柳樹耐鹽的關(guān)鍵因子。

總之，通過轉(zhuǎn)錄組分析，我們比較全面地了解了柳樹響應(yīng)鹽脅迫的轉(zhuǎn)錄水平變化，為進(jìn)一步研究柳樹耐鹽分子機(jī)制，培育柳樹耐鹽品系提供了重要的信息。

參考文獻(xiàn)：

[1] SHABALA S，CUIN T A. Potassium transport and plant salt tolerance ：[ J]. Physiol Plant，2008，133（4）：651-669.

[2] NAN G，ZHANG Y，LI S，et al. NaCl stress-induced transeriptomics analysis of Salix linearistipularis （Salix mongolica）[J]. J Biol Res （Thessalon），2016，23：1.

[3] CHINNUSAMY V，ZHU J，ZHU J K. Salt stress signaling and mechanisms of plant salt tolerance[ J]. Genet Eng，2006，27：141-177.

[4] MUNNS R，TESTER M. Mechanisms of salinity tolerance[J]. An-nu Rev Plant Biol，2008，59：651-681.

[5] CHEN S，POLLE A. Salinity tolerance of Populus[J]. Plant Biol（Stuttg），2010，12（2）：317-333.

[6] NISHIYAMA R，LE D T，WATANABE Y，et al. Transcriptome analyses of a salt-tolerant cytokinin-deficient mutant reveal differ-ential regulation of salt stress response by cytokinin deficiency[J].PLoS ONE，2012，7（2）：e32124.

[7] WARD J M，HIRSCHI K D，SZE H. Plants pass the salt[J].Trends Plant Sci，2003，8（5）：200-201.

[8] NISHIYAMA R，WATANABE Y，F(xiàn)UJITA Y，et al. Analysis of cytokinin mutants and regulation of cytokinin metabolic genes re-veals important regulatory roles of cytokinins in drought，salt andabscisic acid responses ，and abscisic acid biosynthesis[J]. PlantCell，2011 ，23（6）：2169-2183.

[9] KRASENSKY J，JONAK C. Drought，salt，and temperature stress-induced metabolic rearrangements and regulatory networks[ J ]. JExp Bot，2012，63（4）：1593- 1608.

[ 10] SANCHEZ-BARRENA MJ，MARTINEZ-RIPOLL M，ZHU J K，et al. The structure of the Arabidopsis thaliana SOS3 ：molecularmechanism of sensing calcium for salt stress response[ J]. J MolBiol，2005，345（5）：1253- 1264.

[ 11] NAKAGAMI H，PITZSCHKE A，HIRT H. Emerging MAP kinase pathways in plant stress sigalling[ J]. Trends Plant Sci，2005，10（7）：339-346.

[12] XUP，LIU z，F(xiàn)AN X，et al. De novo transcriptome sequencing and comparative analysis of differentially expressed genes in Gos-sypium aridum under salt stress[J]. Gene，2013，525（1）：26-34.

[13] BUYUK I，INAL B，ILHAN E，et al. Genome-wide identification of salinity responsive HSP70s in common bean[J]. Mol Biol Rep，2016，43（11）：1251-1266.

[ 14] TSUKAGOSHI H，SUZUKI T，NISHIKAWA K，et al. RNA-seq analysis of the response of the halophyte，Mesembryanthemum crys-tallinum （ice plant）to high salinity[J]. PLoS ONE，2015，10（2）：e0118339.

[ 15] BUSHMAN B S，AMUNDSEN K L，WARNKES E，et al. Tran-scriptome profiling of Kentucky bluegrass （Poa pratensis L.）ac-cessions in response to salt stress[ J]. BMC Genomics，2016，17：48.

[16] LONG W，ZOU X，ZHANG X. Transcriptome analysis of canola （Brassica napus）under salt stress at the germination stage [ J ].PLoS ONE，2015，10（2）：e0116217.

[17] ZHANGJ，JIANG D，LIU B，et al. Transcriptome dynamics of a desert poplar （Populus pruinosa）in response to continuous salinitystress[J]. Plant Cell Rep，2014 ，33（9）：1565-1579.

[18] GOYALE，AMITSK，SINGH RS，et al. Transcriptome profiling of the salt-stress response in Triticum aestivum cv. Kharchia Local[J]. Sci Rep，2016，6：27752.

[19] KUMARIS，SABHARWAL V P，KUSHWAHA H R，et al. Tran-scriptome map for seedling stage specific salinity stress response in-dicates a specific set of genes as candidate for saline tolerance inOryza sativa L.[J]. Funct Integr Genomics，2009，9（1）：109-123.

[ 20] OHTSU K，SMITH M B，EMRICHSJ，et al. Global gene expression analysis of the shoot apical meristem of maize （Zea mays L. ）[J]. Plant J，2007，52（3）：391-404.

[21] WEBER A P，WEBER K L，CARR K，et al. Sampling the Arabi-dopsis transcriptome with massively parallel pyrosequencing[J]. Plant Physiol，2007，144（1）：32-42.

[22]李敏，王瑩，李玉娟，等.2個(gè)柳樹親本耐鹽性比較試驗(yàn)[J].浙江農(nóng)業(yè)科學(xué)，2017（7）：1220- 1222.

[23] WICKER T，SCHLAGENHAUF E，GRANER A，et al. 454 sequencing put to the test using the complex genome of barley[J].BMC Genomics，2006，7：275.

[ 24] MOROZOVA O，HIRST M，MARRA M A. Applications of new sequencing technologies for transcriptome analysis[J]. Annu RevGenomics Hum Genet，2009，10：135-151.