楊 闖 王 玲 全成滔 余良倩 戴 成 郭 亮 傅廷棟 馬朝芝

華中農(nóng)業(yè)大學作物遺傳改良全國重點實驗室 / 國家油菜工程技術(shù)研究中心 / 洪山實驗室, 湖北武漢 430070

土壤鹽漬化是威脅全球農(nóng)業(yè)生產(chǎn)的主要環(huán)境挑戰(zhàn)之一[1]。世界上大約有20%的灌溉農(nóng)田會受到土壤鹽漬化的不利影響[2]。自然環(huán)境惡化、灌溉方式不當以及氣候變化等因素也正在進一步加劇土壤鹽漬化問題[3]。預計到2050年, 鹽漬化土壤面積將達到全球可利用耕地的50%以上[4]。近些年來, 我國的農(nóng)業(yè)耕地也在不斷受到鹽漬化侵蝕, 鹽漬地總面積已達到3.69×107hm2, 給作物生產(chǎn)和農(nóng)業(yè)可持續(xù)發(fā)展帶來了極大的挑戰(zhàn)[5-6]。油菜是世界四大油料作物之一, 具有重要的油用、飼用、菜用以及旅游觀光價值[7]。作為食用植物油和植物蛋白的主要來源, 油菜在我國農(nóng)產(chǎn)品中占有重要地位, 對于保障我國油料供給安全尤為重要[8]。然而, 鹽脅迫目前是影響油菜生長發(fā)育的主要環(huán)境條件之一, 可能會造成油菜減產(chǎn)、品質(zhì)下降甚至死亡。因此, 解析甘藍型油菜響應(yīng)鹽脅迫的機制, 并從中挖掘油菜響應(yīng)鹽脅迫的核心基因, 對于油菜的耐鹽性改良具有重要意義。

鹽脅迫對植物的毒害作用主要包括滲透脅迫和離子脅迫2種[9]。鹽脅迫的滲透作用會導致土壤溶液的水勢低于根細胞水勢, 從而抑制根系的吸水作用; 此外, 滲透脅迫還會導致葉片氣孔關(guān)閉從而使光合效率下降, 最終影響植物的正常生長發(fā)育。離子脅迫主要是由細胞中的Na+和Cl-過量積累導致的。Na+的毒性主要是通過抑制新陳代謝途徑中關(guān)鍵酶活性從而影響植物的生長發(fā)育; 同時,過量的Na+也會影響植物對氮、磷、鉀、鈣等其他礦質(zhì)元素的吸收利用。由于植物對NO3-、SO42-以及Cl-的吸收都是通過相同的非選擇性陰離子通道, 因此過量的Cl-可能會導致氮素以及硫素的缺乏。

植物為了響應(yīng)鹽脅迫, 也進化出了互相關(guān)聯(lián)的調(diào)節(jié)通路并引起細胞內(nèi)的廣泛變化, 其中許多變化是適應(yīng)性反應(yīng), 可導致抗逆性增加, 這也是植物在長期進化過程中保留下來的適應(yīng)性機制[10]。鹽脅迫早期響應(yīng)階段, 植物可以通過高滲感受器OSCA1以及低滲感受器MSL8感受滲透壓的變化[11-12], 同時植物還可以通過鈉離子感受器MOCA1感受鹽濃度變化[13]。植物感受到鹽脅迫后, 可通過滲透調(diào)節(jié)、組織特異性離子轉(zhuǎn)運、激素信號傳導、以及活性氧清除等機制來應(yīng)對鹽脅迫, 其中也涉及到一系列基因的表達, 包括相關(guān)功能基因和調(diào)節(jié)基因, 如轉(zhuǎn)錄因子、蛋白激酶等。SOS信號通路是近20年來發(fā)現(xiàn)的典型的鹽脅迫響應(yīng)通路,SOS1、SOS2以及SOS3編碼構(gòu)成SOS信號通路的蛋白質(zhì), 在多種植物的鹽脅迫響應(yīng)過程中發(fā)揮著重要功能[14]。AtNHX1和AtNHX2編碼液泡膜上的Na+/H+逆向轉(zhuǎn)運蛋白, 主要在根和葉中表達, 并選擇性地將Na+運輸?shù)揭号? 其突變體對鹽脅迫更加敏感[15]。其中,轉(zhuǎn)錄因子作為關(guān)鍵的調(diào)控因子, 其介導的基因調(diào)控網(wǎng)絡(luò)在植物響應(yīng)鹽脅迫時起著不可或缺的作用。擬南芥AtNIG1是植物中首先被鑒定出來參與耐鹽性的BHLH轉(zhuǎn)錄因子家族成員, 過表達AtNIG1表現(xiàn)出更高的耐鹽性[16]。玉米轉(zhuǎn)錄因子ZmWRKY114作為負調(diào)控因子, 通過ABA依賴性途徑參與鹽脅迫響應(yīng)[17]。水稻轉(zhuǎn)錄因子基因OsDOF15通過介導乙烯生物合成途徑, 在鹽脅迫抑制初生根伸長的過程中起著關(guān)鍵作用[18]。棉花轉(zhuǎn)錄因子基因GhABF3可以通過降低細胞氧化程度顯著提高陸地棉的耐鹽性[19]。大豆轉(zhuǎn)錄因子基因GmNAC06可以通過控制毛狀根Na+/K+比值維持離子穩(wěn)態(tài), 從而增強大豆的耐鹽性[20]。這些研究結(jié)果表明, 相關(guān)轉(zhuǎn)錄因子的鑒定及功能解析對于理解植物響應(yīng)鹽脅迫的分子機制至關(guān)重要。

隨著高通量測序技術(shù)的快速發(fā)展, RNA-seq技術(shù)已經(jīng)被廣泛地用于多種植物的研究中, 同時也為轉(zhuǎn)錄組學研究提供了大量的分析方法。權(quán)重基因共表達網(wǎng)絡(luò)分析(WGCNA)是最常用的分析RNA-seq數(shù)據(jù)的方法之一。相比于其他分析方法, WGCNA在產(chǎn)生具有生物學意義的結(jié)果方面也更加穩(wěn)健[21]。該方法可以通過構(gòu)建基因共表達網(wǎng)絡(luò)將上萬個基因劃分成十幾個具有相似表達模式的基因模塊, 并通過分析與性狀高度相關(guān)的關(guān)鍵模塊從中挖掘核心基因。Zheng等[22]利用324個RNA-Seq數(shù)據(jù)構(gòu)建了玉米表皮蠟質(zhì)合成相關(guān)基因的共表達網(wǎng)絡(luò), 發(fā)現(xiàn)幾乎所有已知的有關(guān)蠟質(zhì)合成的基因都在同一個模塊, 這為后續(xù)深入挖掘玉米表皮蠟質(zhì)合成相關(guān)基因提供了重要的參考依據(jù)。李旭凱等[23]利用WGCNA預測出與多種脅迫高度相關(guān)的25個抗逆關(guān)鍵基因, 為水稻的多抗逆基因挖掘以及利用提供了參考依據(jù)。Ye等[24]利用鹽梯度處理的RNA-seq數(shù)據(jù), 闡明了互花米草高耐鹽性的轉(zhuǎn)錄組動力學, 并通過WGCNA分析發(fā)現(xiàn)SaOST1等蛋白激酶編碼基因是耐鹽調(diào)控網(wǎng)絡(luò)中的核心基因。Zhao等[25]以花生為材料, 利用WGCNA篩選出一個與干旱脅迫生理數(shù)據(jù)密切相關(guān)的重要模塊, 并從中預測到9個抗旱候選基因。Yang等[26]基于剪接異構(gòu)體的WGCNA分析, 篩選出與不同非生物脅迫相關(guān)的23個共表達模塊, 并發(fā)現(xiàn)許多重要的核心基因在單一應(yīng)激或多種應(yīng)激反應(yīng)中起關(guān)鍵作用, 其中主要包括轉(zhuǎn)錄因子、RNA結(jié)合蛋白等。

利用WGCNA構(gòu)建基因共表達網(wǎng)絡(luò), 已經(jīng)被廣泛地應(yīng)用于多種植物的轉(zhuǎn)錄組研究以及核心基因篩選。本研究通過測得的90份RNA-seq數(shù)據(jù), 獲得了甘藍型油菜鹽脅迫處理前后的高分辨率時間動態(tài)轉(zhuǎn)錄表達譜, 進一步通過差異基因分析和WGCNA分別構(gòu)建油菜根系以及葉片組織響應(yīng)鹽脅迫的基因共表達網(wǎng)絡(luò), 并從中挖掘油菜響應(yīng)鹽脅迫的核心基因。研究結(jié)果可為甘藍型油菜的鹽脅迫共表達網(wǎng)絡(luò)研究以及油菜的耐鹽性分子育種提供理論支持。

1 材料與方法

1.1 植物材料和鹽脅迫處理

選取甘藍型油菜ZS11作為試驗材料, 采用水培法進行培養(yǎng), 培養(yǎng)條件為: 光照16 h、黑暗8 h, 溫度24℃左右, 相對濕度為75%。培養(yǎng)27 d后, 選取長勢一致的幼苗,使用200 mmol L–1氯化鈉模擬鹽脅迫處理, 分別在鹽脅迫處理后0、0.25、0.5、1、3、6、12、24 h進行取樣, 同時各個時間點均用正常培養(yǎng)條件下的幼苗作為對照組, 取樣時分離根部組織和葉片組織, 立即將其包入錫箔紙并投入–80℃的液氮中儲存, 每個處理3個生物學重復。取樣完畢后送公司進行RNA提取以及轉(zhuǎn)錄組測序。鹽脅迫處理條件如圖1所示, 其中0 h的樣本組未顯示。本研究所使用的原始數(shù)據(jù)尚未公開發(fā)表, 但單個基因受到鹽脅迫處理后的轉(zhuǎn)錄表達量已上傳并可通過BnIR[27]網(wǎng)站(http://yanglab.hzau.edu.cn/BnIR/eFP_single_gene)查詢。

圖1 鹽脅迫樣品處理條件Fig. 1 Treatment conditions of salt stress samples

1.2 RNA-seq數(shù)據(jù)處理與分析

1.2.1 RNA-seq數(shù)據(jù)的質(zhì)控與定量 使用FastQC和MultiQC軟件[28]對RNA-seq數(shù)據(jù)進行質(zhì)控, 使用Trim Galore過濾掉低質(zhì)量的reads片段, 獲得clean reads用于后續(xù)分析。使用Hisat2[29]將其批量比對到ZS11參考基因組, 統(tǒng)計所有樣本的比對率。使用Samtools[30]將sam文件轉(zhuǎn)換成bam文件, 使用FeatureCounts[31]對回帖到參考基因上的CDS序列進行計數(shù), 獲得原始counts值。根據(jù)counts值以及CDS序列長度, 使用TBtools[32]計算出基因的TPM值。同時使用原始counts值, 利用R包DESeq2[33]進行標準化, 對鹽脅迫根系以及葉片組織在各個時間點的處理組樣本與其對照組樣本進行兩兩差異基因分析,最終獲得組織在各個處理時間點的差異基因以及上下調(diào)基因數(shù)目。篩選差異基因的閾值設(shè)置為|log2(Fold Change)|>1 & FDR<0.05。

1.2.2 利用WGCNA構(gòu)建基因共表達網(wǎng)絡(luò) WGCNA(weighted gene co-expression network analysis), 即權(quán)重基因共表達網(wǎng)絡(luò)分析, 可以通過構(gòu)建共表達網(wǎng)絡(luò)尋找協(xié)同表達的基因模塊, 探究基因模塊與樣本性狀之間的相關(guān)性, 從中挖掘核心基因, 經(jīng)常適用于非生物逆境脅迫不同時間點應(yīng)答、病原菌侵染后不同時間點應(yīng)答等研究方向[34]。本研究考慮到低表達基因一般不具備生物學意義, 將Deseq2分析得到的鹽脅迫組織響應(yīng)差異基因賦予TPM值,進一步篩選得到至少在3個樣本中TPM值大于1的組織響應(yīng)差異表達基因用于WGCNA的輸入值。最小模塊的基因數(shù)不低于100個, 其余參數(shù)按照默認設(shè)置。使用R包Cluster Profiler[35]對模塊進行GO富集分析, 閾值設(shè)置為FDR<0.05。本研究定義模塊MEs與樣本時間點之間的相關(guān)性絕對值大于0.6,P值小于0.05為顯著相關(guān); 以權(quán)重值(weight)大于0.2作為基因間的有效連接, 連通度前10%的基因為核心基因。使用Cytoscape軟件對關(guān)鍵模塊的共表達網(wǎng)絡(luò)進行可視化分析[36]。

2 結(jié)果與分析

2.1 鹽脅迫轉(zhuǎn)錄組數(shù)據(jù)的產(chǎn)生及初步驗證

本研究于鹽脅迫處理0、0.25、0.5、1、3、6、12、24 h對油菜葉片組織和根系組織進行轉(zhuǎn)錄組測序, 獲得90份RNA-seq數(shù)據(jù)。通過質(zhì)檢和過濾, 共獲得52.78億條clean reads。將其比對到ZS11參考基因組, 發(fā)現(xiàn)所有樣本的比對率都超過95%, 表明本研究所獲得的RNA-seq數(shù)據(jù)具備較高質(zhì)量。為了驗證鹽脅迫處理的有效性, 利用擬南芥已報道參與鹽脅迫響應(yīng)的marker基因, 觀察其對應(yīng)的油菜同源基因在鹽脅迫處理前后的表達模式。SOS1、WRKY33通常在鹽脅迫早期響應(yīng)階段被激活[37-38], 油菜根系組織中的BnaC09.SOS1以及BnaC04.WRKY33在鹽脅迫1 h內(nèi)表達量顯著上升(圖2-A, B)。HKT1作為Na+運載體, 可以裝載葉片中的鈉離子到韌皮部再循環(huán)到根部[39],油菜葉片組織中的BnaA09.HKT1在鹽脅迫處理3 h時出現(xiàn)顯著差異表達(圖2-C)。ABI1可以負調(diào)控ABA促進氣孔關(guān)閉從而影響植物耐鹽性[40], 油菜葉片以及根系組織中的BnaA03.ABI1在鹽脅迫處理3 h時表達量顯著上升(圖2-D)。以上結(jié)果均與擬南芥中的報道結(jié)果類似, 這表明本研究對甘藍型油菜幼苗進行的鹽脅迫處理是有效的,也進一步驗證了RNA-seq數(shù)據(jù)的可靠性。

圖2 鹽脅迫響應(yīng)marker基因的表達模式Fig. 2 Relative expression pattern of marker genes under salt stress

2.2 鹽脅迫的相關(guān)性分析

為了探究樣本處理之間的相關(guān)性, 本研究進行了主坐標分析以及層次聚類分析。從鹽脅迫的主坐標分析結(jié)果來看, 第1主坐標主要解釋的是組織表達的差異性, 其貢獻度較大, 約占75.08%, 表明后續(xù)在利用WGCNA構(gòu)建基因表達網(wǎng)絡(luò)時, 應(yīng)將脅迫的根系和葉片組織分開分析,避免組織差異性覆蓋脅迫處理差異性; 第2主坐標主要解釋的是對照組以及鹽脅迫處理之間的差異性, 鹽脅迫處理3 h、6 h、12 h、24 h的樣本表達與對照組形成明顯差異,這也表明脅迫處理的有效性和數(shù)據(jù)的良好性(圖3-A)。從層次聚類結(jié)果來看, 根和葉由于組織表達的差異性較大聚為2類, 這與主坐標分析結(jié)果相對一致。另外, 無論是根系組織還是葉片組織, 1 h處理前后的組織樣本均形成了明顯的聚類差異。根據(jù)組織響應(yīng)鹽脅迫的樣本聚類結(jié)果進行初步劃分, 可以將鹽脅迫處理1 h前后劃分為早期響應(yīng)階段以及后期響應(yīng)階段(圖3-B)。

圖3 鹽脅迫的相關(guān)性分析Fig. 3 Correlation analysis of salt stress

2.3 差異基因分析

本研究利用DESeq2包進行差異基因分析, 獲得了油菜根系以及葉片組織在鹽脅迫各個處理時間點的差異基因以及上下調(diào)差異基因數(shù)目, 并分析了差異基因數(shù)目隨鹽脅迫處理時間的變化趨勢。鹽脅迫處理1 h內(nèi), 葉片組織在0.25 h、0.5 h和1 h的差異基因數(shù)目始終小于根系組織, 表明根系組織在早期階段的響應(yīng)程度大于葉片組織; 鹽脅迫處理1 h后, 葉片組織在任意一個處理時間點的差異基因數(shù)目都大于根系組織, 表明葉片組織在后期階段的響應(yīng)程度大于根系組織(圖4-A)。其中, 根系組織在鹽脅迫處理3 h時達到最高的響應(yīng)程度, 此時差異基因為10,883個; 而葉片組織在處理24 h時達到最高的響應(yīng)程度, 此時差異基因數(shù)目為18,414個。進一步分析鹽脅迫上下調(diào)差異基因數(shù)目的變化趨勢, 發(fā)現(xiàn)鹽脅迫處理24 h內(nèi), 根系組織中的上調(diào)基因數(shù)目始終大于下調(diào)基因, 而葉片組織僅在處理3 h內(nèi)的上調(diào)基因數(shù)目大于下調(diào)基因, 而隨著鹽脅迫時間的增長,葉片組織中的上調(diào)基因數(shù)開始降低并逐漸少于下調(diào)基因(圖4-B)。表明油菜幼苗對鹽脅迫的響應(yīng)程度存在組織差異性。鹽脅迫處理12 h內(nèi), 根系組織和葉片組織的響應(yīng)程度均呈現(xiàn)出先升高后下降的變化趨勢, 變化拐點均出現(xiàn)在處理3 h附近, 表明油菜幼苗對鹽脅迫的響應(yīng)程度也存在組織相似性。本研究定義至少在任意一個時間點存在顯著差異的基因作為鹽脅迫組織響應(yīng)差異基因集, 通過去重取并集, 得到根系組織以及葉片組織響應(yīng)差異基因分別是20,462個和29,334個, 表明全程鹽脅迫處理24 h內(nèi), 葉片組織整體上的響應(yīng)程度比根系更劇烈。

圖4 鹽脅迫處理24 h內(nèi)的差異基因數(shù)目及變化趨勢Fig. 4 Number and trend of differential genes under salt stress treatment for 24 hours

2.4 權(quán)重基因共表達網(wǎng)絡(luò)的構(gòu)建

為了構(gòu)建油菜幼苗響應(yīng)鹽脅迫的基因共表達網(wǎng)絡(luò),并從中挖掘核心基因, 基于鹽脅迫處理24 h的轉(zhuǎn)錄組表達矩陣用于WGCNA分析。同時, 為了防止組織差異性過大影響分析結(jié)果, 在保證不低于WGCNA所要求的樣本數(shù)的情況下, 分別構(gòu)建了根系組織以及葉片組織響應(yīng)鹽脅迫的基因共表達網(wǎng)絡(luò)。過濾低表達基因, 篩選26,839個葉片組織差異表達基因, 使用動態(tài)樹切割方法, 選取最佳軟閾值β=14構(gòu)建葉片組織響應(yīng)鹽脅迫的基因共表達網(wǎng)絡(luò)(圖5-A), 共識別出14個不同的共表達模塊, 最大和最小模塊分別是turquoise和tan模塊, 分別含有8651個和260個基因(附表1)。過濾低表達基因, 篩選18,952個根系組織響應(yīng)差異表達基因, 選取最佳軟閾值β=18構(gòu)建根系組織響應(yīng)鹽脅迫的基因共表達網(wǎng)絡(luò)(圖5-B), 共識別出11個不同的共表達模塊, 最大和最小模塊分別是turquoise以及greenyellow模塊, 分別含有4058個和234個基因(附表2)。

圖5 最佳軟閾值的選擇Fig. 5 Selection of optimal soft threshold (power)

2.5 早期及后期響應(yīng)模塊的鑒定

在葉片組織響應(yīng)鹽脅迫的共表達網(wǎng)絡(luò)中, 通過計算模塊與樣本處理時間點之間的相關(guān)性發(fā)現(xiàn), tan模塊、salmon模塊分別與鹽脅迫早期響應(yīng)階段0.5 h、1 h顯著正相關(guān); red模塊、cyan模塊、turquiose模塊以及green模塊分別與鹽脅迫后期響應(yīng)階段3 h、6 h、12 h、24 h顯著正相關(guān)(圖6-A)。在根系組織響應(yīng)鹽脅迫的共表達網(wǎng)絡(luò)中,通過計算模塊與樣本處理時間點之間的相關(guān)性發(fā)現(xiàn), yellow模塊與鹽脅迫早期響應(yīng)階段0.25 h以及0.5 h顯著正相關(guān); blue模塊和turquiose模塊分別與鹽脅迫處理1 h、3 h顯著正相關(guān); black模塊與鹽脅迫處理12 h顯著正相關(guān);purple模塊以及red模塊均與鹽脅迫處理24 h顯著正相關(guān)(圖6-B)。綜合考慮模塊與樣本處理時間點的相關(guān)性、模塊內(nèi)的轉(zhuǎn)錄因子所占比例以及模塊的GO富集結(jié)果, 選取葉片組織中的tan模塊和根系組織中的yellow模塊作為組織早期響應(yīng)鹽脅迫的關(guān)鍵模塊, 葉片組織中的green模塊和根系組織的red模塊作為組織后期響應(yīng)鹽脅迫的關(guān)鍵模塊。葉片組織中的早期響應(yīng)tan模塊包含260個基因, 其中有75個轉(zhuǎn)錄因子, 該模塊內(nèi)的絕大多數(shù)基因在鹽脅迫處理0.5 h的表達量顯著高于其他時間點(圖7-A)。根系組織中的早期響應(yīng)yellow模塊包含2322個基因, 其中含有276個轉(zhuǎn)錄因子, 該模塊內(nèi)的絕大多數(shù)基因在鹽脅迫處理0.25 h以及0.5 h顯著表達(圖7-B)。葉片中的后期響應(yīng)green模塊包含1654個基因, 其中含有138個轉(zhuǎn)錄因子;根系中的后期響應(yīng)red模塊包含1238個基因, 含有99個轉(zhuǎn)錄因子, 2個模塊內(nèi)的絕大多數(shù)基因均在鹽脅迫處理24 h時的表達量顯著高于其他時間點(圖7-C, D)。

圖6 基因模塊與樣本時間點的相關(guān)性熱圖Fig. 6 Correlation heat map between gene modules and sample time point

圖7 早期及后期響應(yīng)模塊的基因表達譜Fig. 7 Gene expression profiles of early and late response modules

2.6 早期及后期響應(yīng)模塊的GO富集分析

為了分析組織響應(yīng)鹽脅迫的共表達模塊特征, 本研究分別對早期響應(yīng)階段的tan模塊和yellow模塊, 后期響應(yīng)階段的green和red模塊進行GO富集分析。早期響應(yīng)階段, 根系最先受到鹽脅迫的刺激, 根系組織中的yellow模塊富集到了對生物刺激反應(yīng)的調(diào)控、對外部刺激反應(yīng)的調(diào)控、磷酸化信號轉(zhuǎn)導、對鹽的響應(yīng)等生物學過程(圖8-B)。葉片中的tan模塊富集到與葉片衰老、鹽脅迫響應(yīng)的正向調(diào)控、自噬調(diào)節(jié)等生物學過程(圖8-A)。通過對比分析發(fā)現(xiàn), 早期響應(yīng)模塊均富集到了脫落酸、乙烯、茉莉酸等激素信號通路, 這可能表明油菜幼苗早期響應(yīng)鹽脅迫就會有多種激素參與全身器官的系統(tǒng)性反應(yīng)。早期響應(yīng)階段,葉片中還富集到了對細菌的防御反應(yīng)和對昆蟲的防御反應(yīng), 根系也富集到了對細菌防御反應(yīng)的調(diào)控以及對細菌來源分子的響應(yīng), 這可能表明油菜幼苗早期響應(yīng)鹽脅迫可能與早期響應(yīng)細菌等生物脅迫有相似之處。除此之外,葉片中也富集到了抗毒素合成通路, 根系中富集到了類黃酮生物合成過程的調(diào)控。這些早期響應(yīng)階段的生物學通路富集, 可能表明油菜幼苗在早期響應(yīng)鹽脅迫時也會通過合成抗毒素等物質(zhì)來響應(yīng)鹽脅迫。鹽脅迫后期響應(yīng)階段,green和red模塊顯著富集到了對過氧化氫的響應(yīng)、細胞對鹽脅迫的反應(yīng)、對脫水反應(yīng)的調(diào)控等鹽脅迫響應(yīng)相關(guān)通路(圖8-C, D)。其中, 葉片中的green模塊也富集到了與早期響應(yīng)相似的生物學通路, 如脫落酸激活信號通路;green模塊也特異富集到了有機羥基化合物合成、肌醇生物合成等通路(圖8-C)。根系中的red模塊特異富集到了有機酸轉(zhuǎn)運、植物型次生細胞壁的發(fā)生、苯丙烷生物合成過程、硫代葡萄糖苷代謝過程、木質(zhì)素代謝過程等通路(圖8-D)。總體來看, 無論是早期響應(yīng)模塊, 還是后期響應(yīng)模塊, 都富集到了鹽脅迫響應(yīng)的相關(guān)通路, 表明利用WGCNA可以構(gòu)建具有生物學意義的共表達模塊, 這些模塊可成為本研究的重點關(guān)注對象。

(圖8)

圖8 早期及后期響應(yīng)模塊的GO富集分析Fig. 8 GO enrichment of early and late response modules

2.7 核心基因的鑒定與分析

為了識別關(guān)鍵模塊中的核心基因,利用cytoscape進行共表達網(wǎng)絡(luò)可視化。以權(quán)重值大于0.2作為閾值, 篩選連通度最大的前10%基因作為核心基因, 早期響應(yīng)tan模塊篩選到23個核心基因, 包含15個核心轉(zhuǎn)錄因子(圖9-A)。其中連通度最高的3個核心基因分別為BnaA04G0009600ZS(AGP20)、BnaC04G0304800ZS(CAMBP25)、BnaA10G 0150300ZS(RING1)。BnaA04G0009600ZS(AGP20)編碼調(diào)節(jié)細胞壁擴張的阿拉伯半乳糖蛋白, 其表達主要由鹽脅迫誘導[41]。早期響應(yīng)yellow模塊篩選到221個核心基因,包含26個核心轉(zhuǎn)錄因子 (圖9-B)。其中連通度最高的3個核心基因分別為BnaC09G0584800ZS(NHL3)、BnaC08G 0479500ZS(WAKL8)、BnaA09G0091000ZS(SUPA)。BnaA09G 0091000ZS對應(yīng)的擬南芥基因SUPA定位于過氧化物酶體中, 鹽脅迫處理15 min內(nèi)迅速上調(diào)表達, 過表達可導致擬南芥中活性氧類物質(zhì)積累增加以及抗逆性改變[42]。后期響應(yīng)green模塊篩選到核心基因137個, 包含14個核心轉(zhuǎn)錄因子(圖9-C)。其中連通度最高的3個核心基因分別為BnaA03G0371900ZS(AOX1A)、BnaC07G0256700ZS(PBL19)、BnaA02G0027300ZS(DGK1)。BnaA03G 0371900ZS對應(yīng)的擬南芥基因AOX1A編碼線粒體電子傳遞鏈的替代氧化酶, 會限制脯氨酸誘導的氧化應(yīng)激, 從而利于擬南芥鹽度恢復[43]。后期響應(yīng)red模塊篩選到核心基因115個, 包含12個核心轉(zhuǎn)錄因(圖9-D)。其中連通度最高的3個核心基因分別為BnaC05G0203500ZS(PME6)、BnaA02G0180600ZS(WRKY57)、BnaC06G0284600ZS(KTI1)。BnaC05G0203500ZS對應(yīng)的擬南芥PME6編碼果膠甲酯酶,在保衛(wèi)細胞中高度表達, 并且是氣孔功能所必需的, 基因恢復可以挽救保衛(wèi)細胞壁果膠甲酯化狀態(tài)、氣孔功能和植物生長[44], 該基因也在油菜鹽脅迫后期顯著差異表達,可能通過氣孔調(diào)節(jié)應(yīng)對鹽脅迫。

圖9 早期及后期響應(yīng)模塊的基因共表達網(wǎng)絡(luò)可視化Fig. 9 Co-expression network visualization of early and late response modules

轉(zhuǎn)錄因子作為關(guān)鍵調(diào)控因子, 其介導的基因調(diào)控網(wǎng)絡(luò)經(jīng)常在植物響應(yīng)鹽脅迫時發(fā)揮著重要作用。本研究著重關(guān)注核心轉(zhuǎn)錄因子, 通過雙向blast將核心基因比對至擬南芥, 并借助TAIR數(shù)據(jù)庫注釋核心轉(zhuǎn)錄因子的功能(附表3), 發(fā)現(xiàn)油菜同源基因WRKY33、DDF1、AZF2、ARR1、ZFHD1、DREB2B等是擬南芥中已知參與編碼鹽脅迫響應(yīng)的核心轉(zhuǎn)錄因子。其中, 擬南芥中的同源基因WRKY33對應(yīng)油菜中的多個拷貝BnaA04G0247200ZS、BnaC03G 0217300ZS、BnaA03G0185200ZS、BnaC04G0562300ZS, 4個拷貝均在早期響應(yīng)tan模塊中占據(jù)較高連通度, 并參與多種非生物脅迫響應(yīng), 尤其是對鹽脅迫的響應(yīng)[37]。Yellow模塊中的早期響應(yīng)基因BnaA06G0081500ZS、BnaC05G 0100700ZS在擬南芥中的同源基因DDF1編碼ERF/AP2轉(zhuǎn)錄因子家族DREB亞家族成員, 過表達可增強對高水平鹽的耐受性[45]。Yellow模塊中的早期響應(yīng)基因BnaC05G 0398500ZS在擬南芥中的同源基因AZF2編碼鋅指蛋白轉(zhuǎn)錄因子, 響應(yīng)ABA、高鹽和輕度脫水反應(yīng)[46]。Green模塊中的后期響應(yīng)基因BnaA01G0191600ZS在擬南芥中的同源基因ARR1屬于細胞分裂素信號成分擬南芥響應(yīng)調(diào)節(jié)因子1, ARR1蛋白可以通過MPK3/6負調(diào)控作用保持穩(wěn)定從而增強耐鹽性[47]。Red模塊中的后期響應(yīng)基因BnaA07G 0311300ZS在擬南芥中的同源基因ZFHD1編碼鋅指轉(zhuǎn)錄因子家族成員, 可由干旱、高鹽和脫落酸誘導表達[48]。無論是早期響應(yīng)模塊還是后期響應(yīng)模塊, 均有已知的核心轉(zhuǎn)錄因子參與擬南芥的鹽脅迫響應(yīng), 這表明本研究所篩選的關(guān)鍵模塊以及核心轉(zhuǎn)錄因子具有較高的可靠性。據(jù)此初步推測, 這些核心基因可能也在甘藍型油菜的耐鹽性過程中發(fā)揮著重要作用。

2.8 核心轉(zhuǎn)錄因子在505份鹽脅迫油菜群體變異數(shù)據(jù)中的SNPs及單倍型分析

進一步利用已發(fā)表的505份鹽脅迫處理油菜群體變異數(shù)據(jù)庫, 分析核心轉(zhuǎn)錄因子在甘藍型油菜群體變異中的SNPs及單倍型情況(http://yanglab.hzau.edu.cn/BnIR/single_locus)[49]。結(jié)果表明, 多個核心轉(zhuǎn)錄因子在505份油菜群體中存在較為豐富的SNPs變異以及單倍型類型。其中,BnaC04G0015700ZS(BnWRKY46)是根系組織早期響應(yīng)核心基因, 其在鹽脅迫處理0.5 h時可上調(diào)表達約6倍。BnaA07G0271900ZS(BnWRKY57)是根系組織后期響應(yīng)核心基因, 其在鹽脅迫后期處理24 h可達到最大上調(diào)倍數(shù)約7倍。BnWRKY46和BnEWRKY57均由2個內(nèi)含子和3個外顯子組成, 經(jīng)在群體內(nèi)序列變異分析發(fā)現(xiàn), 這2個核心基因均呈現(xiàn)出較為豐富的SNPs變異(圖10-A, C)。單倍型分析發(fā)現(xiàn),BnWRKY46可初步劃分為3種單倍型, 不同單倍型在鹽脅迫處理后的根長變短, 其中hap.B和hap.C型在鹽脅迫處理前后的根長相對值顯著小于hap.A型(圖10-B)。BnWRKY57可劃分為3種單倍型, 不同單倍型的根長在鹽脅迫處理后顯著縮短, 其中hap.C型在鹽脅迫處理前后的根長相對值顯著大于hap.A和hap.B型(圖10-D)。表明本研究所挖掘到的核心轉(zhuǎn)錄因子中很有可能存在參與鹽脅迫響應(yīng)的關(guān)鍵候選基因, 深入鑒定鹽脅迫相關(guān)的極端單倍型也將有利于為耐鹽型油菜提供優(yōu)質(zhì)的種質(zhì)資源。

圖10 核心基因BnWRKY46和BnWRKY57在505份群體數(shù)據(jù)中的SNPs及單倍型分析Fig. 10 SNPs and haplotype analysis of hub genes BnWRKY46 and BnWRKY57 in 505 population data

3 討論

3.1 甘藍型油菜響應(yīng)鹽脅迫的時間動態(tài)轉(zhuǎn)錄表達譜分析

油菜作為世界上重要的油料作物, 在生長過程中經(jīng)常受到鹽脅迫的影響, 可能會導致減產(chǎn)、品質(zhì)下降甚至死亡。而甘藍型油菜被認為是中等耐鹽作物, 因此解析甘藍型油菜響應(yīng)鹽脅迫的調(diào)控網(wǎng)絡(luò), 對于油菜的耐鹽性改良以及鹽堿地的開發(fā)和利用都具有重要意義。本研究利用轉(zhuǎn)錄組測序獲得了甘藍型油菜幼苗在鹽脅迫處理前后的90份RNA-seq數(shù)據(jù), 將其比對到ZS11參考基因組, 發(fā)現(xiàn)所有樣本的比對率都超過95%, 鹽脅迫響應(yīng)marker基因的表達模式也與擬南芥報道結(jié)果類似, 表明鹽脅迫處理是有效的, 數(shù)據(jù)質(zhì)量是可靠的。

由于差異基因數(shù)目可以在一定程度上反映脅迫響應(yīng)程度, 本研究發(fā)現(xiàn)油菜幼苗對鹽脅迫的響應(yīng)程度存在組織差異性和相似性。全程鹽脅迫處理24 h內(nèi), 葉片組織整體上的響應(yīng)程度比根系更劇烈。根系組織24 h內(nèi)響應(yīng)鹽脅迫的差異基因動態(tài)變化趨勢與2007年發(fā)表的擬南芥鹽脅迫24 h內(nèi)的轉(zhuǎn)錄動態(tài)極為類似, 均呈現(xiàn)出先上升后下降的趨勢。并且, 根系組織在24 h內(nèi)的任意處理時間點的上調(diào)基因數(shù)目始終大于下調(diào)基因數(shù)目, 這一點也與擬南芥中的研究結(jié)果相對一致[50]。然而, 不同之處在于, 油菜中的組織響應(yīng)差異基因數(shù)目遠高于擬南芥, 這可能是由于異源四倍體的油菜經(jīng)歷過自然加倍等事件具有更多的同源基因。綜上所述, 本研究獲得了甘藍型油菜響應(yīng)鹽脅迫的高分辨率時間動態(tài)轉(zhuǎn)錄表達譜, 可為后續(xù)研究人員提供可靠的數(shù)據(jù)參考。

研究人員曾在擬南芥中發(fā)現(xiàn)非生物逆境脅迫早期反應(yīng)與后期反應(yīng)的轉(zhuǎn)錄差異, 并認為早期反應(yīng)可能是非特異性的, 主要是在脅迫處理1 h內(nèi), 而基因轉(zhuǎn)錄的特異性應(yīng)激反應(yīng)是在1～3 h后被檢測出來[51]。根據(jù)鹽脅迫樣本層次聚類結(jié)果, 本研究也發(fā)現(xiàn)油菜響應(yīng)鹽脅迫存在明顯的早期響應(yīng)和后期響應(yīng)的聚類差異。對于全程鹽脅迫處理24 h的樣本, 處理1 h前后的根系以及葉片組織均形成了明顯的聚類差異。

3.2 利用WGCNA構(gòu)建鹽脅迫響應(yīng)基因共表達網(wǎng)絡(luò)

WGCNA分析經(jīng)常被用于各種植物的轉(zhuǎn)錄組數(shù)據(jù)分析, 但多數(shù)限制于樣本數(shù)目或處理時間點有限而采用多脅迫或多組織混合分析的方式。本研究利用鹽脅迫處理24 h內(nèi)的多時間點樣本數(shù)據(jù), 可以更好地反映基因響應(yīng)鹽脅迫的表達模式, 從而最大程度地排除其他干擾因素,將不同表達模式的基因歸類到相應(yīng)模塊。同時, 為了避免組織差異性覆蓋脅迫處理的差異性, 本研究分別利用根系以及葉片組織響應(yīng)差異表達基因進行WGCNA分析,構(gòu)建組織響應(yīng)鹽脅迫的基因共表達網(wǎng)絡(luò)。通過分析模塊與樣本時間點的相關(guān)性可以獲得組織在不同處理時間點響應(yīng)鹽脅迫的顯著模塊, 并著重關(guān)注早期響應(yīng)階段的tan模塊和yellow模塊以及后期響應(yīng)階段的green和red模塊。

已有研究表明, 擬南芥中的自噬水平在鹽脅迫處理30 min內(nèi)可達到峰值, 并且這些鹽誘導的PCD事件是由Na+特異引起的[52]。本研究在葉片組織中的tan模塊顯著富集到了自噬反應(yīng)的調(diào)控通路, 該結(jié)果與擬南芥中的報道相對一致。油菜在鹽脅迫處理0.5 h左右, 體內(nèi)自噬調(diào)節(jié)通路就會顯著富集, 可能表明此時油菜幼苗的葉片不單受到鹽處理的滲透作用, 鹽脅迫的離子作用已經(jīng)可以通過根系傳導到葉片, 并引起葉片中鹽脅迫特異的離子脅迫反應(yīng)。植物可以通過多種串擾途徑協(xié)調(diào)各種激素的合成、信號傳導以及新陳代謝來應(yīng)對鹽脅迫, 從而建立有效的防御系統(tǒng)[53]。本研究在早期響應(yīng)模塊中顯著富集到了脫落酸、乙烯、茉莉酸、水楊酸等激素信號通路, 這可能表明油菜在鹽脅迫1 h內(nèi)就會有多種激素信號參與早期脅迫反應(yīng)的調(diào)控。除激活多種激素信號之外, 鹽脅迫早期響應(yīng)階段的根系和葉片組織中也顯著富集到了植物抗毒素合成、植保素代謝、類黃酮合成以及對細菌的防御反應(yīng)等通路。已有研究表明, 植物抗毒素主要包括倍半萜、類黃酮和植保素等物質(zhì), 其中植保素是擬南芥中一種主要的植物抗毒素, 可以由多種植物病原體誘導產(chǎn)生[54]。這可能表明鹽脅迫早期響應(yīng)階段可能和病原菌等生物脅迫具有一定的相似性, 同時也會引起機體合成抗毒素等物質(zhì)進行初級的防御反應(yīng)。

鹽脅迫造成的生理干旱會引起脫落酸積累從而調(diào)控氣孔關(guān)閉。本研究在早期響應(yīng)模塊和后期響應(yīng)模塊均富集到由不同共表達基因參與的脫落酸信號激活通路。這可能表明鹽脅迫全程響應(yīng)過程中, 會有不同的共表達基因在不同的響應(yīng)階段參與調(diào)控相似的生物學通路。與早期響應(yīng)階段相比, 后期響應(yīng)階段也具有特異性。后期響應(yīng)階段的green模塊特異富集到了肌醇合成相關(guān)通路。已有研究表明, 原本是受到鹽脅迫影響的植物, 在肌醇生產(chǎn)保持不間斷時會表現(xiàn)出正常的生長情況, 推測肌醇代謝途徑的選擇性控制可能是提高耐鹽性的方法之一[55]。后期響應(yīng)階段的red模塊特異富集到了植物型次生細胞壁的發(fā)生以及木質(zhì)素代謝過程的通路。已有研究表明, 植物可以通過半纖維素和木質(zhì)素沉積來加強次生細胞壁來增加細胞壁厚度, 從而維持鹽脅迫的穩(wěn)態(tài)[56]。Red模塊也特異富集到了有機酸轉(zhuǎn)運以及苯丙烷合成通路。苯丙烷代謝作為植物重要的次級代謝途徑之一, 其代謝產(chǎn)物如木質(zhì)素以及有機酸等, 在調(diào)控植物適應(yīng)性生長的過程中發(fā)揮著重要功能[57]。以上結(jié)果表明本研究所挑選的關(guān)鍵模塊具有較高的可靠性, 同時也說明利用WGCNA可以構(gòu)建出具有生物學意義的共表達模塊和通路, 這為下一步挖掘核心基因奠定了基礎(chǔ)。

在篩選核心基因時, 本研究發(fā)現(xiàn)BnaA04G0247200ZS由于具有較高的連通度而被篩選出來, 其對應(yīng)擬南芥同源基因WRKY33。已有研究表明, 過表達WRKY33可增強擬南芥的耐鹽性, 該基因可通過控制根細胞凋亡屏障的形成, 賦予根組織的耐鹽能力[58]。油菜中WRKY33可通過增強植保素合成基因的表達增強油菜對菌核病的抗性[59]。本研究分析發(fā)現(xiàn),WRKY33在油菜ZS11中有7個同源拷貝,其中4個拷貝均在鹽脅迫早期響應(yīng)階段的同一個tan模塊中被篩選到, 這說明WGCNA對于挖掘具有多拷貝的核心基因, 也具有較高的適用性和可靠性, 同時推測這些核心轉(zhuǎn)錄因子可能也在油菜響應(yīng)鹽脅迫的過程中發(fā)揮著重要功能。進一步結(jié)合505份鹽脅迫處理的油菜群體變異數(shù)據(jù)庫, 發(fā)現(xiàn)一些核心轉(zhuǎn)錄因子在甘藍型油菜群體中存在較為豐富的SNPs變異以及單倍型類型。其中, 早期響應(yīng)基因BnaA08G0015900ZS在鹽脅迫處理0.5 h時可上調(diào)表達約18倍。其在擬南芥中的同源基因ERF8參與ABA以及免疫信號傳導, 并介導植物對病原體的防御反應(yīng)[60]。單倍型分析發(fā)現(xiàn), 該基因可初步劃分為5種單倍型, 且不同單倍型在鹽脅迫處理前后的根長方面存在極顯著差異。后期響應(yīng)模塊中的核心轉(zhuǎn)錄因子基因BnaA05G0156300ZS(NFYA5)在3 h前與對照組基本無差異, 3 h表達量顯著上升并在24 h達到最大差異倍數(shù)。該基因有3種單倍型, 不同單倍型在幼苗期的植株總干重比(鹽脅迫/對照組)方面存在極顯著差異。如BnaA07G0311300ZS(ZFHD1)也是后期響應(yīng)基因, 有2種單倍型, 不同單倍型在幼苗期的植株總干重比(鹽脅迫/對照組)方面存在極顯著差異。如BnaC01G0510700ZS(NAC047), 有2種單倍型, 不同單倍型在幼苗期的根長比(鹽脅迫/對照組)方面存在極顯著差異。這些基因在模式植物擬南芥中有所報道, 但在油菜中報道還相對較少, 可為后續(xù)甘藍型油菜耐鹽性改良提供重要的候選基因資源。

附表請見網(wǎng)絡(luò)版: 1) 本刊網(wǎng)站http://zwxb.chinacrops.org/; 2) 中國知網(wǎng)http://www.cnki.net/; 3) 萬方數(shù)據(jù)http://c.wanfangdata.com.cn/Periodical-zuowxb.aspx。