傅明川 李 浩 陳義珍 柳展基 劉任重 王立國

利用WGCNA鑒定棉花抗黃萎病相關(guān)基因共表達(dá)網(wǎng)絡(luò)

傅明川*李 浩 陳義珍 柳展基 劉任重 王立國

山東棉花研究中心 / 農(nóng)業(yè)部黃淮海棉花遺傳改良與栽培生理重點實驗室, 山東濟南 250100

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(Weighted Gene Co-expression Network Analysis, WGCNA)作為一種典型的系統(tǒng)生物學(xué)方法, 可用來鑒定協(xié)同表達(dá)的基因模塊, 探索模塊與目標(biāo)性狀之間的生物學(xué)相關(guān)性, 并挖掘網(wǎng)絡(luò)中的核心基因。黃萎病(wilt)可嚴(yán)重降低棉花(spp.)的產(chǎn)量及品質(zhì), 因此研究抗病基因及抗病機制, 對于棉花抗病育種工作具有重要意義。本研究以黃萎病菌侵染不同時間點的海島棉()幼苗根尖轉(zhuǎn)錄組數(shù)據(jù)為基礎(chǔ), 分析其差異表達(dá)基因, 并選取在不同樣本之間表達(dá)水平變異最大的前50%基因(共35,647個)進(jìn)行WGCNA。結(jié)果表明, 在黃萎病菌侵染條件下, 共有22,850個基因差異表達(dá), 其中4685個為所有侵染時間點共有的差異基因。共鑒定到18個基因共表達(dá)模塊, 其中5個為抗黃萎病相關(guān)特異性模塊(black、mediumpurple3、darkolivegreen、plum3模塊分別與侵染2 h、6 h、48 h、72 h時間點正相關(guān); mediumpurple2與侵染2 h時間點負(fù)相關(guān))。GO和KEGG富集分析表明, 特異性模塊可以富集到有生物學(xué)意義的富集結(jié)果, 如刺激響應(yīng)、鈣離子結(jié)合、黃酮類化合物合成代謝等抗逆相關(guān)通路。通過計算模塊內(nèi)基因的連通性, 挖掘出了網(wǎng)絡(luò)中的核心基因, 功能預(yù)測表明, 這些基因可能在逆境脅迫響應(yīng)中發(fā)揮重要作用。本研究為進(jìn)一步理解棉花抗黃萎病的分子機制、選育棉花抗病新品種提供了理論指導(dǎo)。

棉花; 黃萎病; 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析; 核心基因

棉花(spp.)是世界范圍內(nèi)重要的纖維及油料作物之一, 其在我國經(jīng)濟發(fā)展中也占有十分重要的地位。黃萎病(wilt)作為土壤傳播的維管束系統(tǒng)性病害, 是棉花生產(chǎn)中最嚴(yán)重的病害之一。過去通常認(rèn)為黃萎病菌的致病機制是“導(dǎo)管堵塞”; 目前研究表明, 黃萎病菌分泌的毒素可能是其主要致病因素, 其中細(xì)胞壁降解酶、脂多糖蛋白復(fù)合物以及核蛋白可能是重要的致病成分。Wang等[1]研究證實, 黃萎病菌VdNEP蛋白是一種重要的黃萎病致病因子。通過研究高致病性菌株Vd080發(fā)現(xiàn),基因與黃萎病菌的致病力和微菌核形成有關(guān)[2]。構(gòu)成了黃萎病菌的RGS基因家族, 其中在孢子生成、菌絲生長、微菌核形成過程中發(fā)揮重要作用[3]。由于缺乏有效的黃萎病防治藥劑, 選育和種植抗病品種是目前最經(jīng)濟有效、綠色環(huán)保的防治措施。但抗病基因資源缺乏、抗病機制不明, 導(dǎo)致棉花抗黃萎病育種工作進(jìn)展緩慢[4]。因此, 挖掘棉花抗病基因、探究其抗病機制, 對于棉花抗病育種工作具有重要意義。Kawchuk等[5]利用圖位克隆在番茄()位點中分離出2個緊密相連的抗黃萎病主效基因和,研究證實基因可介導(dǎo)感病馬鈴薯品種()對黃萎病菌的抗性。進(jìn)一步研究表明, 由于棉花黃萎病菌缺少基因, 導(dǎo)致在棉花中過表達(dá)并不能明顯增強其抗病性[6]。基因組學(xué)、轉(zhuǎn)錄組學(xué)及蛋白質(zhì)組學(xué)的迅速發(fā)展, 為研究棉花抗黃萎病的分子機制提供了有利基礎(chǔ)。Zhang等[7]通過對黃萎病菌侵染下的海島棉()“Pima 90-53”幼苗根系進(jìn)行轉(zhuǎn)錄組分析, 首次報道了“SA→NPR1→TGA→PR-1→抗病”代謝通路。目前, 在棉花中已鑒定出一些抗病相關(guān)基因。如是海島棉中一個ERF轉(zhuǎn)錄因子基因, 轉(zhuǎn)的煙草()植株內(nèi)抗病相關(guān)基因表達(dá)量升高, 并可提高轉(zhuǎn)基因植株的抗病性[8];為一類絲氨酸/蘇氨酸蛋白激酶基因, 參與調(diào)控多個抗逆相關(guān)信號通路, 從而在抵御病原菌侵染脅迫及氧化應(yīng)激響應(yīng)中發(fā)揮作用[9];參與調(diào)控活性氧含量及JA信號通路相關(guān)基因的表達(dá), 轉(zhuǎn)基因的擬南芥()植株對黃萎病菌的抗性提高[10]。此外,[11]、[12]、[13]等基因也均在抗病響應(yīng)中發(fā)揮重要作用。

目前, 棉花基因組測序工作已陸續(xù)完成, 這為從海量生物數(shù)據(jù)中挖掘有用信息提供了堅實基礎(chǔ)。相比傳統(tǒng)分子生物學(xué)方法, 借助生物信息學(xué)手段, 可更加快捷地定位目標(biāo)基因, 挖掘與性狀相關(guān)的關(guān)鍵基因。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(Weighted Gene Co-expression Network Analysis, WGCNA)以芯片或RNA-seq表達(dá)數(shù)據(jù)為基礎(chǔ), 通過冪指數(shù)加權(quán)構(gòu)建無尺度拓?fù)渲丿B矩陣, 用以描述基因之間的相互關(guān)系, 并依此將表達(dá)模式相近的基因劃分到一個基因表達(dá)模塊中[14]。WGCNA多用來研究協(xié)同表達(dá)的基因模塊與目標(biāo)性狀之間的生物學(xué)相關(guān)性, 并探索基因共表達(dá)網(wǎng)絡(luò)中的核心基因。作為一種典型的系統(tǒng)生物學(xué)方法, WGCNA已廣泛應(yīng)用于植物學(xué)研究中。如Tan等[15]通過分析鎘處理不同時間點的17個水稻 ()轉(zhuǎn)錄組數(shù)據(jù), 鑒定得到22個基因模塊, 結(jié)合差異表達(dá)分析, 共挖掘到164個鎘脅迫響應(yīng)相關(guān)基因; Zou等[16]通過對2個棉花品系不同發(fā)育時期的纖維轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行WGCNA, 共鑒定得到5個纖維發(fā)育相關(guān)特異性模塊, 并挖掘出模塊內(nèi)的核心基因; 通過對14份玉米()不同發(fā)育階段的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行WGCNA, 研究者鑒定得到14個組織特異性模塊, 并對其中2個模塊的基因互作網(wǎng)絡(luò)進(jìn)行了進(jìn)一步研究, 從中挖掘到了、、等開花相關(guān)核心基因[17]。

本研究以黃萎病菌侵染不同時間點的棉花轉(zhuǎn)錄組數(shù)據(jù)為材料, 對其進(jìn)行差異表達(dá)分析; 通過構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò), 劃分基因模塊, 并篩選出抗病相關(guān)特異性模塊; 經(jīng)GO及KEGG富集分析, 探究模塊功能; 根據(jù)基因在相應(yīng)網(wǎng)絡(luò)中的連通性, 鑒定出模塊內(nèi)的核心基因。本研究可為進(jìn)一步理解棉花抗黃萎病的分子機制提供理論基礎(chǔ), 并為棉花抗病育種工作提供新的基因資源。

1 材料與方法

1.1 數(shù)據(jù)獲取及差異表達(dá)分析

棉花在黃萎病菌侵染不同時間點的轉(zhuǎn)錄組數(shù)據(jù)下載自NCBI數(shù)據(jù)庫(PRJNA234454)[18]。該數(shù)據(jù)以高抗黃萎病海島棉品種“海7124”為材料, 對生長2周的幼苗根系使用Vd8孢子懸液侵染。分別選取侵染2、6、12、24、48和72 h的根尖, 提取總RNA, 并以未侵染材料(0 h)為對照。首先, 利用fastq-dump軟件(https://ncbi.github.io/sra-tools/fastq-dump.html)將下載得到的SRA文件轉(zhuǎn)換為fastq格式文件, 然后利用FastQC (http://www.bioinformatics.babraham. ac.uk/projects/fastqc/)對測序結(jié)果進(jìn)行質(zhì)量評估, 并通過Trimmomatic軟件[19]進(jìn)行質(zhì)控, 將得到的clean data用于后續(xù)分析。以海島棉基因組為參考基因組[20], 使用HISAT2[21]進(jìn)行序列比對后, 利用featureCounts[22]計數(shù), 得到每個基因在各個樣本中的raw counts。將數(shù)據(jù)導(dǎo)入R中, 利用DESeq2[23]進(jìn)行差異表達(dá)分析。本研究選取|log2FC| > 1 (log2FC代表處理與對照表達(dá)量比值的對數(shù)值)及adj<0.001 (adj代表校正后的值)的基因作為差異表達(dá)基因。

1.2 加權(quán)基因共表達(dá)網(wǎng)絡(luò)構(gòu)建

利用R程序中的WGCNA軟件包[14]進(jìn)行加權(quán)基因共表達(dá)網(wǎng)絡(luò)構(gòu)建。以標(biāo)準(zhǔn)化后的基因表達(dá)矩陣做為輸入, 共21個轉(zhuǎn)錄組樣本(7個時間點, 各3次重復(fù))。通過計算每個基因在各個樣本間表達(dá)水平的變異程度, 選取變異最大的前50%基因進(jìn)行WGCNA。經(jīng)過閾值篩選, 最終選擇=18對原始有尺度關(guān)系矩陣進(jìn)行冪處理, 得到無尺度化鄰接矩陣。為更好評估基因間表達(dá)模式的相關(guān)性, 進(jìn)一步將鄰接矩陣轉(zhuǎn)化為拓?fù)渲丿B矩陣(Topological Overlap Matrix, TOM), 并利用拓?fù)湎喈惥仃?dissTOM=1-TOM), 采用動態(tài)剪切算法進(jìn)行基因聚類及模塊劃分。模塊內(nèi)最少基因數(shù)為30 (minModuleSize=30), 相似模塊合并閾值為0.25 (cutHeight=0.25), 網(wǎng)絡(luò)類型為“signed” (type="signed"或networkType="signed", 依不同函數(shù)而定)。

1.3 特異性模塊篩選

對每個模塊中的所有基因進(jìn)行主成分分析(Principle Component Analysis, PCA), 將主成分1(PC1)的值稱為該模塊的模塊特征向量(Module Eigengene, ME)。為篩選抗病相關(guān)特異性模塊, 分別計算每個模塊的模塊特征向量與不同性狀之間(此處為不同侵染時間)的相關(guān)系數(shù)及相應(yīng)值。>0代表正相關(guān),<0代表負(fù)相關(guān)。本研究選擇||>0.70及<0.001的模塊作為特異性模塊進(jìn)一步分析。

1.4 富集分析及代謝通路分析

利用R程序中的clusterProfiler軟件包[24]進(jìn)行GO (Gene Ontology)和KEGG (Kyoto Encyclopedia of Genes and Genomes)富集分析。閾值為<0.01及<0.05。

1.5 轉(zhuǎn)錄因子分析

將各模塊中的蛋白序列提交到PlantTFDB數(shù)據(jù)庫(http://planttfdb.cbi.pku.edu.cn/prediction.php#)[25]分析預(yù)測, 從而得到每個模塊中的轉(zhuǎn)錄因子。

1.6 RNA提取及qRT-PCR

選取生長2周、長勢一致的“海7124”幼苗, 采取浸根法用濃度1×107個mL-1的黃萎病菌Vd8孢子懸液侵染10 min。選取侵染后2、6、12、24、48、72 h和相應(yīng)時間點未侵染幼苗根系, 利用TRIzol法提取總RNA, 取1 μg RNA反轉(zhuǎn)錄為cDNA。

qPCR使用Thermo Fisher Scientific的QuantStudio 5實時熒光定量PCR系統(tǒng), 選用Aidlab公司的2×Sybr Green qPCR Mix試劑盒, 熒光染料為SYBR Green, 內(nèi)參基因為。反應(yīng)程序95℃預(yù)變性3 min; 95℃ 15 s, 60℃ 15 s, 40個循環(huán); 熔點曲線程序為95℃ 15 s, 60℃ 1 min, 95℃ 15 s。反應(yīng)體系包含2×SYBR qPCR Mix 10 μL、DNA Template (稀釋10倍) 0.8 μL、正向引物(10 μmol L-1) 0.4 μL、反向引物(10 μmol L-1) 0.4 μL、ddH2O 8.4 μL。使用2–DDCt法分析基因相對表達(dá)量[26], 設(shè)置3次生物學(xué)重復(fù)。所用基因引物見表1。

2 結(jié)果與分析

2.1 差異表達(dá)分析

通過分析棉花在黃萎病菌侵染不同時間點(0 h、2 h、6 h、12 h、24 h、48 h、72 h)的轉(zhuǎn)錄組數(shù)據(jù), 最終共得到22,850個差異表達(dá)基因(|log2FC| > 1,adj< 0.001)。其中, 顯著上調(diào)基因9398個, 顯著下調(diào)基因13,171個, 另有281個基因在不同時間點表現(xiàn)出不同的上下調(diào)趨勢。差異表達(dá)基因最多的時間點出現(xiàn)在侵染后24 h (14,679個), 最少的為2 h (10,729個), 所有時間點共有的差異表達(dá)基因為4685個(圖1)。

2.2 加權(quán)基因共表達(dá)網(wǎng)絡(luò)構(gòu)建

對表達(dá)矩陣中變異較低的基因進(jìn)行過濾, 最終選取35,647個基因進(jìn)行WGCNA。當(dāng)=18時, 無尺度網(wǎng)絡(luò)擬合指數(shù)2>0.8, 平均連通性趨近于0, 表明用此值進(jìn)行冪處理可以得到符合要求的無尺度網(wǎng)絡(luò), 因此選擇=18構(gòu)建無尺度網(wǎng)絡(luò)(圖2)。采用動態(tài)剪切算法對基因進(jìn)行聚類及模塊劃分, 通過計算每個模塊的模塊特征向量, 合并相似模塊, 最終共得到18個基因共表達(dá)模塊(圖3)。模塊內(nèi)基因數(shù)量為61 (palevioletred2) ~ 8431 (black)個。各模塊內(nèi)差異表達(dá)基因所占比例為11% (thistle2) ~ 81% (paleturquoise)。值得注意的是, 除black及mediumpurple2模塊外, 同一模塊內(nèi)的差異表達(dá)基因通常表現(xiàn)出相似的上下調(diào)模式, 如paleturquoise、turquoise模塊內(nèi)的差異基因均為顯著下調(diào), 而palevioletred2、thistle2模塊內(nèi)均為顯著上調(diào)。

圖1 不同侵染時間點的差異表達(dá)基因

只有一個黑點的列代表某個數(shù)據(jù)集特有的差異基因; 有2個或以上用實線連接的黑點的列代表相應(yīng)數(shù)據(jù)集之間共有的差異基因。

The column with one black point represents the DEGs only in the corresponding set; the column with two or more black points linked by a solid line represents the intersection among the sets.

轉(zhuǎn)錄因子是生物過程中一類重要的調(diào)控蛋白, 本研究也對每個模塊內(nèi)的轉(zhuǎn)錄因子進(jìn)行了進(jìn)一步分析。結(jié)果表明, 模塊內(nèi)轉(zhuǎn)錄因子數(shù)量分布為3 (mediumpurple2) ~ 590 (black)個, 在相應(yīng)模塊內(nèi)所占比例為2% (mediumpurple2) ~ 19% (plum3), 不同模塊內(nèi)轉(zhuǎn)錄因子的分布類型存在差異, 但主要集中在ERF (9.27%)、MYB (8.09%)、bHLH (7.94%)、WRKY (7.03%)、C2H2 (5.74%)、NAC (5.36%)、bZIP (5.28%)等轉(zhuǎn)錄因子家族, 報道表明這些轉(zhuǎn)錄因子均參與調(diào)控植物的抗逆過程[27]。每個模塊內(nèi)差異表達(dá)的轉(zhuǎn)錄因子數(shù)占該模塊內(nèi)所有轉(zhuǎn)錄因子的比例為20% (thistle2) ~ 89% (navajowhite1), 數(shù)量顯著偏高(χ2=92.24,=2e-12), 進(jìn)一步表明轉(zhuǎn)錄因子在抗逆調(diào)控過程中可能發(fā)揮重要作用。

圖2 軟閾值的選擇

a: 不同軟閾值下的無尺度網(wǎng)絡(luò)擬合指數(shù)(2), 紅線代表2=0.8。b: 不同軟閾值下的平均連通性。

a: scale-free topology fit index as a function of the soft-thresholding power, the red line represents that2is equal to 0.8. b: mean connectivity as a function of the soft-thresholding power.

圖3 基因聚類樹及模塊劃分

a: 基于拓?fù)湎喈惥仃嚇?gòu)建的基因聚類樹。b: 使用動態(tài)剪切算法得到的基因模塊, 不同顏色代表不同模塊。c: 合并相似模塊后的模塊劃分結(jié)果。

a: gene clustering on TOM-based dissimilarity. b: module division by dynamic tree cut, different colors represent different modules. c: module division after merging similar modules.

2.3 抗病相關(guān)特異性模塊鑒定

在18個基因模塊中, 有5個與黃萎病菌侵染存在高度特異性(||>0.70,<0.001)。其中, black (=0.72,=4e-04)、mediumpurple3 (=0.78,=5e-05)、darkolivegreen (=0.76,=9e-05)、plum3 (=0.97,=6e-12)模塊分別與侵染2 h、6 h、48 h、72 h時間點正相關(guān); mediumpurple2 (=-0.82,=8e-06)與侵染2 h時間點負(fù)相關(guān)(圖4)。

2.4 特異性模塊富集分析

GO通路可分為生物過程(Biological Process, BP)、分子功能(Molecular Function, MF)及細(xì)胞組分(Cellular Component, CC) 3類。分析結(jié)果表明, black和mediumpurple3模塊分別富集到41個和36個生物過程, 其中black模塊主要富集到氮化合物轉(zhuǎn)運(GO:0071705)、蛋白定位(GO:0045184)、蛋白轉(zhuǎn)運(GO:0015031)等調(diào)控通路, 以及一些信號相關(guān)通路,如信號轉(zhuǎn)導(dǎo)(GO:0007165)、信號轉(zhuǎn)導(dǎo)調(diào)控(GO:1902531)等(圖5-a); mediumpurple3模塊主要富集到一些代謝相關(guān)通路, 如氨基酸代謝(GO:0006520)、硫化物代謝(GO:0006790)、甘氨酸代謝(GO:0006544)等生物學(xué)過程(附表1)。此外, black、mediumpurple3、darkolivegreen、plum3及mediumpurple2模塊分別富集到21、7、7、1和2個分子功能, 如蛋白質(zhì)轉(zhuǎn)運蛋白活性(GO:0008565)、鈣離子結(jié)合(GO:0005509)(附表2); black和mediumpurple2模塊分別富集到44和4個細(xì)胞組分, 最顯著的分別為膜衣(GO:0030117)和蛋白酶體核心復(fù)合物(GO:0005839)(附表3)。

圖4 模塊與性狀相關(guān)性熱圖

每行代表一個模塊, 每列代表一種性狀。矩形框里的數(shù)字代表模塊與性狀之間的相關(guān)系數(shù)及相應(yīng)值?？共∠嚓P(guān)特異性模塊用紅色標(biāo)示, 其基因及轉(zhuǎn)錄因子數(shù)量標(biāo)于左側(cè)模塊矩形框內(nèi)。

Each row corresponds to a module, and each column corresponds to a trait. The correlation coefficient and the correspondingvalue are shown in each cell. The specific modules which significantly associated withinfection are colored in red, and the corresponding numbers of genes and transcription factors are shown in the left cells.

KEGG富集分析表明, black、mediumpurple3和darkolivegreen模塊分別富集到8、10和1個KEGG代謝通路, 其中black模塊主要富集到mRNA監(jiān)測通路、內(nèi)吞作用和蛋白質(zhì)輸出等通路(圖5-b); mediumpurple3主要富集到氨基酸生物合成、檸檬酸循環(huán)和碳代謝等通路; darkolivegreen富集到黃酮類化合物生物合成代謝通路(附表4)。

圖5 Black模塊的GO和KEGG富集分析

縱坐標(biāo)代表GO條目或KEGG代謝通路, 橫坐標(biāo)代表富集到的基因在模塊內(nèi)所占比率。點的大小代表富集到的基因數(shù)量, 點的顏色代表多重校驗后的值大小。

The vertical axis represents the enriched GO term or KEGG pathway, and the horizontal axis represents the ratio of enriched genes in the module. The point size represents the gene number enriched in the pathway, and the color represents the-value corrected for multiple testing.

2.5 核心基因鑒定及基因互作網(wǎng)絡(luò)構(gòu)建

核心基因通常指模塊內(nèi)具有高連通性的基因, 本研究選取每個模塊中kME值(Eigengene Connectivity)最高的前5個基因做為核心基因, 并利用核心基因及其互作基因繪制基因互作網(wǎng)絡(luò)圖(圖6)。通過與擬南芥進(jìn)行同源比對, 對核心基因及其互作基因進(jìn)行了功能注釋(表2, 附表5)。此外, 本研究選取了9個核心基因進(jìn)行了qRT-PCR分析, 結(jié)果表明其表達(dá)模式與轉(zhuǎn)錄組數(shù)據(jù)基本一致(圖7)。

3 討論

本研究以黃萎病菌侵染下海島棉幼苗根系RNA-seq數(shù)據(jù)為基礎(chǔ), 對其進(jìn)行了差異表達(dá)分析, 共鑒定出22,850個差異表達(dá)基因。Chen等[18]利用該轉(zhuǎn)錄組數(shù)據(jù)共得到17,517個差異表達(dá)基因, 結(jié)果產(chǎn)生差異的主要原因是在分析過程中使用了不同的參考基因組及分析軟件。Chen等[18]以雷蒙德氏棉()基因組為參考基因組, 比對率約為76%, 本研究以海島棉基因組為參考基因組, 比對率約為94%, 因此比對效果更好。Zhang等[28]通過分析黃萎病菌侵染24 h下陸地棉()幼苗根的轉(zhuǎn)錄組數(shù)據(jù), 共鑒定到4794個差異表達(dá)基因, 其中一些基因在本研究中具有相似的差異表達(dá)模式(如、、等), 但也有部分基因在2個研究中表現(xiàn)出相反的上下調(diào)變化趨勢(如、、等), 暗示這些基因在海島棉和陸地棉的抗病過程中可能發(fā)揮不同作用。

通過WGCNA方法, 本研究共鑒定到5個(black、mediumpurple3、darkolivegreen、plum3、mediumpurple2)抗病相關(guān)特異性模塊。富集分析結(jié)果表明, 特異性模塊可得到具有相關(guān)生物學(xué)意義的功能富集和代謝通路富集結(jié)果。如GO富集分析中, black模塊可富集到刺激響應(yīng)(GO:0051716)、激素介導(dǎo)的信號通路(GO:0009755)、激素刺激應(yīng)答(GO:0032870)等抗逆相關(guān)細(xì)胞過程; plum3模塊富集到鈣離子結(jié)合分子功能(GO:0005509), 而研究表明鈣離子是植物防御應(yīng)答過程中一類非常重要的第二信使[29]; KEGG富集分析中, darkolivegreen模塊可富集到黃酮類化合物生物合成代謝通路, 此類化合物是植物中所特有的一類多功能復(fù)合物, 在植株抵御生物/非生物脅迫過程中發(fā)揮重要作用[30]。

圖6 抗病相關(guān)特異性模塊內(nèi)的基因共表達(dá)網(wǎng)絡(luò)

點的大小和深淺代表該基因在網(wǎng)絡(luò)中連通性的高低。轉(zhuǎn)錄因子用三角形表示, 其他基因用圓形表示。a: Black模塊內(nèi)的基因共表達(dá)網(wǎng)絡(luò)。b: Mediumpurple3模塊內(nèi)的基因共表達(dá)網(wǎng)絡(luò)。c: Darkolivegreen模塊內(nèi)的基因共表達(dá)網(wǎng)絡(luò)。d: Plum3模塊內(nèi)的基因共表達(dá)網(wǎng)絡(luò)。e: Mediumpurple2模塊內(nèi)的基因共表達(dá)網(wǎng)絡(luò)。

The genes with higher connectivity in the corresponding network are shown with larger circle sizes and darker colors. The transcription factors are indicated by triangles, with other genes being indicated by circles. a: Gene co-expression network of the black module. b: Gene co-expression network of the mediumpurple3 module. c: Gene co-expression network of the darkolivegreen module. d: Gene co-expression network of the plum3 module. e: Gene co-expression network of the mediumpurple2 module.

圖7 核心基因表達(dá)分析

a: qRT-PCR分析, 上方垂直線代表3次重復(fù)的標(biāo)準(zhǔn)差。b: RNA-seq表達(dá)熱圖。

a: expression analysis by qRT-PCR, the bars indicate standard deviation of three replications. b: heat map of the hub genes.

*:< 0.05, **:< 0.01.

表2 抗病相關(guān)特異性模塊中核心基因的功能注釋

通過計算模塊內(nèi)基因的連通性, 可推測該基因在網(wǎng)絡(luò)中的位置及重要程度。結(jié)合已經(jīng)報道過的棉花中抗黃萎病相關(guān)基因發(fā)現(xiàn), 這些基因在相應(yīng)模塊內(nèi)都具有較高的連通性。如black模塊中的基因(, kME=0.89), 可通過激活基因的表達(dá)介導(dǎo)棉花抗病過程[31]; mediumpurple2模塊中的(, kME=0.89)為一類NBS-LRR (nucleotide-binding site–leucine rich repeat)基因, 在棉花抗黃萎病侵染過程中發(fā)揮重要作用[32]。此外, darkolivegreen模塊中的、plum3模塊中的在相應(yīng)模塊內(nèi)都具有較高連通性,其在陸地棉中的同源基因及均為抗病相關(guān)基因[33-34]。其他非抗病相關(guān)特異性模塊中也包含一些高連通性的抗病基因, 如turquoise模塊中的(, kME= 0.90)[35]、(, kME=0.95)[36]及mediumorchid模塊中的(, kME=0.91)[37]基因等。

選取特異性模塊中連通性最高的前5個基因為核心基因, 推測其可能在抗病過程中發(fā)揮重要作用。這些基因在棉花中的功能大多尚不明確, 而它們在擬南芥中的同源基因部分已報道為抗逆相關(guān)基因。如plum3模塊中的為一個ERF轉(zhuǎn)錄因子家族基因, 其在擬南芥中的同源基因為, 該基因可通過調(diào)控JA/ET信號通路響應(yīng)病原菌脅迫[38]; darkolivegreen模塊中的同源基因參與調(diào)控ABA信號轉(zhuǎn)導(dǎo), 在植株非生物脅迫應(yīng)答中發(fā)揮作用[39]; black模塊中在擬南芥中的同源基因為, 該基因通過激活G6PD蛋白調(diào)節(jié)細(xì)胞內(nèi)氧化還原反應(yīng)平衡, 從而在植株響應(yīng)逆境脅迫中發(fā)揮重要作用[40]。

與核心基因連接度較高的基因, 同樣可能在逆境脅迫應(yīng)答中發(fā)揮作用。如在darkolivegreen模塊中,的高連接度基因在擬南芥中的同源基因為, 該基因可通過調(diào)節(jié)活性氧含量、細(xì)胞死亡過程及基因的表達(dá)參與響應(yīng)病原菌侵染脅迫[41]; 在plum3模塊中,的高連接度基因在擬南芥中的同源基因編碼一個細(xì)胞分裂素信號通路相關(guān)轉(zhuǎn)錄因子ARR12, 參與調(diào)控植株的干旱脅迫應(yīng)答[42]; 在black模塊中,的高連接度基因在擬南芥中的同源基因編碼PAD4蛋白, 可與另一個抗病相關(guān)蛋白EDS1組成復(fù)合體共同調(diào)控SA信號通路, 進(jìn)而在病原菌脅迫應(yīng)答中發(fā)揮作用[43]。此外, 結(jié)合差異表達(dá)分析結(jié)果表明, 以上基因基本都發(fā)生了顯著上下調(diào), 從而在表達(dá)水平上也說明它們可能在抗病過程中發(fā)揮作用。

本研究構(gòu)建的網(wǎng)絡(luò)中核心基因在棉花中的具體生物學(xué)功能目前大多尚不明確, 后續(xù)可進(jìn)一步通過過量表達(dá)、VIGS、基因敲除等生物技術(shù)手段對其開展深入研究。

4 結(jié)論

通過構(gòu)建權(quán)重基因共表達(dá)網(wǎng)絡(luò)共鑒定到5個抗病相關(guān)特異性模塊, 并揭示了它們的生物學(xué)意義。以連通性為指標(biāo), 揭示了特異性模塊中可能在抗病過程中發(fā)揮重要作用的關(guān)鍵基因。本研究結(jié)果為進(jìn)一步理解棉花抗病過程的分子機制提供了理論指導(dǎo)。

附表 請見網(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。

[1] Wang J Y, Cai Y, Gou J Y, Mao Y B, Xu Y H, Jiang W H, Chen X Y. VdNEP, an elicitor from, induces cotton plant wilting., 2004, 70: 4989–4995.

[2] Li Z F, Liu Y J, Feng Z L, Feng H J, Klosterman S J, Zhou F F, Zhao L H, Shi Y Q, Zhu H Q., encoding CYC8 glucose repression mediator protein, is required for microsclerotia formation and full virulence in, 2015, 10: e0144020.

[3] Xu J, Wang X Y, Li Y Q, Zeng J G, Wang G L, Deng C Y, Guo W Z. Host-induced gene silencing of a regulator of G protein signaling gene () confers resistance towilt in cotton, 2018, 16: 1629–1643.

[4] 朱荷琴, 李志芳, 馮自力, 馮鴻杰, 魏鋒, 趙麗紅, 師勇強, 劉世超, 周京龍. 我國棉花黃萎病研究十年回顧及展望棉花學(xué)報, 2017, 29(增刊): 37–50. Zhu H Q, Li Z F, Feng Z L, Feng H J, Wei F, Zhao L H, Shi Y Q, Liu S C, Zhou J L. Overview of cottonwilt research over the past decade in china and its prospect in future., 2017, 29(suppl): 37–50 (in Chinese with English abstract).

[5] Kawchuk L M, Hachey J, Lynch D R, Kulcsar F, van Rooijen G, Waterer D R, Robertson A, Kokko E, Byers R, Howard R J, Fischer R, Prüfer D. Tomatodisease resistance genes encode cell surface-like receptors, 2001, 98: 6511–6515.

[6] 劉琳琳, 張文文, 周易, 苗玉煥, 許蓮, 劉敏, 張坤, 張獻(xiàn)龍, 朱龍付. 棉花與番茄抗棉花黃萎病不依賴于中國科學(xué): 生命科學(xué), 2014, 44: 803–814. Liu L L, Zhang W W, Zhou Y, Miao Y H, Xu L, Liu M, Zhang K, Zhang X L, Zhu L F. Resistance of cotton and tomato tofrom cotton is independent on., 2014, 44: 803–814 (in Chinese).

[7] Zhang Y, Wang X F, Ding Z G, Ma Q, Zhang G R, Zhang S L, Li Z K, Wu L Q, Zhang G Y, Ma Z Y. Transcriptome profiling ofinoculated withprovides a resource for cotton improvement, 2013, 14: 637.

[8] Zuo K J, Qin J, Zhao J Y, Ling H, Zhang L D, Cao Y F, Tang K X. Over-expression GbERF2 transcription factor in tobacco enhances brown spots disease resistance by activating expression of downstream genes, 2007, 391: 80–90.

[9] Zhang Y, Wang X F, Li Y Y, Wu L Z, Zhou H M, Zhang G Y, Ma Z Y. Ectopic expression of a novel Ser/Thr protein kinase from cotton (), enhances resistance toinfection and oxidative stress in, 2013, 32: 1703–1713.

[10] Li T G, Wang B L, Yin C M, Zhang D D, Wang D, Song J, Zhou L, Kong Z Q, Klosterman S J, Li J J, Adamu S, Liu T L, Subbarao K V, Chen J Y, Dai X F. TheTIR-NBS-LRR genemediates resistance againstwilt, 2019, 20: 857–876.

[11] Munis M F H, Tu L L, Deng F L, Tan J F, Xu L, Xu S C, Long L, Zhang X L. A thaumatin-like protein gene involved in cotton fiber secondary cell wall development enhances resistance againstand other stresses in transgenic tobacco, 2010, 393: 38–44.

[12] Mo H J, Wang X F, Zhang Y, Zhang G Y, Zhang J F, Ma Z Y. Cotton polyamine oxidase is required for spermine and camalexin signalling in the defence response to, 2015, 83: 962–975.

[13] He X, Zhu L F, Wassan G M, Wang Y J, Miao Y H, Shaban M, Hu H Y, Sun H, Zhang X L. GhJAZ2 attenuates cotton resistance to biotic stresses via the inhibition of the transcriptional activity of GhbHLH171, 2018, 19: 896–908.

[14] Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis, 2008, 9: 559.

[15] Tan M P, Cheng D, Yang Y N, Zhang G Q, Qin M J, Chen J, Chen Y H, Jiang M Y. Co-expression network analysis of the transcriptomes of rice roots exposed to various cadmium stresses reveals universal cadmium-responsive genes, 2017, 17: 194.

[16] Zou X Y, Liu A Y, Zhang Z, Ge Q, Fan S M, Gong W K, Li J W, Gong J W, Shi Y Z, Tian B M, Wang Y L, Liu R X, Lei K, Zhang Q, Jiang X, Feng Y L, Zhang S Y, Jia T T, Zhang L P, Yuan Y L, Shang H H. Co-expression network analysis and hub gene selection for high-quality fiber in upland cotton () using RNA sequencing analysis, 2019, 10: 119.

[17] 楊宇昕, 桑志勤, 許誠, 代文雙, 鄒棖. 利用WGCNA進(jìn)行玉米花期基因共表達(dá)模塊鑒定作物學(xué)報, 2019, 45: 161–174. Yang Y X, Sang Z Q, Xu C, Dai W S, Zou C. Identification of maize flowering gene co-expression modules by WGCNA., 2019, 45: 161–174 (in Chinese with English abstract).

[18] Chen J Y, Huang J Q, Li N Y, Ma X F, Wang J L, Liu C, Liu Y F, Liang Y, Bao Y M, Dai X F. Genome-wide analysis of the gene families of resistance gene analogues in cotton and their response towilt, 2015, 15: 148.

[19] Bolger A M, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data, 2014, 30: 2114–2120.

[20] Wang M J, Tu L L, Yuan D J, Zhu D, Shen C, Li J Y, Liu F Y, Pei L L, Wang P C, Zhao G N, Ye Z X, Huang H, Yan F L, Ma Y Z, Zhang L, Liu M, You J Q, Yang Y C, Liu Z P, Huang F, Li B Q, Qiu P, Zhang Q H, Zhu L F, Jin S X, Yang X Y, Min L, Li G L, Chen L L, Zheng H K, Lindsey K, Lin Z X, Udall J A, Zhang X L. Reference genome sequences of two cultivated allotetraploid cottons,and, 2019, 51: 224–229.

[21] Kim D, Langmead B, Salzberg S L. HISAT: a fast spliced aligner with low memory requirements, 2015, 12: 357.

[22] Liao Y, Smyth G K, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features, 2014, 30: 923–930.

[23] Love M I, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2, 2014, 15: 550.

[24] Yu G C, Wang L G, Han Y Y, He Q Y. clusterProfiler: an R package for comparing biological themes among gene clusters, 2012, 16: 284–287.

[25] Jin J P, Tian F, Yang D C, Meng Y Q, Kong L, Luo J C, Gao G. PlantTFDB 4.0: toward a central hub for transcription factors and regulatory interactions in plants, 2016, 45: D1040–D1045.

[26] Livak K J, Schmittgen T D. Analysis of relative gene expression data using real-time quantitative PCR and the 2?ΔΔCTmethod, 2001, 25: 402–408.

[27] Erpen L, Devi H S, Grosser J W, Dutt M. Potential use of the DREB/ERF, MYB, NAC and WRKY transcription factors to improve abiotic and biotic stress in transgenic plants, 2018, 132: 1–25.

[28] Zhang W W, Zhang H C, Liu K, Jian G L, Qi F J, Si N. Large-scale identification ofgenes asso-ciated withby comparative transcriptomic and reverse genetics analysis, 2017, 12: e0181609.

[29] Lecourieux D, Ranjeva R, Pugin A. Calcium in plant defence-signalling pathways, 2006, 171: 249–269.

[30] Liu J Y, Osbourn A, Ma P D. MYB transcription factors as regulators of phenylpropanoid metabolism in plants, 2015, 8: 689–708.

[31] Li C, He X, Luo X Y, Xu L, Liu L L, Min L, Jin L, Zhu L F, Zhang X L. Cotton WRKY1 mediates the plant defense- to-development transition during infection of cotton byby activatingexpression, 2014, 166: 2179–2194.

[32] Yang J, Ma Q, Zhang Y, Wang X F, Zhang G Y, Ma Z Y. Molecular cloning and functional analysis of, a gene inthat plays an important role in conferring resistance towilt, 2016, 575: 687–694.

[33] Yang J, Wang G N, Ke H F, Zhang Y, Ji L L, Huang L Z, Zhang C Y, Wang X F, Ma Z Y. Genome-wide identification of cyclophilin genes inand functional characterization of a CYP with antifungal activity against, 2019, 19: 272.

[34] Xiong X P, Sun S C, Li Y J, Zhang X Y, Sun J, Xue F. The cotton WRKY transcription factornegatively regulates the defense response against, 2019, 7: 393–402.

[35] Cheng H Q, Han L B, Yang C L, Wu X M, Zhong N Q, Wu J H, Wang F X, Wang H Y, Xia G X. The cotton MYB108 forms a positive feedback regulation loop with CML11 and participates in the defense response againstinfection, 2016, 67: 1935–1950.

[36] Li Y B, Han L B, Wang H Y, Zhang J, Sun S T, Feng D Q, Yang C L, Sun Y D, Zhong N Q, Xia G X. The thioredoxin GbNRX1 plays a crucial role in homeostasis of apoplastic reactive oxygen species in response toinfection in cotton, 2016, 170: 2392–2406.

[37] Duan X P, Zhang Z D, Wang J, Zuo K J. Characterization of a novel cotton subtilase genein response to extracellular stimulations and its role inresistance, 2016, 11: e0153988.

[38] Moffat C S, Ingle R A, Wathugala D L, Saunders N J, Knight H, Knight M R. ERF5 and ERF6 play redundant roles as positive regulators of JA/Et-mediated defense againstin, 2012, 7: e35995.

[39] Kumar D, Kumar R, Baek D W, Hyun T K, Chung W S, Yun D J, Kim J Y.functions as a positive regulator in plant responses to ABA, 2017, 10: 223–243.

[40] Dal Santo S, Stampfl H, Krasensky J, Kempa S, Gibon Y, Petutschnig E, Rozhon W, Heuck A, Clausen T, Jonak C. Stress-induced GSK3 regulates the redox stress response by phosphorylating glucose-6-phosphate dehydrogenase in, 2012, 24: 3380–3392.

[41] Pan Q N, Cui B M, Deng F Y, Quan J L, Loake G J, Shan W X.encodes a novel endoplasmic reticulum (ER)-localized protein inand negatively regulates resistance against biotrophic pathogens, 2016, 209: 1641–1654.

[42] Nguyen K H, Ha C V, Nishiyama R, Watanabe Y, Leyva-González M A, Fujita Y, Tran U T, Li W Q, Tanaka M, Seki M, Schaller G E, Herrera-Estrella L, Tran L-S P.type B cytokinin response regulators ARR1, ARR10, and ARR12 negatively regulate plant responses to drought, 2016, 113: 3090–3095.

[43] Cui H T, Qiu J D, Zhou Y, Bhandari D D, Zhao C H, Bautor J, Parker J E. Antagonism of transcription factor MYC2 by EDS1/PAD4 complexes bolsters salicylic acid defense ineffector-triggered immunity, 2018, 11: 1053–1066.

Identification of co-expressed modules of cotton genes responding toinfection by WGCNA

FU Ming-Chuan*, LI Hao, CHEN Yi-Zhen, LIU Zhan-Ji, LIU Ren-Zhong, and WANG Li-Guo

1Cotton Research Center of Shandong Academy of Agricultural Sciences / Key Laboratory of Cotton Breeding and Cultivation in Huang-Huai-Hai Plain, Ministry of Agriculture, Jinan 250100, Shandong, China

Weighted gene co-expression network analysis (WGCNA) is a classic systematic biological method, which can be used to identify co-expressed modules, investigate relationships between modules and specific traits, and screen hub genes in the networks.wilt, caused by the fungal pathogen, can cause severe fibre quality reduction and yield loss of cotton. Studying on cotton genes and molecular mechanisms related to defense responses againstcan shed light on cotton breeding. In this study, 21 transcriptome data ofseedling roots infected bywere used to investigate differentially expressed genes (DEGs). By filtering out the genes with low variation, 35,647 genes were selected for WGCNA. In total, 22,850 DEGs were identified underinfection, with 4685 in common at all inoculated time points. Co-expression network analysis identified 18 modules, in which five modules significantly associated withinfection (black, mediumpurple3, darkolivegreen, plum3 positively correlated with the 2 h, 6 h, 48 h, and 72 h inoculated time points, respectively; mediumpurple2 negatively correlated with the 2 h inoculated time point). GO and KEGG enrichment analysis were performed on these five specific modules, which could be enriched in GO terms and metabolic pathways, such as cellular response to stimulus, calcium ion binding and flavonoid biosynthesis. The hub genes were screened by calculating gene connectivity in the corresponding networks, and they may play important roles in the resistance against biotic/abiotic stresses. In summary, by WGCNA, a few defense-related co-expressed modules and hub genes were identified. The results of this study would be beneficial for further understanding of the molecular mechanisms of pathogen resistance in cotton, and provide new gene resources for cotton breeding in the future.

cotton;wilt; WGCNA; hub gene

10.3724/SP.J.1006.2020.94124

本研究由山東省農(nóng)業(yè)科學(xué)院農(nóng)業(yè)科技創(chuàng)新工程(CXGC2018B01, CXGC2018E06)項目資助。

This study was supported by the Agricultural Science and Technology Innovation Project of Shandong Academy of Agricultural Sciences (CXGC2018B01, CXGC2018E06).

2019-08-23;

2020-01-15;

2020-02-18.

傅明川, E-mail: fumingchuan@shandong.cn

URL: http://kns.cnki.net/kcms/detail/11.1809.s.20200218.1558.002.html