孫源長，羅 琳，林 輝，林冬梅，林占熺

(福建農(nóng)林大學國家菌草工程技術研究中心，福州 350002)

【研究意義】隨著全球人口的急劇增長以及科學技術的快速發(fā)展，人們對能源的需求量逐漸增大。石油和煤炭等不可再生能源日益減少，因此迫切需要開發(fā)利用可再生能源，而生物質(zhì)能是其重要組成部分。目前,一些作物和生物量較大的植物作為生物能源被廣泛應用，但這兩者競爭生長，不適合一起種植在耕地中，將適應能力強、生物量大的植物種植到鹽漬地等生長環(huán)境惡劣的區(qū)域是很好的選擇。最近，研究表明全球超過10×108hm2的土地受鹽分影響，約占地球陸地表面的7%，并且還在持續(xù)蔓延。近期研究結果表明，伊朗西北部水土流失和鹽漬化問題嚴重[1]，而亞洲(中國、印度、巴基斯坦)、美洲(美國)、非洲(南非)等地的土壤鹽含量也嚴重超標[2]。土壤鹽分是影響植物生長和發(fā)育最嚴重的非生物脅迫之一，鹽度主要以兩種方式影響植物的生長發(fā)育。滲透脅迫是在高濃度鹽環(huán)境下根表面的水勢降低，導致植物吸水量減少，擾亂植物的正常生命過程；離子脅迫是植物細胞內(nèi) Na+/ K+的失衡引起的，許多生化反應將不能進行[3-4]。蘆竹又稱“荻蘆竹”，為禾本科蘆竹亞科蘆竹屬的多年生草本植物，已在亞洲、南歐和北非等地種植上千年，能夠適應各種不同環(huán)境，在溫帶地區(qū)可以長期在寒冷、鹽堿、干旱和澇漬等條件下生存而不會造成重大損失[5-7]。這種植物能夠耐受鹽脅迫，在高鹽度下仍保持較高的光合作用和生物量，不會對光化學反應以及葉片結構和完整性造成較大損傷[8]。因此，探究蘆竹莖的鹽響應方式，鑒定出參與鹽響應的重要基因，對未來分子育種，提高植物耐鹽性，減輕能源壓力具有重大意義?！厩叭搜芯窟M展】動物可以移動遷徙來應對即將和已經(jīng)發(fā)生的脅迫，而植物卻只能通過調(diào)節(jié)自身狀態(tài)來應對環(huán)境的變化[9]。自身狀態(tài)的調(diào)節(jié)主要是通過基因的表達調(diào)控來實現(xiàn)，并且這種表達調(diào)控是動態(tài)的。為應對各種生物和非生物脅迫，植物在長期進化過程中產(chǎn)生了多種自我保護和防御機制。例如當環(huán)境發(fā)生劇烈變化時，植物會調(diào)整一系列信號通路，來提高自身生存率[10-11]；通過關閉氣孔來調(diào)節(jié)葉片的光合過程[10]；通過調(diào)節(jié)體內(nèi)抗氧化系統(tǒng)消除多余的活性氧，維持體內(nèi)氧化還原平衡[10-11]；通過細胞滲透調(diào)節(jié)物質(zhì)改變細胞內(nèi)水勢，維持一定的滲透壓，從而保護植物在鹽脅迫下的持續(xù)生長[10]；植物激素通過協(xié)同或拮抗作用來調(diào)節(jié)自身[12]；轉(zhuǎn)錄因子(TF)與啟動子區(qū)域中順式作用元件發(fā)生特異性結合，進而調(diào)控靶基因的表達[13]。同樣，蘆竹主要也是通過這些方式來減輕外界不利環(huán)境對自身的影響?！颈狙芯壳腥朦c】目前關于蘆竹鹽響應的相關研究并不多，所以迫切需要對其進行深入研究以豐富蘆竹的耐鹽基因資源，為未來的分子育種奠定基礎。RNA測序 (RNA-Seq) 一直是研究應激耐受分子機制的有效方法，被廣泛應用于各個領域[14-18]。【擬解決的關鍵問題】本研究以不同濃度NaCl鹽溶液處理的“綠洲1號”蘆竹的莖為材料，利用RNA-seq手段分析，探究蘆竹莖的鹽響應方式，鑒定出參與鹽響應的重要基因，進一步加深對蘆竹耐鹽分子機制的理解。

1 材料與方法

1.1 供試材料

從福建農(nóng)林大學國家菌草工程技術研究中心選取生長6個月并且長勢一致的“綠洲1號”蘆竹(原產(chǎn)地：萊索托)組培幼苗，之后移栽到溫室(在30 ℃/25 ℃以14 h光照/10 h黑暗的光周期生長，濕度為60%，光照強度為350 μmol/(m2·s)并用1/5霍格蘭營養(yǎng)液水培1個月進行緩苗，之后分別用0、50、100、150和200 mmol/L NaCl鹽濃度溶液處理，處理時間為36 h。

1.2 RNA提取與轉(zhuǎn)錄組測序

分別收集不同NaCl鹽濃度溶液處理后的“綠洲1號”莖，每個重復由2～4株蘆竹自下往上數(shù)第1～8節(jié)混合獲取，共3次重復。樣品用液氮研磨成粉末狀混勻。

稱取樣品粉末0.2 g，用Trizol試劑提取莖的總RNA。取2 μL RNA在1%瓊脂糖凝膠上初步評估RNA降解程度及污染情況。通過Nano Drop 2000?分光光度計測定RNA的純度與濃度。使用Agilent2100(Agilent Technologies)檢測RNA完整性。對合格的樣品委托諾禾致源公司進行轉(zhuǎn)錄組測序。簡要流程如下：用帶有Oligo(dT)的磁珠富集mRNA，隨后加入fragmentation buffer將mRNA打斷成短片段，以mRNA為模板，合成雙鏈cDNA。之后進行末端修復，加A尾并連接測序接頭。最后進行PCR擴增，并用AMPure XP beads純化PCR產(chǎn)物，得到最終文庫。對合格的樣本在Illumina Novaseq 6000平臺上進行雙端測序，讀長為150 bp。

表1 用于qPCR反應的引物

1.3 轉(zhuǎn)錄組數(shù)據(jù)的質(zhì)量控制與組裝

為了得到更加完整的轉(zhuǎn)錄本，此次組裝額外使用了本實驗室同一品種的15個葉片轉(zhuǎn)錄組數(shù)據(jù)。原始數(shù)據(jù)(Raw reads)由FastQC(https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)軟件進行質(zhì)量檢測。任何低質(zhì)量的原始數(shù)據(jù)和接頭序列均通過Trimmomatic[19]軟件過濾，最終得到高質(zhì)量數(shù)據(jù)(Clean data)。使用Trinity-v2.13.2[20]軟件對轉(zhuǎn)錄組數(shù)據(jù)進行從頭組裝，利用CD-HIT[21]軟件去冗余(相似性設置為0.95)得到最終單基因簇(Unigenes)序列。

1.4 基因表達定量和差異表達分析

使用Bowtie2和RSEM對基因進行表達定量[22-23]。在R中使用DESeq2[24]進行差異表達基因(Differentially expressed genes，DEGs)的檢測。DEGs的過濾標準是|log2FoldChange|≥1并且P≤ 0.05。

1.5 GO和KEGG富集分析

使用clusterProfiler[25]R包對DEGs進行GO與KEGG富集分析。參數(shù)設置為校正后的P<0.05, Unigene的注釋來源為eggNOG5.0[26]數(shù)據(jù)庫。

1.6 轉(zhuǎn)錄因子鑒定

所有差異表達基因均通過iTAK(v1.6)[27](http://itak.feilab.net/cgi-bin/itak/online_itak.cgi)在線程序進行轉(zhuǎn)錄因子的預測，比對類型為核酸序列。

1.7 基因的qRT-PCR驗證

隨機選取5個差異表達基因驗證轉(zhuǎn)錄組數(shù)據(jù)是否可靠。根據(jù)Michele等[28]的實驗結果，選取較為穩(wěn)定的RPN6基因作為內(nèi)參基因，所有基因通過Primer 5.0軟件設計引物，qPCR引物見表1。使用TransScript?Uni All-in-One First-Strand cDNA Synthesis SuperMix for qPCR(One-Step gDNA Removal)試劑盒進行反轉(zhuǎn)錄。qPCR反應根據(jù)PerfectStart?Green qPCR SuperMIX試劑盒說明書進行。采用2-△△Ct法計算基因的相對表達量[29]，用GraphPad Prism 8 軟件進行圖像繪制。

2 結果與分析

2.1 測序數(shù)據(jù)的組裝結果

使用Trinity[20]軟件對Clean data進行De novo組裝，共產(chǎn)生971 309個轉(zhuǎn)錄本，GC含量為47.14%。對組裝后的結果文件用CD-HIT[21]軟件去冗余，得到最終的Unigenes。Transcript contigs的N30達到2053 bp, N50達到1228 bp，contig 平均長度為823.47 bp；Unigene contigs的N30達到1915 bp, N50達到1126 bp，contig平均長度為769.52 bp。轉(zhuǎn)錄本和基因的長度數(shù)量分布統(tǒng)計情況見圖1。以上結果表明此次轉(zhuǎn)錄本組裝效果優(yōu)良，可以進行后續(xù)分析。

圖1 轉(zhuǎn)錄本和單基因簇長度數(shù)量分布統(tǒng)計Fig.1 Length and number distribution statistics of transcripts and unigenes

2.2 不同鹽濃度脅迫下的基因表達量化和差異表達分析

在R中通過 DESeq2包分析“綠洲1號”響應鹽脅迫的差異表達基因。在葉片中檢測到的差異表達基因為7305個，其數(shù)量分布情況見圖2。隨著鹽濃度的增加，上調(diào)差異表達基因、下調(diào)差異表達基因與總的差異表達基因的數(shù)量都逐漸增加。上調(diào)的差異表達基因數(shù)目始終大于下調(diào)的差異表達基因數(shù)目。以上結果表明“綠洲1號”莖在響應鹽脅迫時，隨著鹽濃度的增加其調(diào)控過程越來越復雜，并且主要通過基因的上調(diào)起作用。

圖2 不同鹽濃度下莖差異表達基因的數(shù)目Fig.2 Number of differentially expressed genes in stems under different salt concentrations

由圖3可知，在莖中有104個基因是不同濃度鹽處理下的差異表達基因所共有的，將它們作為蘆竹鹽響應核心差異表達基因。隨著鹽濃度的增加，各個組別中特有的差異表達基因數(shù)目也逐漸增加。stem_CK_vs_50 mmol/L到stem_CK_vs_100 mmol/L、stem_CK_vs_150 mmol/L、stem_CK_vs_200 mmol/L的特異性差異表達基因數(shù)目分別從111個上升到292、1141、3347個。這可能是由于隨著脅迫程度的加深，越來越多特定的信號傳導機制被激活。以上結果表明莖存在一些較為保守的基因來響應鹽脅迫，并且隨著鹽濃度的增加越來越多的特異性差異表達基因參與到表達調(diào)控中來。

圖3 不同鹽濃度下莖差異表達基因的Venn圖Fig.3 Venn plot of differentially expressed genes in stems under different salt concentrations

2.3 差異表達基因的GO和KEGG富集分析

由圖4可知，在蘆竹莖響應鹽脅迫過程中，50 mmol/L 鹽濃度下差異表達基因主要富集到鉀離子穩(wěn)態(tài)、硝酸鹽響應、細胞壁修飾、水楊酸合成相關條目。100 mmol/L 鹽濃度下，臭氧、泛素化、滲透壓調(diào)節(jié)相關條目富集到的基因數(shù)目最多。150 mmol/L鹽濃度下富集到對硝酸鹽、無機陰離子、光形態(tài)、端粒酶活性等條目。與150 mmol/L鹽濃度下富集結果類似，200 mmol/L鹽濃度下除了富集到無機陰離子轉(zhuǎn)運條目還富集到與水響應、光反應調(diào)節(jié)相關的條目。

圖4 莖中差異表達基因的GO富集Fig.4 GO enrichment of differentially expressed genes in stems

由圖5可知，在莖中，50 mmol/L鹽濃度只富集到了單萜生物合成途徑；而100、150和200 mmol/L鹽濃度富集到氮代謝、晝夜節(jié)律、萜類合成、玉米素合成和維生素B6代謝等通路。

圖5 莖中差異表達基因的KEGG富集Fig.5 KEGG enrichment of differentially expressed genes in stems

2.4 鹽響應核心差異表達基因分析

由圖6可知，在蘆竹莖中大部分鹽響應核心差異表達基因上調(diào)，只有少數(shù)基因下調(diào)；并且與下調(diào)基因相比上調(diào)基因的差異倍數(shù)更大。說明蘆竹響應鹽脅迫主要是通過基因上調(diào)起作用,并且上調(diào)基因的反應更加強烈。在莖中88315_c0_g1_i1、1200_c2_g1_i1、13330_c0_g1_i8、46571_c0_g1_i7、3434_c0_g1_i11和713_c0_g1_i32的差異倍數(shù)最大。這些基因在“綠洲1號”莖基礎鹽響應過程中發(fā)揮極其重要的作用，注釋的具體情況見表2。

圖6 莖中核心差異表達基因的差異倍數(shù)分析Fig.6 Fold analysis of core genes expression differences in stems based on log2FoldChange

表2 莖中重要的鹽響應核心差異表達基因

2.5 轉(zhuǎn)錄因子家族鑒定

對莖參與鹽響應表達調(diào)控的差異表達基因進行轉(zhuǎn)錄因子預測，一共預測到41個轉(zhuǎn)錄因子基因家族，家族成員數(shù)目最多的是AP2/ERF-ERF家族，共94個；其次是MYB-related家族和WRKY家族，家族成員分別為65、57個(圖7)。

圖7 莖中差異表達基因轉(zhuǎn)錄因子預測數(shù)目Fig.7 Predicted numbers of transcription factors of differentially expressed genes in stems

為探究AP2/ERF-ERF、MYB-related和WRKY轉(zhuǎn)錄因子家族的表達模式，對其基因表達差異倍數(shù)進行分析，AP2/ERF-ERF的log2FoldChange在-1.85～10.88；MYB-related的log2FoldChange在-3.00～7.85；WRKY的log2FoldChange在-1.58～6.13。并且這3類轉(zhuǎn)錄因子在響應鹽脅迫時大部分呈上調(diào)狀態(tài)，說明其主要通過增加基因的表達起作用(圖8)。

A：AP2/ERF-ERF；B：MYB-related；C：WRKY圖8 轉(zhuǎn)錄因子家族的差異倍數(shù)分析Fig.8 Fold analysis of transcription factor families in stems based on log2FoldChange

對前8個差異表達倍數(shù)較大的轉(zhuǎn)錄因子進行分析(表3)，AP2/ERF-ERF、MYB-related和WRKY轉(zhuǎn)錄因子家族差異倍數(shù)最大的轉(zhuǎn)錄因子分別為：18947_c0_g1_i3、7562_c1_g1_i2和24870_c0_g1_i2，平均log2FoldChange達到7.88、5.81和5.28,它們在蘆竹響應鹽脅迫時可能發(fā)揮著重要作用。

表3 差異表達倍數(shù)較大的轉(zhuǎn)錄因子

2.6 實時熒光定量PCR驗證

隨機挑選5個差異表達基因?qū)D(zhuǎn)錄組數(shù)據(jù)進行qRT-PCR驗證，可以看出這些基因的相對表達量和轉(zhuǎn)錄數(shù)據(jù)的基因表達模式基本一致(圖9)，也說明轉(zhuǎn)錄組數(shù)據(jù)較為可靠。

左側縱坐標為qPCR相對表達量，右側縱坐標為TPM值。A：13361_c0_g1_i15；B: 14873_c0_g1_i6；C: 70740_c0_g1_i2；D: 129260_c0_g2_i1；E：122017_c0_g1_i1The left ordinate is the relative expression level of qPCR, and the right ordinate is the TPM value. A: 13361_c0_g1_i15; B: 14873_c0_g1_i6; C: 70740_c0_g1_i2; D: 129260_c0_g2_i1; E: 122017_c0_g1_i1圖9 差異表達基因qRT-PCR 驗證Fig.9 qRT-PCR analysis of differentially expressed genes

3 討 論

如果為植物提供如光照、溫度和水分等適宜的生長環(huán)境，植物就會相對高產(chǎn)，反之有可能導致植物發(fā)育不良甚至死亡。而產(chǎn)量與植物發(fā)育過程中的非生物脅迫條件呈負相關[30]。在全球范圍內(nèi)，由于非生物脅迫(干旱、鹽漬、洪水、霜凍、營養(yǎng)缺乏等因素)導致的損失估計為每年1.2×1011美元[31](http://www.fao.org/docrep/008/y5800e/Y5800E06.htm)。鹽脅迫是降低作物產(chǎn)量最重要的非生物脅迫之一，因此，對植物進行鹽脅迫的相關研究將有利于篩選出耐鹽的關鍵基因，對提高植物的耐鹽性和產(chǎn)量具有重要意義。蘆竹具有較強的抗逆性，能夠在高鹽環(huán)境下生長。本實驗以“綠洲1號”蘆竹的莖為材料，使用不同濃度的NaCl鹽溶液進行處理，利用RNA-seq手段分析蘆竹莖的鹽響應方式，鑒定出重要的鹽響應相關基因，為以后的遺傳改良提供指導。

3.1 蘆竹鹽響應不同鹽環(huán)境的轉(zhuǎn)錄差異

本研究共檢測到7305個鹽脅迫差異表達基因，隨著鹽溶液濃度的增加,差異表達基因的數(shù)目也相繼增加，并且越來越多的特異性差異表達基因參與到表達調(diào)控中來。以上結果表明不同濃度的鹽環(huán)境均可誘導大量的基因來響應鹽脅迫，隨著濃度的增加表達調(diào)控的復雜程度也可能隨之上升。

對差異表達基因進行GO與KEGG富集分析，與代謝產(chǎn)物相關的條目主要富集到與萜類合成、玉米素合成、亞油酸代謝和水楊酸生物合成途徑，這些代謝物質(zhì)作為調(diào)節(jié)劑在植物應激過程中發(fā)揮著重要作用[32-34]，對提高蘆竹的鹽耐受性有重要作用。無論從低鹽還是高鹽環(huán)境都富集到與萜類合成相關的條目，可以推測蘆竹能通過快速調(diào)節(jié)萜類化合物合成來應對鹽環(huán)境的變化。富集到鉀離子穩(wěn)態(tài)、與細胞壁相關[35]、調(diào)節(jié)對滲透壓的反應、對鹽脅迫反應的正調(diào)節(jié)、應對缺水的正調(diào)控、無機陰離子跨膜轉(zhuǎn)運和磷酸根離子跨膜轉(zhuǎn)運等條目，將對蘆竹感應鹽脅迫、維持細胞的滲透壓穩(wěn)態(tài)和調(diào)節(jié)相關的信號通路有所幫助。在鹽濃度達到200 mmol/L時，即較高鹽環(huán)境下，還富集到了維生素B6代謝途徑。維生素B6近些年來被廣泛研究，它是多種代謝酶的必需輔酶，參與植物的生長發(fā)育、生物脅迫和非生物脅迫等一系列過程。它還是一種有效的抗氧化劑，研究表明維生素B6在植物體內(nèi)抗氧化防御中具有作用，在強光下限制1O2積累并防止1O2介導的氧化損傷[36]。

3.2 蘆竹鹽響應過程中的關鍵基因

共篩選到104個鹽響應的核心差異表達基因，其中大部分基因上調(diào)，只有少數(shù)基因下調(diào)，說明這些基因主要是通過增加表達量來響應鹽脅迫。對其中差異倍數(shù)最大的基因進行注釋，其中GATL1編碼糖基轉(zhuǎn)移酶，有助于木聚糖的生物合成。此外，植物糖基轉(zhuǎn)移酶能夠?qū)⑻翘砑拥叫》肿又?，例如植物激素[37]、水楊酸[38]和類黃酮[39]。而木聚糖與細胞壁合成相關，細胞壁參與了鹽分脅迫的感知和響應，對于細胞的壓力保護和反應都至關重要；CRF4是ERF/AP2 轉(zhuǎn)錄因子家族，之前的研究表明它在響應非生物脅迫中起著重要作用[40]；SF3B5A編碼剪接因子3B亞基5/RDS3復合亞基，在選擇性剪接中起作用；SYR1為系統(tǒng)素受體，主要位于質(zhì)膜上，有助于植物防御過程中系統(tǒng)素介導信號傳導的調(diào)節(jié)[41]。

3.3 蘆竹鹽響應過程中起重要作用的轉(zhuǎn)錄因子

植物編碼大量轉(zhuǎn)錄因子，它們主要按照 DNA 結合結構域分類[42]。越來越多的證據(jù)表明，轉(zhuǎn)錄因子可調(diào)節(jié)多種生物過程，例如各種生物及非生物脅迫、光和壓力信號、種子成熟和花發(fā)育等[43]。目前已經(jīng)鑒定到了一系列的轉(zhuǎn)錄因子家族在鹽脅迫過程中發(fā)揮著重要作用，例如bZIP、AP2/ERF-ERF、MYB、WRKY、NAC和bHLH等[44]。

為篩選出在蘆竹鹽響應過程中起重要作用的轉(zhuǎn)錄因子，本實驗對差異基因進行研究，共注釋到41個轉(zhuǎn)錄因子基因家族，推測它們可能在蘆竹對鹽脅迫的反應中起重要作用，其中成員數(shù)目最多的包括AP2/ERF-ERF、MYB-related和WRKY家族。AP2/ERF-ERF家族是植物所特有的一類基因家族，該家族中的每一個成員都具有單個或2個AP2結構域，主要參與植物的生長發(fā)育、細胞周期和生物及非生物脅迫過程，廣泛參與脫落酸、赤霉素和細胞分裂素等信號通路。MYB-related家族是植物最大的轉(zhuǎn)錄因子家族之一，根據(jù) DNA 結合結構域相鄰重復的數(shù)量可分為四大類，在植物的發(fā)育和脅迫反應中起重要作用[45]。WRKY家族是調(diào)節(jié)植物生物過程信號網(wǎng)絡的重要組成部分，通常作為抑制因子和激活因子發(fā)揮作用，在調(diào)節(jié)植物的許多脅迫反應中起關鍵作用。例如，在水稻中，熱休克誘導型HSP101啟動子驅(qū)動的OsWRKY11過表達將會導致耐熱和耐旱性增強[46]。同樣，OsWRKY45的過表達除了提高抗病性外，還導致耐鹽性和耐旱性增強[47]。此外，NAC、bZIP、bHLH等轉(zhuǎn)錄因子家族也在本研究中被發(fā)現(xiàn)，這與前人的研究一致，并且這些轉(zhuǎn)錄因子絕大部分通過上調(diào)表達量來響應鹽脅迫。以上結果表明AP2/ERF-ERF、MYB-related和WRKY轉(zhuǎn)錄因子家族在“綠洲1號”蘆竹響應鹽脅迫的過程中可能具有重要的調(diào)控功能。

4 結 論

對0、50、100、150和200 mmol/L NaCl鹽溶液進行處理的“綠洲1號”蘆竹的莖進行轉(zhuǎn)錄組測序與分析，共篩選到10個差異表達倍數(shù)較大的基因和一些差異表達倍數(shù)較大的轉(zhuǎn)錄因子，它們可能在蘆竹響應鹽脅迫的過程中發(fā)揮著重要作用。本研究豐富了蘆竹的耐鹽基因資源，為分子育種提供了理想的改造靶點。