綜合生理學和轉錄組學揭示蘆竹耐鹽的候選基因

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

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

隨著全球人口的快速增長以及人們生活水平的不斷提高，對能源的需求也顯著增加。但是化石燃料和不可再生能源是有限的，并且在日益減少，因此迫切需要大力發展和利用可再生能源，而生物質能是其重要的組成部分。目前馬鈴薯(Solanumtuberosum)、甘蔗(Saccharumofficinarum)、小麥(Triticumaestivum)等作物和一些生物量較大的草本植物作為生物燃料應用十分廣泛。雖然這些草本植物可以在耕地中正常生長，但是其與農作物競爭生長，不適宜一起種植，因此將適應能力強的草本植物種植到鹽漬地、干旱地區和土壤性質差的區域是個不錯的選擇。其中，土壤鹽漬化是一個嚴重的環境問題，目前世界上有100多個國家的土壤受到影響，面積超過了9億公頃[1]。它會顯著影響作物產量和區域的農業生產。蘆竹隸屬于禾本科蘆竹亞科蘆竹屬，具有高生物量、高蛋白、高抗逆性等特點，在合適的環境下，蘆竹第2年就可以生產出高達40 t·hm-2·a-1的干物質生物質能，對應的能量值約為150 MWh·hm-2·a-1[2]，可見蘆竹的生產潛力巨大，并且其在高鹽度條件下也可實現較為可觀的生物產量，因此將它們種植于鹽漬地中不僅可以改善生態環境而且它們不與糧食作物競爭，有利于克服糧食安全風險。之前Lloyd L. Nackley等人對蘆竹進行過鹽脅迫實驗，使它們在一系列鹽濃度(0～42 dS·m-1NaCl)環境中生長60天，結果顯示，在最高鹽度下，蘆竹整體生長減少了80%以上。然而，死亡率為零，表明蘆竹能夠耐受較高濃度的鹽環境[3]。雖然已知蘆竹具有較強的耐鹽性，但其分子機制和一些重要的鹽響應基因尚未可知，制約了對其耐鹽高產的遺傳改良。

近10年來，隨著高通量測序技術的飛速發展，其成本不斷下降，越來越多的研究人員對植物進行測序并通過生物信息學分析篩選到一些重要的基因與信號通路[4-7]，為相關的育種工作提供了更加有效的方案。本研究對蘆竹采用不同濃度梯度的鹽溶液處理，通過生理學與轉錄組學進行研究，篩選出參與鹽響應的重要基因，豐富蘆竹的耐鹽基因資源，未來可以針對這些基因構建突變體，對這些基因進行深入研究，為耐鹽品種的培育提供科學參考。

1 材料與方法

1.1 試驗材料

從福建農林大學國家菌草工程技術研究中心菌草圃選取了6個月大的長勢一致的蘆竹組培苗，移栽到溫室用1/5霍格蘭營養液培養1個月進行緩苗，使植物在14 h光照/10 h黑暗的光周期下生長，濕度為60%，光照強度為350 μmol·m-2·S-1。待適應環境后分別進行0 mmol·L-1，50 mmol·L-1，100 mmol·L-1，150 mmol·L-1和200 mmol·L-1NaCl鹽濃度溶液的處理，處理時間為36小時。

1.2 生理指標的測定

測定鹽脅迫后葉片的超氧化物歧化酶、過氧化物酶、過氧化氫酶活性，脯氨酸、可溶性糖、丙二醛含量和Fv/F0指標。每個重復由2～4株蘆竹從上至下第3～5片葉片混合得到，共3個重復。數據分析軟件采用GraphPad Prism 8。采用索萊寶公司超氧化物歧化酶(SOD)試劑盒測定SOD活性(型號：BC0175)，過氧化物酶(POD)試劑盒測定POD活性(型號：BC0095)，過氧化氫酶(CAT)試劑盒測定CAT活性(型號：BC0205)，脯氨酸(Pro)含量試劑盒測定Pro含量(型號：BC0295)，植物可溶性糖(Soluble sugar)含量試劑盒測定可溶性糖含量(型號：BC0035)，丙二醛(MDA)含量試劑盒測定MDA含量(型號：BC0025)。用英國Hansatech公司植物效率儀(型號：PocketPEA)測定2-4株蘆竹從上至下第3-5片葉片的中部光系統II的最大光化學效率(Fv/F0)，每片葉子測定4次，每個重復的測定結果取平均值，共3個重復。

1.3 RNA的提取與轉錄組數據的獲得

收集0，50，100，150，200 mmol·L-1NaCl鹽濃度溶液處理過的葉片，葉片中每個重復由2-4株蘆竹的從上至下數第2片葉片混合得到，共3個重復。

取各個樣品各0.2 g，用Trizol試劑盒提取總RNA。在1%瓊脂糖凝膠上評估RNA降解與污染。RNA的純度通過Nano Drop 2000(Thermo scientific)分光光度計檢查。RNA完整性通過Agilent 2100(Agilent Technologies)檢測。對合格的樣本在Illumina Novaseq 6000平臺進行雙端測序，讀長為150 bp。

1.4 轉錄組數據的質量控制及De novo組裝

為了得到更加完整的轉錄本，此次組裝額外使用了本實驗室相同品種的15個蘆竹莖的轉錄組數據。原始數據由FastQC進行質量檢測。任何低質量的原始數據和接頭序列都通過Trimmomatic[8]軟件進行過濾，最終得到高質量的數據(Clean data)。使用Trinity-v2.13.2[9]軟件對轉錄組數據合并進行De novo組裝，組裝得到的轉錄組文件利用CD-HIT[10]軟件去冗余(相似性設置為0.95)得到最終Unigene序列。

1.5 基因表達量化和差異表達分析

使用Bowtie2比對軟件和RSEM進行基因表達定量，以通過Perl腳本align_and_estimate_abundance.pl和來自Trinity軟件的-est_method RSEM獲得讀取計數的數量[11-12]。隨后，使用R包DESeq2[13]進行差異表達基因的分析。DEGs的過濾標準是|log2FoldChange|≥1并且p-value adjusted(padj)≤0.05。

1.6 GO和KEGG富集分析

為了研究DEGs的生物學意義，使用clusterProfiler[14]R包對其進行GO與KEGG富集分析。參數設置為校正后的P值小于0.05,基因的注釋來源為eggNOG5.0[15]數據庫。

1.7 加權基因共表達網絡分析

WGCNA(R包中的WGCNA v1.70-3)[16]用于探索基因與生理指標之間的關系，以及基因與基因之間的關系。參照前人方法，通過絕對中值差異篩選出具有4分位數(Q1)表達和TPM值大于2的基因輸入WGCNA進行網絡的構建[17]。最小模塊大小設置為30，最大模塊大小設置為20000，merge cut height設置為0.25，合并相似度為0.8的模塊，其他參數均為默認值。Cytoscape[18]工具用于構建和可視化交互網絡。Turquoise和green模塊的權重分別設置為0.18與0.24。

1.8 轉錄組數據驗證

通過qRT-PCR驗證轉錄組的準確性。本研究隨機選取5個DEGs進行驗證，以證明轉錄組數據的可靠性。根據Michele Poli等人的實驗結果，選取較為穩定的RPN6基因作為內參基因[19]，從Phytozome13(https://phytozome-next.jgi.doe.gov/)網站獲取高粱的RPN6基因的cds序列(S.bicolor v3.1.1|Sobic.006G092800.2 CDS)。之后從蘆竹最長轉錄本的cds中挑選與高粱RPN6基因相似度最高的作為蘆竹的RPN6 cds序列。所有的基因均采用Primer5.0軟件設計引物。根據制造商的說明，使用TransScript?Uni All-in-One First-Strand cDNA Synthesis SuperMix for qPCR(One-Step gDNA Removal)試劑盒進行反轉錄，使用前將各組分別點甩離心。在PCR管中建立20 μL體系，得到的產物直接用于qPCR反應。qPCR反應根據PerfectStart?Green qPCR SuperMIX試劑盒說明書進行。相對定量結果采用2-△△Ct法計算基因的相對表達量[20]。qPCR引物見表1。

表1 用于qPCR反應的引物Table 1 Primers for qPCR reactions

2 結果與分析

2.1 鹽脅迫對蘆竹生理生化的影響

隨著NaCl鹽濃度的改變，蘆竹的脯氨酸含量與可溶性糖含量發生了顯著變化。其中脯氨酸含量在鹽濃度為200 mmol·L-1處上升至最高值(P<0.01)。而可溶性糖含量在150 mmol·L-1處達到了最高值，在200 mmol·L-1處顯著下降到最低值(P<0.01)。具體情況如圖1所示。

圖1 不同鹽濃度下蘆竹的生理變化Fig.1 Physiological changes of Arundo donax under different salt concentrations

2.2 測序數據的組裝結果

使用Trinity-v2.13.2[9]軟件對所有Clean data進行從頭組裝，共產生971 309個轉錄本，GC含量為47.14%。之后對組裝好的轉錄本用CD-HIT[10]軟件去冗余，得到最終的Unigene。表2與表3是對組裝結果的統計，transcript contigs的N50達到了1 228 bp，unigene contigs的N50達到了1 126 bp，表明組裝效果較好，數據值得進一步分析。

表2 基于所有transcript contigs的統計數據Table 2 Stats based on all transcript contigs

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

在R中通過DESeq2包分析了蘆竹響應鹽脅迫的DEGs。DEGs的過濾標準是|log2FoldChange|≥1并且padj≤0.05。在葉片中檢測到的差異表達基因為8 745個，其分布如圖2所示?？梢钥闯鲭S著鹽濃度的增加，在葉片中無論是上調差異表達基因，下調差異表達基因還是總的差異表達基因數量都是逐漸增加的。不高于100 mmol·L-1鹽濃度時都是上調的差異表達基因數目小于下調的差異表達基因數目，當鹽濃度大于等于150 mmol·L-1時上調的差異表達基因數目大于下調的差異表達基因數目。這些結果表明蘆竹在響應鹽脅迫時，隨著鹽濃度的增加調控過程也越來越復雜。

表3 基于所有unigene contigs的統計數據Table 3 Stats based on all unigene contigs

圖3可以較為直觀地看出葉片各個組間差異表達基因的重疊數量。結果表明在葉片中有166個基因是不同濃度鹽處理下的差異表達基因所共有的，這些基因可能是鹽響應核心差異表達基因。并且隨著鹽濃度的增加，各個組別中特有的差異表達基因數目也逐漸增加。這可能是由于隨著脅迫程度的加深，越來越多特定的信號通路受到調控。例如在葉片中，leaf_0 mmol·L-1_vs_50 mmol·L-1到leaf_0 mmol·L-1_vs_100 mmol·L-1到leaf_0 mmol·L-1_vs_150 mmol·L-1再到leaf_0 mmol·L-1_vs_200 mmol·L-1的特有差異表達基因數目從295個上升到653個、1 226個、3 685個。從以上結果可以看出，葉片存在一些較為保守的差異表達基因來響應鹽脅迫，并且隨著鹽濃度的加深有越來越多的特異性差異表達基因參與到表達調控中來。

圖2 不同鹽濃度下葉片差異表達基因的數目Fig.2 Number of differentially expressed genes in leaves under different salt concentrations

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

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

對166個在葉片中的鹽響應核心差異表達基因進行差異倍數分析。在葉片中鹽響應的核心差異表達基因大部分是上調的，小部分是下調的；并且這些上調基因相對于下調基因的差異倍數更大(圖4)。說明響應鹽脅迫主要是通過基因的上調起作用的，這些上調基因相對于下調基因的反應更加強烈。

對葉片中差異倍數最大的前10個核心差異表達基因進行注釋，這些基因很可能在蘆竹基礎鹽響應過程中發揮著極其重要的作用。其中差異表達倍數最大的414_c1_g2_i5編碼去極化激活的Ca2+通道，與鈣離子轉運相關；0_c1_g3_i2編碼富含亮氨酸的重復蛋白激酶家族成員，與植物的防御相關。具體情況見表4。

圖4 葉片中核心差異表達基因的差異倍數分析Fig.4 Fold analysis of core genes expression differences in leaves based on log2FoldChange

表4 葉片中重要的鹽響應核心差異表達基因Table 4 Important salt-responsive core differentially expressed genes in leaves

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

對葉片在不同鹽濃度下的差異表達基因進行GO與KEGG富集分析。GO富集結果(圖5)表明：在葉片響應鹽脅迫中，50 mmol·L-1與100 mmol·L-1鹽濃度的差異表達基因主要富集到了與硝酸鹽、一氧化氮、活性氮和活性氧相關的條目。與之不同的是隨著鹽濃度的增加達到150 mmol·L-1，一些富集集中在了與光合作用、脯氨酸相關的條目。當鹽濃度上升到200 mmol·L-1時，除了富集到了與光合作用相關的條目還富集到了細胞分裂素、支鏈氨基酸相關途徑。KEGG富集結果(圖6)表明：在葉片中，50 mmol·L-1鹽濃度只富集到了氮代謝通路；而100 mmol·L-1，150 mmol·L-1和200 mmol·L-1鹽濃度還富集到了其他條目：例如類黃酮生物合成、芪類與二芳基庚烷類和姜酚生物合成、光合作用、光合作用-天線蛋白、玉米素生物合成和硒化合物代謝等代謝通路。

圖5 葉片中差異表達基因的GO富集結果Fig.5 GO enrichment results of differentially expressed genes in leaves

圖6 葉片中差異表達基因的KEGG富集結果Fig.6 KEGG enrichment results of differentially expressed genes in leaves

2.6 WGCNA分析

植物在響應鹽脅迫時會發生特定基因的表達調控，這些將反應在某些生理特征中。為了探究與特定性狀相關的基因集的關鍵(hub)基因,研究蘆竹特異性的鹽響應應答過程，采用WGCNA分析蘆竹葉片響應鹽脅迫的表達譜，檢測基因與生理指標以及模塊內基因之間的關系。

2.6.1模塊的篩選 結合生理指標和葉片的15個轉錄組數據，通過絕對中值差異篩選出具有4分位數(Q1)表達和TPM值大于2的共14 013個基因輸入WGCNA進行網絡的構建。通過計算基因間鄰接函數確定軟閾值為12。并且聚類得到38個模塊，其中turquoise模塊與Pro的相關性最高，相關性為0.87；green模塊與Pro的相關性也較高，相關性分別為0.79。具體見圖7。

2.6.2hub基因的鑒定 將2個目標模塊的hub差異表達基因導入Cytoscape構建可視化交互式網絡,結果如圖8所示。Turquoise模塊的240249_c0_g1_i1，104188_c0_g1_i3，219926_c0_g1_i1，18316_c0_g1_i4和7477_c11_g1_i1處于網絡的中心位置；Green模塊的10454_c0_g1_i12，6117_c0_g1_i16，19205_c0_g1_i4，23765_c0_g1_i1和86581_c0_g1_i9處于網絡的中心位置。對每個模塊Degree程度最強的基因進行注釋，結果見表5，這些被認為核心hub基因，其中與Pro相關的turquoise，green模塊的hub基因編碼核堿基抗壞血酸轉運蛋白和外膜孔蛋白等相關蛋白質，還包括一些具有某些特定結構的蛋白例如RING/U-box超家族蛋白和RRM結構域蛋白等等；而對2個模塊的核心hub基因進行表達量分析(各組TPM平均值取log2)，圖9結果顯示turquoise與green模塊的基因表達量隨著鹽濃度的加深呈上調趨勢。

圖7 模塊與生理指標關聯圖Fig.7 Correlation diagram between modules and physiological indicators

圖8 網絡分析2個模塊中的樞紐基因Fig.8 Net-work analysis the hub genes in two modules

表5 兩個模塊中前5個hub基因的注釋結果Table 5 Annotation results of the top 5 hub genes in the two modules

圖9 2個模塊中排名前5個hub基因的表達量分析Fig.9 Expression analysis of the top 5 hub genes in the two modules

2.7 qRT-PCR驗證轉錄組數據

由于RNA-seq具有一定的缺陷，可能會得到假陽性的結果，故需要利用qRT-PCR技術進行檢驗。對葉片的轉錄組數據隨機挑選5個差異基因進行qRT-PCR驗證，可以看出這些基因的相對表達量與轉錄數據的表達模式是基本一致的(圖10A～10E)，并且基因表達變化數據相關系數R2達到了0.82(圖10F)，也說明了轉錄組分析結果與本實驗中qRT-PCR所確定的表達模式具有較好的一致性。

3 討論

3.1 鹽脅迫對蘆竹的生理影響

植物在鹽脅迫的情況下會誘導抗氧化系統的激活，包括抗氧化化合物，例如類胡蘿卜素和抗壞血酸鹽、谷胱甘肽、α-生育酚和幾種參與活性氧(ROS)解毒的酶，例如超氧化物歧化酶、過氧化物酶、過氧化氫酶等[21]。而一些滲透調節物質也在抗逆過程中發揮著重要作用，如脯氨酸[22]與可溶性糖。脯氨酸是一種具有特殊構象的氨基酸，是初級代謝所必需的。壓力環境導致植物中脯氨酸的過量產生，進而通過維持細胞膨脹或滲透平衡來賦予壓力耐受性；穩定膜，從而防止電解質泄漏；使活性氧(ROS)濃度在正常范圍內，防止植物氧化爆發。MDA是脂質過氧化的細胞毒性產物，也是評價自由基產生和組織損傷的指標之一。Fv/F0是Fv與F0的比值，也是Fv/Fm(最大熒光)的一種表達方式，表示光系統II的最大光化學效率。

圖10 差異表達基因qRT-PCR驗證Fig.10 qRT-PCR analysis of differentially expressed genes注：圖A～E中左側縱坐標為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_i1；F,相關性分析圖Note：In Figures A to E,the 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；F，Correlation diagram of validated genes

本研究中，對蘆竹在0 mmol·L-1，50 mmol·L-1，100 mmol·L-1，150 mmol·L-1和200 mmol·L-1NaCl鹽溶液下的超氧化物歧化酶、過氧化物酶、過氧化氫酶活性，脯氨酸、可溶性糖、丙二醛含量和Fv/F0指標進行測定，結果顯示隨著鹽濃度的增加，脯氨酸含量在200 mmol·L-1處達到最高值。可溶性糖含量在150 mmol·L-1處達到了最高，200 mmol·L-1處下降到最低值。而其他生理指標變化并不明顯，說明其對鹽環境可能并不是十分敏感，并沒有對其造成較大的損傷，具有一定的耐受性，也可能是由于脅迫時間與濃度不足導致的。

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

為了確定蘆竹響應鹽脅迫的相關基因，利用DESeq2進行基因的差異表達分析。結果顯示這些差異表達基因主要與硝酸鹽、一氧化氮、活性氮、活性氧、類黃酮生物合成、芪類與二芳基庚烷類和姜酚生物合成、光合作用、玉米素生物合成和硒化合物代謝等條目相關。這些過程可能在蘆竹響應鹽脅迫的過程中發揮著重要作用。而414_c1_g2_i5(TPC1)，4389_c1_g1_i14，4238_c0_g1_i12，0_c1_g3_i2，8776_c0_g1_i37(GSH2)和37109_c0_g1_i12可能是蘆竹響應鹽脅迫過程中的關鍵基因。414_c1_g2_i5(TPC1)是從植物克隆的第一個TPC通道[23]，定位于液泡膜，負責產生慢速液泡(SV) 電流，參與各種生理過程的調節，例如發芽、氣孔開放、茉莉酸的生物合成以及高鹽濃度誘導的長距離鈣信號傳播[24]。4389_c1_g1_i14在細胞分化，細胞壁生物發生，化學穩態，發育生長等方面起作用；4238_c0_g1_i12是與AtMSU81 同源的無功能假基因；0_c1_g3_i2是富含亮氨酸的重復蛋白激酶家族的成員，可以與FRD3相互作用，參與檸檬酸鹽外流運輸；8776_c0_g1_i37(GSH2)編碼谷胱甘肽生物合成酶，三肽谷胱甘肽(GSH)是真核細胞中主要的非蛋白質硫醇化合物，參與細胞氧化還原狀態的調節，這在植物暴露于氧化應激期間尤為重要，已被證實能夠提高抵抗逆境脅迫的能力[25-26]。37109_c0_g1_i12為HAD超家族，IIIB酸性磷酸酶亞家族成員。

結合生理生化指標對差異基因進行WGCNA分析，以探求與生理變化對應的相關基因，從不同角度研究蘆竹的核心鹽響應基因。與Pro相關的turquoise，green模塊的hub基因編碼核堿基抗壞血酸轉運蛋白和外膜孔蛋白等相關蛋白質，還包括一些具有某些特定結構的蛋白例如RING/U-box超家族蛋白和RRM結構域蛋白等等。其中核堿基抗壞血酸轉運蛋白在植物中運輸黃嘌呤和尿酸，已有相關研究證明黃嘌呤和尿酸在緩解鹽脅迫方面具有潛在用途，并且過表達核堿基抗壞血酸轉運蛋白相關基因可增強植物的耐鹽性[27]。而一些與DUF1645[28]，RING/U-box[29]，RRM[30]和MYB[31]相關的基因家族也被證實在植物應對非生物脅迫過程中發揮著重要作用。

4 結論

本文綜合生理學與轉錄組學來探究蘆竹對不同程度鹽脅迫的響應變化。結果表明隨著鹽濃度的增大，蘆竹脯氨酸和可溶性糖的含量會發生顯著變化。但其對鹽環境并不是十分敏感，具有一定的耐受性。而在響應鹽脅迫的過程中，差異表達基因主要被富集到了與硝酸鹽、一氧化氮、活性氮、活性氧和類黃酮生物合成相關的通路。一些變化倍數較大的核心差異表達基因和一些與脯氨酸含量變化相關的基因可能具有重要作用，可作為候選的耐鹽基因。本研究結果可為蘆竹耐鹽的分子機制研究提供線索，也為今后蘆竹的分子育種提供參考。