基于WGCNA 的馬鈴薯根系抗旱相關(guān)共表達(dá)模塊鑒定和核心基因發(fā)掘

秦天元 孫 超 畢真真 梁文君 李鵬程 張俊蓮白江平,*

1 甘肅省干旱生境作物學(xué)重點(diǎn)實(shí)驗(yàn)室 / 甘肅省作物遺傳改良與栽培種創(chuàng)新重點(diǎn)實(shí)驗(yàn)室, 甘肅蘭州 730070; 2 甘肅農(nóng)業(yè)大學(xué)農(nóng)學(xué)院, 甘肅蘭州 730070

外界脅迫(生物脅迫和非生物脅迫)是農(nóng)業(yè)生產(chǎn)面臨的主要問題之一, 了解脅迫對(duì)植物的影響以及植物對(duì)脅迫信號(hào)的感知、傳遞和應(yīng)答對(duì)解決和提高作物產(chǎn)量具有重要的意義。馬鈴薯作為四大糧食作物之一, 在中國主要種植于年平均降水量小于500 mm 的干旱和半干旱地區(qū), 長(zhǎng)期的或季節(jié)性的干旱脅迫會(huì)嚴(yán)重影響馬鈴薯的植株長(zhǎng)勢(shì)、塊莖產(chǎn)量和商品性[1-2]。馬鈴薯屬淺根系作物, 對(duì)干旱脅迫較為敏感, 因此抗旱性研究一直是馬鈴薯產(chǎn)業(yè)關(guān)注的問題。研究人員從增大根系、減少蒸騰和提高水肥利用率以提高馬鈴薯耐旱性開展了大量研究, 也對(duì)馬鈴薯的根系表型和抗逆相關(guān)生理生化指標(biāo)對(duì)干旱脅迫的響應(yīng)情況進(jìn)行了多種分析, 但針對(duì)馬鈴薯抗旱基因的發(fā)掘和響應(yīng)干旱脅迫的系統(tǒng)性分子機(jī)制研究還相對(duì)較少, 有待進(jìn)一步加強(qiáng)[3-4]。隨著馬鈴薯雙單倍體參考基因組DM 的構(gòu)建完成, 利用轉(zhuǎn)錄組、蛋白組和代謝組等生物信息學(xué)方法研究馬鈴薯的抗逆性已經(jīng)成為一個(gè)有力的手段, 然而前人的研究大多限于單個(gè)調(diào)控因子(如順式表達(dá)元件,反式作用因子), 鮮有利用網(wǎng)絡(luò)分析的方法解析馬鈴薯抗逆機(jī)理的報(bào)道。

近年來, 隨著高通量測(cè)序技術(shù)的飛速發(fā)展與成本不斷降低, 已能采用多樣本的轉(zhuǎn)錄組測(cè)序系統(tǒng)研究生物學(xué)問題。傳統(tǒng)的2 個(gè)樣本的比較分析已經(jīng)不能有效處理海量數(shù)據(jù), 因此網(wǎng)絡(luò)分析應(yīng)運(yùn)而生。其中, 權(quán)重基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene coexpression network analysis, WGCNA)應(yīng)用較為廣泛,已在多種動(dòng)物(如恒河猴[5])和植物(草莓[6]、油菜[7]、玉米[8]等)的多樣本轉(zhuǎn)錄組數(shù)據(jù)的分析中發(fā)揮重要作用。與其他的網(wǎng)絡(luò)分析不同, WGCNA 可以特異地篩選出與目標(biāo)性狀相關(guān)聯(lián)的基因, 并進(jìn)行模塊化分類,得到具有高度生物學(xué)意義的共表達(dá)模塊, 再從中篩選核心基因。因此, 選擇一個(gè)可量化的目標(biāo)性狀是WGCNA 分析首要考慮的問題。本研究選用清除活性氧的抗氧化酶活性等生理生化指標(biāo)作為目標(biāo)性狀。不僅因?yàn)榭寡趸甘窃u(píng)價(jià)植物耐旱性的傳統(tǒng)指標(biāo), 而且在近年來模式植物和模式作物的抗旱性研究中, 活性氧(reactive oxygen species, ROS)可作為第二信使, 促使氧化還原信號(hào)與其他多種信號(hào)通路相互作用, 直接或間接地參與植物的多種代謝過程,共同調(diào)節(jié)植物對(duì)干旱等逆境的響應(yīng)[9-11]。

本試驗(yàn)以5個(gè)時(shí)間梯度的干旱脅迫處理, 及根系轉(zhuǎn)錄組測(cè)序, 結(jié)合抗氧化酶活性和根活力的生理生化數(shù)據(jù), 利用WGCNA 分析構(gòu)建共表達(dá)基因網(wǎng)絡(luò),將基因模塊與生理生化數(shù)據(jù)關(guān)聯(lián)分析, 挖掘出與干旱脅迫調(diào)控通路和抗逆指標(biāo)相關(guān)的核心基因, 以期為進(jìn)一步研究馬鈴薯根系抗旱的分子遺傳機(jī)制提供新的線索。

1 材料與方法

1.1 試驗(yàn)材料

由國際馬鈴薯研究中心引進(jìn)2 個(gè)生育期相同、根系差異明顯的馬鈴薯栽培種C16 (CIP 397077.16)和C119 (CIP 398098.119 (表1), 其組培苗由甘肅農(nóng)業(yè)大學(xué)作物遺傳改良與栽培種創(chuàng)新重點(diǎn)試驗(yàn)室提供。

1.2 試驗(yàn)方法

1.2.1 試驗(yàn)設(shè)計(jì) 試驗(yàn)容器是體積250 mL 的罐頭瓶, 用150 mmol L-1甘露醇模擬干旱環(huán)境。將馬鈴薯試管苗的莖段接種在正常MS 培養(yǎng)基上, 每瓶接生長(zhǎng)狀態(tài)基本一致的5 株莖段, 為了防止污染,每個(gè)處理接5 瓶, 待長(zhǎng)至25 d 時(shí), 取3 瓶長(zhǎng)勢(shì)一致的試管苗轉(zhuǎn)入含有150 mmol L-1甘露醇的液體MS培養(yǎng)基中, 處理0 h、2 h、6 h、12 h 和24 h 后, 對(duì)根系取樣(其中0 h 的樣為對(duì)照), 樣重量不少于2 g,每樣3 個(gè)生物學(xué)重復(fù), 取樣后立即置液氮中速凍,-80℃冰箱保存根, 用于RNA 的提取和后續(xù)轉(zhuǎn)錄組測(cè)序, 以及抗逆相關(guān)生理生化指標(biāo)的測(cè)定。

表1 試驗(yàn)材料Table 1 Experimental materials

1.2.2 轉(zhuǎn)錄組數(shù)據(jù)的獲取與分析 在不同培養(yǎng)條件下, 共采集30 個(gè)樣品[2 個(gè)材料(C16、C119) 5 個(gè)時(shí)間點(diǎn)(0、2、6、12 和24 h) ×3 次重復(fù)]進(jìn)行轉(zhuǎn)錄組測(cè)序。以馬鈴薯基因組及其注釋文件作為參考序列,分析轉(zhuǎn)錄組測(cè)序數(shù)據(jù), 所有數(shù)據(jù)均從Ensembl 網(wǎng)站下 載(http：//plants.ensembl.org/Solanum_tuberosum/Info/Index)。首先使用Hisat2 (v2.2.10)將來自干旱處理和對(duì)照文庫的測(cè)序數(shù)據(jù)和馬鈴薯參考基因組比對(duì),然后利用 Samtools (v0.1.19)將 SAM 文件轉(zhuǎn)換為BAM 文件, 并對(duì)BAM 文件重新排序, 排序后的文件利用cufflinks (v1.2)來組裝轉(zhuǎn)錄本, 估計(jì)這些轉(zhuǎn)錄本的豐度, 并且檢測(cè)樣本間的差異表達(dá)及可變剪切,再利用python 的一個(gè)用于處理高通量數(shù)據(jù)的HTSeq包來計(jì)算樣本的表達(dá)量。根據(jù)基因長(zhǎng)度和匹配到該基因的片段數(shù)目(readcount 值)計(jì)算每一個(gè)基因的FPKM (ragments per kilobase of exon per million fragments mapped), 即每1 百萬個(gè)map 上的reads 中map 到外顯子的每1 千個(gè)堿基上的reads 的個(gè)數(shù)。一般情況下, 當(dāng)基因的FPKM 絕對(duì)值大于1 時(shí)認(rèn)為其表達(dá), 否則認(rèn)為不表達(dá)。將表達(dá)的基因用DESeq2 R package (v3.5.2)軟件和WGCNA 進(jìn)行不同樣本間的差異表達(dá)分析、基于抗逆相關(guān)生理生化數(shù)據(jù)的模塊劃分及核心基因篩選。為了精確定位核心基因, 本試驗(yàn)利用P值小于0.01 的差異表達(dá)基因進(jìn)行后續(xù)篩選。采用 AgriGo (http：//systemsbiology.cau.edu.cn/agriGOv2/)[12-13]在線軟件進(jìn)行GO 富集分析, 利用Cytoscape[14]軟件對(duì)基因調(diào)控網(wǎng)絡(luò)進(jìn)行可視化處理,其中顯著富集的GO 類別P值應(yīng)小于0.05。

2 結(jié)果與分析

2.1 差異基因的分析

通過二代轉(zhuǎn)錄組測(cè)序, 30 個(gè)樣本中共獲得792.64 Gb 和3,801,046,731 個(gè)原始測(cè)序片段, 對(duì)這些原始片段進(jìn)行過濾后得到Q30 大于95.09%的clean reads 3,720,269,696 個(gè), 它們的GC 含量范圍均在41.02%～43.49%之間。30 個(gè)樣本中, clean reads 與基因組的匹配率皆大于85%, 說明所有試驗(yàn)樣品的采集及測(cè)序結(jié)果高度可信。為了篩選不同處理相比于對(duì)照和不同品種間的差異表達(dá)基因, 基于每個(gè)轉(zhuǎn)錄本的readcount 值, 利用DESeq2 軟件進(jìn)行目的樣本間的差異表達(dá)基因分析, 以校正后的P值小于0.01為閾值, 以差異倍數(shù)|log2FC| ＞ 1 為條件, 篩選出差異極顯著的轉(zhuǎn)錄本用于后續(xù)分析。

將C16 和C119 在不同干旱處理時(shí)間點(diǎn)的數(shù)據(jù)分別與未處理的數(shù)據(jù)比對(duì), 篩選出極顯著上調(diào)差異表達(dá)(圖 1-A, B)和極顯著下調(diào)差異表達(dá)的基因(圖1-C, D)。干旱處理2 h 時(shí), C16 中上調(diào)差異表達(dá)基因?yàn)?7 個(gè), 遠(yuǎn)遠(yuǎn)少于C119 中的421 個(gè)(圖1-A,B)。在持續(xù)脅迫到24 h 時(shí), 上調(diào)差異表達(dá)基因數(shù)達(dá)到最大值, C16 中為606 個(gè), 而C119 中為843 個(gè),在處理2 h、6 h、12 h 和24 h 四個(gè)時(shí)間點(diǎn), 相對(duì)于未處理都上調(diào)的基因在C16 中有49 個(gè), 而C119有78 個(gè), 說明在短期干旱處理后C16 根系中上調(diào)表達(dá)基因數(shù)目遠(yuǎn)少于C119。而下調(diào)差異表達(dá)的基因, 脅迫2 h 時(shí) C16 中為12 個(gè), 少于C119 中的50 個(gè)(圖1-C, D)。在持續(xù)脅迫到24 h 時(shí), C16 中下調(diào)表達(dá)的差異基因?yàn)?315 個(gè), 而C119 中為1718個(gè), C16 在持續(xù)處理2 h、6 h、12 h 和24 h 四個(gè)時(shí)間點(diǎn)相對(duì)于0 h 都下調(diào)的基因?yàn)?2 個(gè), 而C119 為21 個(gè), 對(duì)比發(fā)現(xiàn), 在短期處理后C16 根系中上下調(diào)表達(dá)基因數(shù)目遠(yuǎn)少于 C119。以上結(jié)果說明, 相比于C16, C119 在甘露醇模擬的干旱脅迫下有更多的基因積極上下調(diào)表達(dá), 以調(diào)整體內(nèi)代謝途徑的變化來防御逆境傷害。

2.2 差異基因的WGCNA 分析

WGCNA 分析, 即權(quán)重基因共表達(dá)網(wǎng)絡(luò)分析,是一種多個(gè)樣本基因表達(dá)模式的分析方法。通過WGCNA 分析可將表達(dá)模式相似的基因聚類而形成不同的模塊, 分析模塊與特定性狀或表型之間的關(guān)聯(lián)性, 通過基因聚類和關(guān)聯(lián)分析結(jié)果構(gòu)建共表達(dá)調(diào)控網(wǎng)絡(luò)。相比于只關(guān)注差異表達(dá)基因的傳統(tǒng)分析方法, WGCNA 利用數(shù)千或近萬個(gè)變化最大的基因或全部基因分類為不同的基因集, 并與表型進(jìn)行關(guān)聯(lián)分析, 既充分利用了整體組學(xué)信息, 也把數(shù)千個(gè)基因與表型的關(guān)聯(lián)轉(zhuǎn)換為數(shù)個(gè)基因集與表型的關(guān)聯(lián),免去了多重假設(shè)檢驗(yàn)校正。位于調(diào)控網(wǎng)絡(luò)中心的基因被稱為核心基因, 即hub gene, 這類基因通常是關(guān)鍵的調(diào)控基因, 是值得我們深入挖掘和分析的對(duì)象, WGCNA 分析在組學(xué)研究中被廣泛應(yīng)用。本研究將C119 和C16 的根系中響應(yīng)干旱脅迫的差異基因結(jié)合其根系抗逆生理生化指標(biāo)數(shù)據(jù)(附表1), 通過WGCNA 分析和基因集-表型關(guān)聯(lián)分析, 深入挖掘出馬鈴薯根系中抗旱調(diào)控網(wǎng)絡(luò)的核心基因。

圖1 不同處理時(shí)間下C16 與C119 根系中上調(diào)和下調(diào)的差異表達(dá)基因數(shù)Fig. 1 Upregulated and down-regulated genes in roots of C16 and C119 under different treatments

2.2.1 去除離群樣本 由于網(wǎng)絡(luò)模塊分析結(jié)果容易受到離群樣本的影響, 因此在構(gòu)建網(wǎng)絡(luò)模塊之前須去除離群樣本, 以保證結(jié)果的準(zhǔn)確性。本研究共有30 個(gè)樣本, 通過計(jì)算各個(gè)樣本表達(dá)水平的相關(guān)系數(shù)后聚類(圖2-A)。相關(guān)性較低的樣本, 或者在樹形圖上無法聚類的樣本即為離群樣本, 如root.C119.6h.rep1、root.C119.12h.rep3、root.C16. 2h.rep1、root.C16.6h.rep3、root.C16.0h.rep1、root.C16. 12h.rep3。離群樣本被去除后, 剩余樣本的聚類樹如圖 2-B所示。

2.2.2 基因共表達(dá)網(wǎng)絡(luò)中軟閾值的確定 選取C16 和C119 兩個(gè)品種的根系在未處理和短期干旱處理2、6、12 和24 h 得到的基因表達(dá)譜矩陣來構(gòu)建差異基因的共表達(dá)網(wǎng)絡(luò), 需要計(jì)算出每個(gè)基因兩兩之間的相關(guān)系數(shù), 進(jìn)而采用公式Smn= cor(xm,xn),S=[Smn]得到基因之間的相似矩陣, 式中Smn表示基因m和基因n 之間的皮爾森相關(guān)系數(shù),S表示相似矩陣。利用R 軟件中的WGCNA 包來構(gòu)建權(quán)重基因共表達(dá)網(wǎng)絡(luò)。為了使網(wǎng)絡(luò)符合無尺度網(wǎng)絡(luò)分布, 利用WGCNA 包中的函數(shù)pick Soft Threshold 來計(jì)算權(quán)重值。根據(jù)圖3 的結(jié)果, 選擇軟閾值β = 9 來構(gòu)建共表達(dá)網(wǎng)絡(luò)。

2.2.3 基因共表達(dá)網(wǎng)絡(luò)基因聚類和模塊切割 確定軟閾值β=9 之后, 根據(jù)公式為Amn= [(1+Smn)/2]β將相似矩陣轉(zhuǎn)換為鄰接矩陣, 再將鄰接矩陣轉(zhuǎn)為拓?fù)渲丿B矩陣(Topological Overlap Matrix, TOM),同時(shí)利用函數(shù)diss Tom = 1-TOM 對(duì)拓?fù)渲丿B矩陣取逆得到相異性矩陣, 這樣做的目的是為了清除背景噪音和偽關(guān)聯(lián)帶來的誤差。最后利用函數(shù)hclust 對(duì)相異矩陣進(jìn)行層次聚類, 采用動(dòng)態(tài)切割法(Dynamic Tree Cut)將產(chǎn)生的聚類樹切割, 該過程可以把表達(dá)模式相近的基因合并在同一個(gè)分支上,每一個(gè)分支代表一個(gè)共表達(dá)模塊, 不同的顏色代表不同的模塊。差異基因根據(jù)其FPKM 值進(jìn)行相關(guān)度分析并聚類, 相關(guān)度較高的基因被分配到同一個(gè)模塊中(圖4)。

圖2 樣本的聚類樹Fig. 2 Dendrogram of samples

圖3 基因共表達(dá)網(wǎng)絡(luò)軟閾值的確定Fig. 3 Soft threshold determination of gene co-expression network

圖4 基因共表達(dá)網(wǎng)絡(luò)基因聚類數(shù)和模塊切割Fig. 4 Gene co-expression network gene clustering number and modular cutting

圖5 基因共表達(dá)網(wǎng)絡(luò)模塊與生理生化性狀的關(guān)聯(lián)熱圖Fig. 5 Association analysis of gene co-expression network modules with physiological and biochemical traits

2.2.4 基因共表達(dá)網(wǎng)絡(luò)中模塊與生理生化性狀矩陣的關(guān)聯(lián)分析 將上一步所獲模塊中按照75%的相關(guān)性兼并融合, 結(jié)合樣本的4 個(gè)生理生化性狀[超氧化物歧化酶(superoxide dismutase, SOD)、過氧化物酶(peroxidase, POD)、過氧化氫酶(catalase, CAT)和根活力(root vigor, RV)]在不同處理時(shí)間下的矩陣得到15 個(gè)與這些性狀相關(guān)聯(lián)的模塊, 不同的模塊分別以不同的顏色表示(圖5)。其中, 部分模塊與生理生化性狀高度關(guān)聯(lián)。例如, Yellow 模塊內(nèi)的基因表達(dá)模式隨著處理時(shí)間的增加與生理生化性狀的相關(guān)度存在由負(fù)到正的趨勢(shì), 并在24 h 時(shí)與SOD、POD、CAT和RV 均呈極顯著正相關(guān)(分別為r=0.97,P=7E-15;r=0.95,P=2E-12;r=0.87,P=4E-08 和r=0.98,P=3E-16), Red 模塊與與生理生化性狀的關(guān)聯(lián)情況與Yellow 模塊類似; Turquoise 模塊內(nèi)的基因表達(dá)模式則隨著處理時(shí)間的增加與生理生化性狀的相關(guān)度存在由正到負(fù)的趨勢(shì), 并在與處理24 h 時(shí)與SOD、POD、CAT 和根活力(RV)均呈極顯著負(fù)相關(guān)(分別為r=-0.82,P=1E-06;r=-0.81,P=1E-06;r=-0.65,P=6E-04 和r= -0.82,P=1E-06); Blue 模塊與生理生化性狀的關(guān)聯(lián)情況與Turquoise 模塊類似; Pink 模塊與處理2 h 時(shí)的SOD、POD、CAT 和根活力(RV)存在高度正相關(guān)(分別為r=0.80,P=2E-06;r=0.81,P=2E-06;r=0.95,P=3E-12 和r=0.80,P=3E-06)(圖5)。

2.2.5 通過特征向量基因分析確定目標(biāo)基因模塊

基因模塊是指一系列高度相關(guān)基因的集合, 同一個(gè)模塊中的基因可能具有相似的生物學(xué)功能。由于具有相似表達(dá)模式的基因可能具有相似的生物學(xué)功能, 通過WGCNA 構(gòu)建的模塊在很大程度上與特定的性狀相聯(lián)系, 而特定的性狀又與特定的生命活動(dòng)存在高度協(xié)同性。因此, 為了進(jìn)一步解析生理生化性狀特異性模塊的生物學(xué)功能, 根據(jù)基因共表達(dá)網(wǎng)絡(luò)模塊與時(shí)間性狀的關(guān)聯(lián)熱圖, 挑選出與性狀關(guān)聯(lián)度較高的模塊, 并進(jìn)行下一步研究。

由于模塊中的基因數(shù)目龐大, 可利用降維的思路來研究基因模塊, 即把模塊中具有代表性的基因—特征向量基因(module eigengene, ME)挑選出來。用ME 代表模塊中成千上萬的基因進(jìn)行相關(guān)性分析,利用這種方法可以進(jìn)一步理清基因模塊與模塊之間的相互關(guān)系及與不同處理時(shí)間的表型性狀的相關(guān)性, 篩選出目標(biāo)基因模塊。對(duì)所有模塊的ME 進(jìn)行聚類分析后發(fā)現(xiàn), 部分ME 之間存在很高的相關(guān)性(圖6)。ME 之間的相關(guān)性越高, 其所在的模塊的相關(guān)性也就越高, 通過兩兩之間的ME 相關(guān)性分析,發(fā)現(xiàn)Red 模塊與Blue 模塊中ME 之間的相關(guān)性達(dá)到0.87, Yellow 模塊與Turquoise 模塊中ME 之間的相關(guān)性達(dá)到0.88 (圖7), 且Red、Yellow 模塊內(nèi)基因隨著干旱處理時(shí)間的增加, 基因表達(dá)模式與抗逆相關(guān)生理生化性狀由負(fù)相關(guān)逐漸轉(zhuǎn)為極顯著正相關(guān),Tuqurise、Blue 模塊內(nèi)基因隨著處理時(shí)間的增加, 基因表達(dá)模式與目標(biāo)性狀由正相關(guān)逐漸轉(zhuǎn)為極顯著負(fù)相關(guān), 因此這4 個(gè)模塊可作為目標(biāo)基因模塊。

將這4 個(gè)模塊中的所有基因和ME 基因在所有樣本中分別進(jìn)行了表達(dá)水平分析, 發(fā)現(xiàn)每個(gè)模塊內(nèi)的基因表達(dá)水平高度相關(guān), ME 表達(dá)水平與模塊整體表達(dá)水平也高度相關(guān), 說明目標(biāo)模塊的ME 可充分代表其模塊中的整體基因(圖8)。

圖6 ME 聚類樹Fig. 6 ME cluster tree

2.2.6 目標(biāo)基因模塊的GO 注釋 為進(jìn)一步挖掘目標(biāo)基因模塊的功能, 利用AgriGo 在線工具對(duì)篩選出的4 個(gè)目標(biāo)基因模塊Red、Yellow、Turquoise 和Blue 進(jìn)行GO 功能富集分析。結(jié)果表明這4 個(gè)模塊都可以顯著富集到生物學(xué)過程、分子功能以及細(xì)胞組分三大分類下的若干GO 通路(圖9)。如表2 所示,Red 模塊可以富集到氧化還原相關(guān)基因(GO：0016705)和內(nèi)肽酶抑制劑活性基因(GO：0004866)等調(diào)控通路; Yellow 模塊可以富集到與蛋白質(zhì)結(jié)合相關(guān)的通路(GO：0005515); Turquoise 模塊可以富集到光合作用(GO：0015979), 也富集到與氧化還原酶活性相關(guān)的通路(GO：0016705), 還富集到類固醇合成途徑中催化類固醇脫氫反應(yīng)(以類固醇羥基為氫供體, 以NAD 或NADP 為氫受體)的脫氫酶及相關(guān)通路(GO：0006694、GO：0003854、GO： 0033764、GO：0016229、GO： 0016614、GO： 0016616 和GO：0051287);Blue 模塊可以富集到與氧化還原酶活性相關(guān)的基因(GO：0016705) 和 與 細(xì) 胞 骨 架 相 關(guān) 的 基 因(GO：0007010 和GO：0015629)。總體來看, 4 個(gè)模塊都富集到了氧化還原酶活性及細(xì)胞骨架組織相關(guān)的通路, 說明WGCNA 可以有效構(gòu)建到具有生物學(xué)意義的共表達(dá)模塊, 這些模塊可成為本研究的重點(diǎn)關(guān)注對(duì)象。

2.2.7 目標(biāo)模塊中核心基因的篩選和功能分析

圖7 不同模塊兩兩之間ME 的相關(guān)性Fig. 7 ME correlation between different modules

為了獲得上述 4 個(gè)模塊中的核心基因, 利用Cytoscape 軟件對(duì)基因調(diào)控網(wǎng)絡(luò)進(jìn)行可視化處理, 篩選出模塊中高連通性的基因, 這些基因則作為該模塊中的核心基因(圖10 和圖11)。在這些網(wǎng)絡(luò)中, 每一個(gè)節(jié)點(diǎn)代表一個(gè)基因, 節(jié)點(diǎn)通過線相聯(lián)系, 處于連線兩端的基因通常被認(rèn)為具有相同的生物功能。其中在Red 模塊中篩選出了5 個(gè)核心基因, Yellow模塊中篩選出了9 個(gè)核心基因, Turquoise 模塊中篩選出4 個(gè)核心基因, Blue 模塊中篩選出了2 個(gè)核心基因,為了獲取這些核心基因的功能信息, 我們利用馬鈴薯基因數(shù)據(jù)庫(http：//solanaceae.plantbiology.msu.edu/index.shtml)和 NCBI 數(shù)據(jù)庫(https：//www.ncbi.nlm.nih.gov/)查詢這些核心基因的相關(guān)報(bào)道, 并借助TAIR 數(shù)據(jù)庫(https：//www.arabidopsis.org/)注釋這些核心基因在擬南芥中的同源基因的功能(表3), 通過已報(bào)道的文獻(xiàn), 本研究發(fā)現(xiàn)這些篩選出的核心基因基本都與干旱脅迫相關(guān), 其中PGSC0003 DMG400026572在擬南芥中的同源基因AT3G51730基因參與了植物液泡的形成, 而液泡在植物生長(zhǎng)、發(fā)育和應(yīng)激反應(yīng)中起著中心作用。PGSC0003DMG 400029350在擬南芥中的同源基因AT1G05680基因通過調(diào)節(jié)活性氧和氧化還原信號(hào)與植物激素發(fā)生協(xié)同或拮抗作用, 引起植物對(duì)生物和非生物脅迫的保護(hù)反應(yīng)。PGSC0003DMG400006704在擬南芥中的同源基因AT5G35200基因通過誘導(dǎo)啟動(dòng)蛋白磷酸化調(diào)節(jié)下游的信號(hào)包括穿過質(zhì)膜的離子通量、MAPK 的活化和活性氧的形成。PGSC0003 DMG400024687在擬南芥中的同源基因AT1G55000基因與泛素化的特異性相關(guān), 能促進(jìn)泛素轉(zhuǎn)移到合適的靶點(diǎn)進(jìn)而對(duì)細(xì)胞活動(dòng)進(jìn)行調(diào)控。PGSC0003 DMG400003792在擬南芥中的同源基因AT3G14110基因與植物體內(nèi)的活性氧的變化相關(guān), 在逆境環(huán)境下, 該基因可激活多種應(yīng)激反應(yīng)。

圖8 各樣本中不同模塊的所有基因與相應(yīng)ME 的表達(dá)水平Fig. 8 Expression levels of all genes and corresponding ME in different modules of each sample

圖9 目標(biāo)模塊內(nèi)基因GO 注釋Fig. 9 Gene GO annotation in target module

表2 目標(biāo)模塊的部分GO 富集分析結(jié)果Table 2 Partial GO enrichment analysis of target modules

(續(xù)表2)

圖10 Red 和Yellow 模塊內(nèi)的基因共表達(dá)網(wǎng)絡(luò)及其核心基因Fig. 10 Gene co-expression network and hub genes in Red and Yellow modules

圖11 Turquoise 和Blue 模塊內(nèi)的基因共表達(dá)網(wǎng)絡(luò)及其核心基因Fig. 11 Gene co-expression network and hub genes in Turquoise and Blue modules

表3 不同模塊中核心基因的功能注釋Table 3 Functional annotations of hub genes in different modules

2.2.8 核心基因的RT-qPCR 驗(yàn)證 從相關(guān)性最高的4 個(gè)基因模塊中篩選到20 個(gè)可能與干旱脅迫相關(guān)的核心基因(表3), 從中挑選了14 個(gè)基因進(jìn)行實(shí)時(shí)熒光定量PCR 驗(yàn)證。將用于轉(zhuǎn)錄組測(cè)序的備份RNA 反轉(zhuǎn)錄, 通過q-PCR 檢測(cè)這些基因在C16 和C119 中對(duì)干旱脅迫的響應(yīng)情況。最終結(jié)果的變化趨勢(shì)與轉(zhuǎn)錄組結(jié)果基本一致(表4 和圖12), 進(jìn)一步證明了轉(zhuǎn)錄組測(cè)序結(jié)果的可靠性, 并驗(yàn)證了這些核心基因確實(shí)響應(yīng)干旱脅迫。

(續(xù)表3)

表4 字母與基因的對(duì)應(yīng)表Table 4 Correspondence between letters and genes

3 討論

干旱脅迫是影響自然界植物地理分布, 限制農(nóng)業(yè)生產(chǎn)力和威脅糧食安全的主要環(huán)境因素之一。全球氣候變化更加劇了干旱脅迫的不利影響, 預(yù)計(jì)還將導(dǎo)致極端天氣的頻率增加[20-21]。馬鈴薯作為第四大糧食作物, 在我國的主產(chǎn)區(qū)也處在氣候變化的高風(fēng)險(xiǎn)區(qū)域, 干旱脅迫是主要的風(fēng)險(xiǎn)之一[22]。因此, 提高馬鈴薯的抗旱性, 解析其在干旱脅迫下的感知和適應(yīng)機(jī)制對(duì)提高農(nóng)業(yè)生產(chǎn)力和節(jié)約環(huán)境資源均具有重要意義。

植物根系是植物汲取和轉(zhuǎn)化水分和營養(yǎng)的主要器官, 也是植物面臨土壤脅迫環(huán)境時(shí)首先感知脅迫信號(hào)的器官, 在植物應(yīng)對(duì)非生物逆境尤其是干旱脅迫時(shí)起關(guān)鍵作用[23-25]。雖然馬鈴薯根系抗逆遺傳機(jī)理的解析取得了一定的進(jìn)展, 例如通過數(shù)量性狀位點(diǎn)定位和關(guān)聯(lián)分析等方法鑒定了許多重要的根系抗逆基因[26-27], 但是這些研究手段都有一定的局限性,例如數(shù)量性狀定位不能反映遺傳背景復(fù)雜群體的遺傳異質(zhì)性, 此外, 大多數(shù)的農(nóng)藝性狀還與整個(gè)群體結(jié)構(gòu)存在很高的相關(guān)性, 而利用現(xiàn)在生物信息技術(shù)和數(shù)據(jù)分析方法則可以通過共表達(dá)模塊化處理, 將成千上萬具有相似生物學(xué)功能的基因分配到同一個(gè)模塊中, 通過研究模塊的生物學(xué)意義來進(jìn)一步挖掘模塊內(nèi)基因的功能, 較好地解決了對(duì)遺傳背景復(fù)雜群體的性狀解析不足的問題, 因此被廣泛地應(yīng)用于大數(shù)據(jù)組學(xué)的研究之中[28-30]。因此, 本研究干旱處理馬鈴薯幼苗后, 收集其根系材料進(jìn)行轉(zhuǎn)錄組分析,以期獲得馬鈴薯根系抗逆相關(guān)的新的研究線索。與傳統(tǒng)的僅對(duì)比處理前后的轉(zhuǎn)錄組差異分析不同, 本研究針對(duì)2 個(gè)品種, 分別設(shè)置了4 個(gè)處理時(shí)間梯度,加上未處理的對(duì)照, 每樣3 個(gè)生物學(xué)重復(fù), 共計(jì)30個(gè)樣本進(jìn)行了轉(zhuǎn)錄組測(cè)序。與此同時(shí), 本研究還檢測(cè)了所有處理時(shí)間下每個(gè)樣本的脅迫相關(guān)生理生化指標(biāo)(SOD、POD、CAT 和RV)。

圖12 核心基因的RT-qPCR 驗(yàn)證Fig. 12 RT-qPCR validation of Hub genes

大量的轉(zhuǎn)錄組數(shù)據(jù)和生理數(shù)據(jù)獲得后, 如何從中挖掘蘊(yùn)含的生物學(xué)意義成為首要考慮的問題。網(wǎng)絡(luò)分析在基因組、轉(zhuǎn)錄組和代謝組等高通量組學(xué)數(shù)據(jù)的發(fā)掘中應(yīng)用較為廣泛。為了發(fā)掘與抗旱相關(guān)的基因或調(diào)控通路, 本研究采用了WGCNA 分析方法,可特異地篩選出與目標(biāo)性狀相關(guān)的基因, 并進(jìn)行模塊化分類, 得到具有高度生物學(xué)意義的共表達(dá)模塊,已經(jīng)被證明是一種高效的數(shù)據(jù)挖掘手段[31]。ROS 的產(chǎn)生是植物受到干旱脅迫后最早的反應(yīng)之一, ROS是許多脅迫相關(guān)生物合成、信號(hào)轉(zhuǎn)導(dǎo)和代謝相關(guān)途徑的重要調(diào)節(jié)劑, 其過量積累對(duì)細(xì)胞具有毒害作用[32], 因此維持其動(dòng)態(tài)平衡是植物應(yīng)對(duì)脅迫的重要調(diào)節(jié)方式, 而SOD、POD 和CAT 作為清除ROS的抗氧化酶, 可較為直接地反映植物對(duì)脅迫的響應(yīng)情況[33], 可作為WGCNA 分析中的目標(biāo)性狀。因此,我們利用WGCNA 分析方法, 富集到了15 個(gè)與生理性狀相關(guān)聯(lián)的共表達(dá)模塊, 其中4 個(gè)模塊與這些抗逆性狀高度相關(guān), 通過進(jìn)一步的網(wǎng)絡(luò)分析, 找到了這4 個(gè)模塊的核心基因。表明WGCNA 可以構(gòu)建到具有生物學(xué)意義的共表達(dá)模塊和通路, 為下一步深入研究打下基礎(chǔ)。

本研究通過WGCNA 富集到的15 個(gè)與脅迫相關(guān)生理生化指標(biāo)(SOD、POD、CAT 和RV)相關(guān)聯(lián)的共表達(dá)模塊, 由于模塊較多, 因此挑選了4 個(gè)與生理性狀關(guān)聯(lián)度最高的模塊進(jìn)行下一步研究。模塊內(nèi)的基因表達(dá)模式隨著處理時(shí)間的增加, 存在與生理性狀由負(fù)相關(guān)逐漸轉(zhuǎn)為極顯著正相關(guān)的Yellow 模塊和Red 模塊, 以及存在相反趨勢(shì)的Turquoise 模塊和Blue 模塊。由于每個(gè)模塊的基因數(shù)量龐大, 本研究利用各個(gè)模塊的特征向量基因?qū)δK兩兩之間進(jìn)行相關(guān)性分析, 結(jié)果表明, 這4 個(gè)模塊間的確存在極顯著的相關(guān)性且關(guān)聯(lián)度最高, 進(jìn)一步印證這4 個(gè)模塊可作為根系抗逆性狀機(jī)理研究的目標(biāo)模塊。通過對(duì)這4 個(gè)模塊的GO 功能富集分析發(fā)現(xiàn), 均可以得到與抗逆相關(guān)的具有生物學(xué)意義的調(diào)控途徑。例如,Red 模塊、Turquoise 模塊和Blue 模塊都富集到與氧化還原活性相關(guān)的通路(GO：0016705), Turquoise 模塊還富集到類固醇合成途徑中催化類固醇脫氫反應(yīng)的脫氫酶及相關(guān)通路(GO：0006694)。有研究表明, 在擬南芥和玉米中過量表達(dá)油菜素固醇合成途徑的脫氫酶基因, 可導(dǎo)致ABA 的積累, 增強(qiáng)植物抗性, 提高玉米在逆境下的產(chǎn)量[34-35]。說明WGCNA 分析的確可富集到與目標(biāo)性狀相關(guān)的具有生物學(xué)意義的共表達(dá)模塊。

本研究針對(duì)4 個(gè)與目標(biāo)性狀關(guān)聯(lián)度最高的模塊分析表明, 其核心基因包含干旱脅迫相關(guān)基因, 例如Red 模塊中的參與木質(zhì)部分化的基因(PGSC0003 DMG400020562)、Yellow 模塊中的 MAP 激酶(PGSC0003DMG400020683)、Turquoise 模塊可能參與葉綠素合成的基因(PGSC0003DMG400003792)。其中Yellow 模塊中PGSC0003DMG400029350基因在擬南芥中的同源基因AT1G05680可編碼一種UDP-葡萄糖基轉(zhuǎn)移酶UGT74E2, 在擬南芥中的研究表明, 其可被過氧化氫或非生物脅迫誘導(dǎo)表達(dá),通過對(duì)生長(zhǎng)素吲哚-3-丁酸(IBA)和吲哚-3-乙酸(IAA)的動(dòng)態(tài)平衡調(diào)節(jié), 參與氧化還原信號(hào)與植物激素信號(hào)的整合互作, 增加植物在非生物逆境中的保護(hù)反應(yīng)[36]。這些基因在模式植物擬南芥和其他作物如玉米、水稻等中有較多的研究, 在馬鈴薯中的報(bào)道還相對(duì)較少, 可進(jìn)行深入探究。

以上結(jié)果表明, 利用WGCNA 分析, 可挖掘到與目標(biāo)性狀的生物學(xué)意義高度關(guān)聯(lián)的基因模塊和核心基因, 為挖掘目標(biāo)基因提供新的研究思路, 為解析復(fù)雜的農(nóng)藝性狀提供參考。在本研究中, 重點(diǎn)關(guān)注了Red、Yellow、Turquoise 和Blue 四個(gè)相關(guān)性高的調(diào)控模塊, 其余的基因模塊雖然沒有被詳細(xì)論述,但也可能包含與根系抗逆相關(guān)的通路, 可進(jìn)一步挖掘其蘊(yùn)含的生物學(xué)意義。

4 結(jié)論

本研究構(gòu)建了與抗逆生理性狀相關(guān)聯(lián)的權(quán)重基因共表達(dá)網(wǎng)絡(luò), 得到了15 個(gè)與根系抗旱密切相關(guān)的基因模塊并進(jìn)行了GO 通路富集。挑選出4 個(gè)與目標(biāo)性狀關(guān)聯(lián)度最佳的基因模塊進(jìn)行深入分析, 其中在Red 模塊中挑選出了5 個(gè)核心基因, Yellow 模塊中挑選出了9 個(gè)核心基因, Turquoise 模塊中挑選出4個(gè)核心基因, Blue 模塊中挑選出了2 個(gè)核心基因, 通過功能注釋, 發(fā)現(xiàn)部分核心基因與已報(bào)道非生物脅迫調(diào)控通路緊密相關(guān), 還有一部分功能僅在模式植物中有一些研究。本研究結(jié)果可為馬鈴薯根系抗旱的分子機(jī)制研究提供線索, 也為耐旱性馬鈴薯新品種培育提供理論支持。

附表 1 C16 和C119 在不同處理時(shí)間下的生化指標(biāo)數(shù)據(jù)Supplementary table 1 Biochemical indicator data for C16 and C119 at different treatment times