基于加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析探究桑樹(shù)光合相關(guān)基因的調(diào)控模式

李 勇，黃 錦，于 翠，莫榮利，朱志賢，董朝霞，胡興明，鄧 文

（湖北省農(nóng)業(yè)科學(xué)院經(jīng)濟(jì)作物研究所，武漢 430064）

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（Weighted correlation networkanalysis，WGCNA）主要是描述多個(gè)樣品之間基因關(guān)聯(lián)模式的系統(tǒng)生物學(xué)方法。根據(jù)基因的表達(dá)量鑒定高度協(xié)同變化的基因模塊，表達(dá)模式相似的基因被聚為一個(gè)模塊，通過(guò)分析基因模塊連通性和進(jìn)一步探索基因模塊與特定性狀或表型之間的關(guān)聯(lián)性，來(lái)挖掘表達(dá)網(wǎng)絡(luò)中的關(guān)鍵基因［1］。普通轉(zhuǎn)錄組分析只能關(guān)注差異表達(dá)基因信息，而WGCNA可以對(duì)上萬(wàn)個(gè)變化的基因或全部基因信息進(jìn)行聚類，并與目標(biāo)表型進(jìn)行關(guān)聯(lián)分析，來(lái)識(shí)別并候選重要的基因模塊，進(jìn)而關(guān)注重要模塊內(nèi)與目標(biāo)表型相關(guān)的關(guān)鍵基因。共表達(dá)網(wǎng)絡(luò)分析已經(jīng)成功應(yīng)用于人和動(dòng)物中復(fù)雜性狀和疾病的研究［2，3］。同時(shí)在植物中也有廣泛應(yīng)用，如通過(guò)WGCNA分析擬南芥根中基因表達(dá)和表型的動(dòng)態(tài)變化［4］，在番茄代謝中的研究［5］以及轉(zhuǎn)錄組數(shù)據(jù)與WGCNA方法結(jié)合研究玉米子粒大小、產(chǎn)量和株高等性狀和組織特異性模塊［6-8］。WGCNA通常適用于多樣本數(shù)據(jù)模式分析，一般要求的樣本數(shù)至少為15個(gè)，樣本數(shù)越多，結(jié)果越穩(wěn)定。本試驗(yàn)基于前期已有的18份樣本轉(zhuǎn)錄組數(shù)據(jù)和3個(gè)桑樹(shù)品種光合作用的表型數(shù)據(jù)以及樣品數(shù)據(jù)［9］，通過(guò)WGCNA分析研究不同桑樹(shù)品種不同時(shí)間段基因的表達(dá)變化，從而更深入探索基因模塊與目標(biāo)表型間的關(guān)聯(lián)性，進(jìn)而研究不同桑樹(shù)品種不同時(shí)間點(diǎn)光合作用相關(guān)基因的調(diào)控網(wǎng)絡(luò)，找出導(dǎo)致不同桑樹(shù)品種光合作用差異的關(guān)鍵基因，挖掘影響光合作用的重要基因，為后續(xù)驗(yàn)證和確定與光合作用差異相關(guān)的關(guān)鍵基因提供更夯實(shí)的證據(jù)，可進(jìn)一步闡明影響桑樹(shù)產(chǎn)量差異的光合作用機(jī)制。

1 材料與方法

1.1 數(shù)據(jù)材料

基于團(tuán)隊(duì)前期轉(zhuǎn)錄組數(shù)據(jù)和鄂桑1號(hào)（E1）、鄂桑2號(hào)（E2）、湖桑32號(hào)（H32）3個(gè)桑樹(shù)品種，10時(shí)（T10）和12時(shí)（T12）2個(gè)時(shí)間點(diǎn)，桑樹(shù)的光合生理參數(shù)、光合作用酶活性、過(guò)氧化防御系統(tǒng)參數(shù)的29個(gè)表型數(shù)據(jù)作為本試驗(yàn)的性狀表型；其中，光合生理參數(shù)：凈光合速率（Pn）、氣孔導(dǎo)度（Gs）、胞間二氧化碳濃度（Ci）、蒸騰速率（Tr）、表觀量子效率（AQY）、光補(bǔ)償點(diǎn)（LCP）、飽和光強(qiáng)（LSP）、羧化效率（CE）、CO2補(bǔ)償點(diǎn)（CCP）與飽和點(diǎn)（CSP）、光合能力（Pm）；葉綠素?zé)晒鈪?shù)：PSII最大量子效率（Fv/Fm）、PSII實(shí)際量子效率（ΦPSII）、表觀電子傳遞速率（ETR）、光化學(xué)淬滅系數(shù)（qP）和非光化學(xué)淬滅（NPQ）；生理生化指標(biāo)：超氧化物酶（SOD）活性、過(guò)氧化物酶（POD）酶活、丙二醛（MDA）含量、過(guò)氧化氫酶（CAT）酶活、RuBP羧化酶活性；葉綠素含量等。對(duì)所有轉(zhuǎn)錄組數(shù)據(jù)中的所有基因以及光合作用的表型數(shù)據(jù)進(jìn)行WGCNA分析。

1.2 基因間相關(guān)系數(shù)的計(jì)算

WGCNA又叫權(quán)重基因共表達(dá)網(wǎng)絡(luò)分析。通過(guò)基因表達(dá)模式構(gòu)建基因共表達(dá)網(wǎng)絡(luò)的前提需要基因之間的相關(guān)系數(shù)及構(gòu)建共表達(dá)網(wǎng)絡(luò)的相似矩陣。對(duì)于無(wú)尺度網(wǎng)絡(luò)，基因i和基因j的相關(guān)系數(shù)為Sij=∣cor（i，j）∣；對(duì)于有尺度網(wǎng)絡(luò)，基因i和基因j的相關(guān)系數(shù)為Sij=∣1+cor（i，j）/2∣。根據(jù)兩兩基因的相似度，計(jì)算基因相似表達(dá)矩陣S=［Sij］。

1.3 軟閾值的選擇

一般需要設(shè)置閾值來(lái)篩選2個(gè)基因之間是否具有相似的表達(dá)模式，高于閾值的可認(rèn)為2個(gè)基因相似。WGCNA其實(shí)是采用軟閾值的方法，對(duì)基因表達(dá)值之間的相關(guān)系數(shù)取n次冪，使網(wǎng)絡(luò)中基因之間的連接服從無(wú)尺度網(wǎng)絡(luò)分布，這是和普通聚類的最大不同，直接結(jié)果是把基因間相關(guān)性的強(qiáng)弱差異放大，這種算法更具有生物學(xué)意義［10］。對(duì)每?jī)蓪?duì)基因（i，j）之間的相關(guān)系數(shù)的冪取某特定的β值，β值即為軟閾值，以此來(lái)計(jì)算所有基因之間的相關(guān)性，也就是adjacency矩陣：aij=|cor（i，j）|β。由這些相關(guān)性系數(shù)可以構(gòu)建網(wǎng)絡(luò)，其中，基因作為網(wǎng)絡(luò)中的節(jié)點(diǎn)（nodes），而aij作為邊的權(quán)重，高相關(guān)性代表強(qiáng)連接，反之亦然。每個(gè)基因連接度的大小反映了與該基因相連的基因的多少。adjacency的結(jié)果直接取決于β的取值，其取值的高低直接影響模塊的構(gòu)建和模塊內(nèi)基因的劃分。根據(jù)接近無(wú)尺度網(wǎng)絡(luò)（Scale-free network）的最低值來(lái)確定β值，所以需要選擇合適的構(gòu)建無(wú)尺度網(wǎng)絡(luò)的β值。WGCNA給出了如何選擇這個(gè)閾值的方法，其函數(shù)為pickSoft-Threshold。

1.4 節(jié)點(diǎn)之間的相關(guān)性計(jì)算

基因之間的相關(guān)性包含了直接相關(guān)和間接相關(guān)。WGCNA中通過(guò)拓?fù)渲丿B（Topological overlapmatrix，TOM）來(lái)計(jì)算2個(gè)基因之間相似性。也就是說(shuō)，如果基因i和j有很多相同的鄰接基因，那么TOM（i，j）就很高，這意味著基因有相似的表達(dá)模式。TOM（i，j）=0，表示基因i和j的網(wǎng)絡(luò)沒(méi)有共同的鄰接基因，而TOM（i，j）=1，表示基因i和j有完全一樣的網(wǎng)絡(luò)鄰接基因。TOM作為相似度的度量可以被轉(zhuǎn)換為相異度，disTOM=1-TOM。

1.5 構(gòu)建分層聚類樹(shù)

根據(jù)基因間的拓?fù)渲丿B計(jì)算基因之間的相異性系數(shù)，并以相異性系數(shù)為基礎(chǔ)構(gòu)建基因間的系統(tǒng)聚類樹(shù)，在系統(tǒng)聚類樹(shù)中，表達(dá)模式相似的基因?qū)儆谝粋€(gè)樹(shù)分支，本分析通過(guò)動(dòng)態(tài)剪切樹(shù)法將所有差異基因根據(jù)不同的表達(dá)模式分割為不同的模塊［1］。每一個(gè)模塊內(nèi)的基因具有表達(dá)模式的連通性和相似性。

1.6 基因模塊與表型的關(guān)聯(lián)

將基因模塊與表型數(shù)據(jù)進(jìn)行關(guān)聯(lián)，可檢測(cè)與光合作用顯著相關(guān)的模塊，再通過(guò)模塊內(nèi)的關(guān)鍵基因篩選與光合作用相關(guān)的重要基因；計(jì)算每一個(gè)模塊內(nèi)的模塊特征值（Module eigengene，ME），再計(jì)算模塊特征向量與表型值的相關(guān)系數(shù)，可以很直觀地反映出每一個(gè)模塊與光合作用的相關(guān)性。

1.7 模塊內(nèi)基因的GO和KEGG富集分析

對(duì)分析獲得的模塊內(nèi)基因進(jìn)行GO和KEGG富集分析，探索這些模塊內(nèi)的基因參與哪些生物學(xué)功能和富集在哪些通路，可挖掘哪些基因與模塊內(nèi)基因具有高度的連通性，進(jìn)一步探索光合作用相關(guān)基因的調(diào)控網(wǎng)絡(luò)。

2 結(jié)果與分析

2.1 用于WGCNA分析的基因和表型矩陣

通過(guò)對(duì)鄂桑1號(hào)（E1）、鄂桑2號(hào)（E2）和湖桑32號(hào)（H32）3個(gè)桑樹(shù)品種，10時(shí)和12時(shí)2個(gè)時(shí)間點(diǎn)的18份樣本的轉(zhuǎn)錄組分析，得到所有差異表達(dá)基因的并集，累計(jì)5 416個(gè)差異表達(dá)基因。在用于構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)的基因表達(dá)矩陣中，以FPKM值為分析標(biāo)準(zhǔn)值，以表達(dá)差異倍數(shù)FoldChange≥2且錯(cuò)誤發(fā)現(xiàn)率FDR＜0.05作為統(tǒng)計(jì)標(biāo)準(zhǔn)，鑒定出DEGs共計(jì)5 291個(gè)。用于基因模塊與表型關(guān)聯(lián)的性狀矩陣來(lái)自3個(gè)桑樹(shù)品種不同的光合作用參數(shù)和生理生化指標(biāo)，另外加上3個(gè)桑樹(shù)品種、2個(gè)時(shí)間點(diǎn)。累計(jì)獲取29個(gè)表型數(shù)據(jù)，并以此作性狀矩陣與基因模塊進(jìn)行關(guān)聯(lián)，進(jìn)一步挖掘與目標(biāo)性狀相關(guān)的光合作用基因。

2.2 軟閾值確定

按照無(wú)尺度網(wǎng)絡(luò)的標(biāo)準(zhǔn)選擇合適的軟閾值（β），軟閾值（β）與無(wú)尺度網(wǎng)絡(luò)評(píng)價(jià)系數(shù)（R2）的關(guān)系以及軟閾值（β）與平均連通性的關(guān)系如圖1所示。構(gòu)建基因加權(quán)共表達(dá)網(wǎng)絡(luò)，確定軟閾值時(shí)，一般選擇無(wú)尺度網(wǎng)絡(luò)評(píng)價(jià)系數(shù)（R2）在0.9左右對(duì)應(yīng)的軟閾值，本研究中以無(wú)尺度網(wǎng)絡(luò)評(píng)價(jià)系數(shù)（0.85）對(duì)應(yīng)的軟閾值（16）來(lái)構(gòu)建基因共表達(dá)網(wǎng)絡(luò)。

圖1 軟閾值確定

2.3 基因聚類和基因模塊的確定

通過(guò)計(jì)算基因之間的相關(guān)系數(shù)、連接系數(shù)和拓?fù)渲丿B，得到每個(gè)基因之間的相異性系數(shù)，利用基因間相異性系數(shù)構(gòu)建基因的系統(tǒng)聚類數(shù)，最后根據(jù)動(dòng)態(tài)剪切樹(shù)算法構(gòu)建劃分基因模塊，如圖2、表1所示，本研究中，一共獲得10個(gè)基因模塊，所有模塊累計(jì)有5 291個(gè)基因，其中，紅色模塊包含的基因數(shù)目最少，只有39個(gè)基因，品紅基因模塊包含的基因數(shù)量最多，有2 414個(gè)基因。

表1 基因模塊和對(duì)應(yīng)的基因數(shù)目

圖2 基因模塊聚類

2.4 基因模塊與表型值的關(guān)聯(lián)

為了進(jìn)一步探索基因模塊與不同性狀的關(guān)系，使用每個(gè)基因模塊的模塊特征值（ME）與本研究中的29個(gè)性狀進(jìn)行關(guān)聯(lián)分析，來(lái)確定某一個(gè)模塊基因是否與某個(gè)特定性狀高度相關(guān)聯(lián)，本試驗(yàn)重點(diǎn)關(guān)注與性狀正相關(guān)的模塊基因。如圖3所示，發(fā)現(xiàn)MEblue和MEpink模塊分別與H32、AQY、LCP、CCP、CSP、NPQ、ETR和POD 8個(gè)性狀高度正相關(guān)，其中,MEblue模塊與H32、NPQ和ETR呈正相關(guān)（相關(guān)系數(shù)為0.97，顯著性分別為9e-11、8e-11和7e-11），MEpink模塊與LCP和NPQ為顯著正相關(guān)（相關(guān)系數(shù)為0.87，顯著性為3e-06和4e-06）；模塊MEmagenta與E2、H32、E、AQY、LCP、LSP、NPQ、ETR、RuBP、SOD、POD和MDA 12個(gè)性狀高度相關(guān)，其中與RuBP相關(guān)性最顯著（相關(guān)系數(shù)為0.83，顯著性為2e-05）；在MEturquoise、MEgreen和MEyellow模塊中，與E1、E、Pn、Gs、Ci、Tr、Fv/Fm、Ca、Cb、Ct、SOD和MDA高度正相關(guān)，其中，MEturquoise模塊與Fv/Fm相關(guān)性最高（相關(guān)系數(shù)為0.87，顯著性為3e-06）；MEgreen模塊與E1、Ci和Ca正相關(guān)性最高（相關(guān)系數(shù)為0.9，顯著性為分別為3e-07、4e-07和5e-07）；MEyellow模塊與Fv/Fm正相關(guān)性最高（相關(guān)系數(shù)為0.89，顯著性為6e-07）；T10最顯著正相關(guān)的模塊是MEpurple（相關(guān)系數(shù)為0.91，顯著性2e-07）；T12最顯著相關(guān)的模塊是MEbrown（相關(guān)系數(shù)為0.93，顯著性為2e-08）；MEblack最顯著相關(guān)的是E2（相關(guān)系數(shù)為0.60，顯著性為0.009）。RuBP羧化酶是植物進(jìn)行光合作用過(guò)程中最重要的一個(gè)酶，而Fv/Fm即PSⅡ最大量子效率是一個(gè)重要的葉綠素?zé)晒鈪?shù)，而且本試驗(yàn)主要探索不同桑樹(shù)品種不同時(shí)間點(diǎn)的光合作用相關(guān)基因。因 此，將MEmagenta、MEgreen、MEyellow、MEpurple和MEblack共5個(gè)模塊作為重點(diǎn)研究模塊。

圖3 基因模塊與性狀相關(guān)性

2.5 基因模塊的GO和KEGG富集分析

在每一個(gè)基因模塊中聚集著表達(dá)相似的基因簇，而表達(dá)模式相似的基因在生物體內(nèi)可能執(zhí)行相似的功能和富集在相同的信號(hào)通路中。選取5個(gè)重點(diǎn)模塊進(jìn)行GO分析，進(jìn)一步探索模塊內(nèi)基因參與的生物功能和信號(hào)通路，分析結(jié)果如表2所示，5個(gè)重要模塊內(nèi)的差異基因顯著富集在碳水化合物代謝、脂質(zhì)代謝、各種離子結(jié)合、與光合作用相關(guān)的各種酶活性以及葉綠體組織、晝夜節(jié)律、葉綠體基質(zhì)、葉綠體類囊體膜和葉綠體膜等GO term中。KEEG富集分析中，如圖4和表3所示，模塊內(nèi)共有48個(gè)差異表達(dá)基因在次生代謝產(chǎn)物的生物合成、碳代謝和光合生物的固碳等通路顯著富集。根據(jù)基因的功能描述和表達(dá)情況進(jìn)一步篩選到14個(gè)與光合作用顯著相關(guān)的基因，以模塊中連通性排在前30的基因作為樞紐基因（Hub gene），再結(jié)合GO和KEGG富集分析結(jié)果，共有12個(gè)與光合作用高度相關(guān)的基因富集在光合系統(tǒng)，光合系統(tǒng)I和II等GO Term，而且主要集聚在MEblack模塊，表明MEblack模塊與光合作用高度相關(guān)。其中，MEblack模塊內(nèi)的L484_021838基因在E2-12_vs_H32-12顯著下調(diào)，與L484_000839具有共同的擬南芥同源基因，蛋白同源性為94.4%，推測(cè)該基因可能與L484_000839基因在生物體內(nèi)執(zhí)行相似的功能，也是參與光合作用的主要基因；MEpurple模塊中的L484_025354基因在E1-10_vs_E1-12和H32-10_vs_H32-12顯著上調(diào)，描述為核糖二磷酸羧化酶/加氧酶激活酶2，由RCA編碼，RCA是一種AAA+家族蛋白，通過(guò)依賴于ATP的方式促進(jìn)從Rubisco活性位點(diǎn)去除抑制位點(diǎn)而介導(dǎo)Rubisco的激活，RCA通過(guò)調(diào)節(jié)CO2固定來(lái)提高農(nóng)作物的CO2同化率、生長(zhǎng)和產(chǎn)量。

表3 基于模塊候選的光合作用相關(guān)基因

圖4 5個(gè)模塊KEGG富集通路

表2 模塊DEGs顯著富集的光合作用GO term

通過(guò)以上分析，最終篩選了3個(gè)與光合作用相關(guān)的基因，L484_000029被劃分到MEmagenta模塊，主要顯著富集在電子轉(zhuǎn)運(yùn)蛋白、PSII和PSII中的光合電子傳輸?shù)菺O Term中；L484_000836基因被劃分到MEblack模塊中，富集在核糖二磷酸羧化酶活性、固碳和葉綠體類囊體膜等GO Term中；L484_025354基因被劃分在MEpurple模塊，主要參與ATP結(jié)合過(guò)程。以上結(jié)果均表明，這3個(gè)候選基因與光合作用密切相關(guān)。

3 小結(jié)與討論

為了研究不同桑樹(shù)品種不同時(shí)間點(diǎn)光合作用相關(guān)基因的調(diào)控網(wǎng)絡(luò)，本研究利用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（WGCNA）方法，結(jié)合基因模塊與表型關(guān)聯(lián)分析，發(fā)掘和篩選不同桑樹(shù)品種光合作用差異的關(guān)鍵基因。

以E1、E2、H32這3個(gè)桑樹(shù)品種的10時(shí)和12時(shí)2個(gè)時(shí)間點(diǎn)轉(zhuǎn)錄組分析，得到5 416個(gè)差異表達(dá)基因，構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)的基因表達(dá)矩陣，共獲得10個(gè)基因模塊，篩選了MEmagenta、MEgreen、MEyellow、MEpurple和MEblack 5個(gè)光合相關(guān)基因模塊，GO和KEGG富集分析表明，其差異基因顯著富集在光合作用相關(guān)酶活性、葉綠體類囊體膜、碳代謝和光合生物碳固定等方面。

以3個(gè)桑樹(shù)品種的光合熒光參數(shù)和生理生化指標(biāo)等29個(gè)表型數(shù)據(jù)作為性狀矩陣與基因模塊進(jìn)行關(guān)聯(lián)分析，最后候選了12個(gè)與光合作用高度相關(guān)的基因，其顯著富集在光合系統(tǒng)、PSI和PSII等GO Term，并主要集聚在MEblack模塊。在MEblack模塊內(nèi)的L484_021838基因顯著下調(diào)，MEpurple模塊中的L484_025354基因在E1-10_vs_E1-12和H32-10_vs_H32-12顯著上調(diào)，由RCA編碼，是固定CO2的關(guān)鍵酶。

通過(guò)以上分析，最終候選3個(gè)與光合作用相關(guān)的基因，其中，MEmagenta模塊中顯著富集在電子轉(zhuǎn)運(yùn)蛋白、PSII和PSII中的光合電子傳輸?shù)菺O Term的L484_000029基因功能注釋為光系統(tǒng)Q（B）蛋白；MEblack模塊中富集在核糖二磷酸羧化酶活性、碳固定和葉綠體類囊體膜等GO Term中的L484_000836基因功能注釋為核糖二磷酸羧化酶；MEpurple模塊中主要參與ATP結(jié)合過(guò)程的L484_025354基因功能注釋為核糖二磷酸羧化酶/加氧酶激活酶2。L484_000029、L484_000836和L484_025354基因與桑樹(shù)光合作用密切相關(guān)。