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

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

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

1 材料與方法

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

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

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

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

1.3 軟閾值的選擇

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

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

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

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

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

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

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

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

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

2 結(jié)果與分析

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

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

2.2 軟閾值確定

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

圖1 軟閾值確定

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

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

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

圖2 基因模塊聚類

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

為了進一步探索基因模塊與不同性狀的關(guān)系，使用每個基因模塊的模塊特征值（ME）與本研究中的29個性狀進行關(guān)聯(lián)分析，來確定某一個模塊基因是否與某個特定性狀高度相關(guān)聯(lián)，本試驗重點關(guān)注與性狀正相關(guān)的模塊基因。如圖3所示，發(fā)現(xiàn)MEblue和MEpink模塊分別與H32、AQY、LCP、CCP、CSP、NPQ、ETR和POD 8個性狀高度正相關(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個性狀高度相關(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羧化酶是植物進行光合作用過程中最重要的一個酶，而Fv/Fm即PSⅡ最大量子效率是一個重要的葉綠素熒光參數(shù)，而且本試驗主要探索不同桑樹品種不同時間點的光合作用相關(guān)基因。因 此，將MEmagenta、MEgreen、MEyellow、MEpurple和MEblack共5個模塊作為重點研究模塊。

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

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

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

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

圖4 5個模塊KEGG富集通路

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

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

3 小結(jié)與討論

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

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

以3個桑樹品種的光合熒光參數(shù)和生理生化指標等29個表型數(shù)據(jù)作為性狀矩陣與基因模塊進行關(guān)聯(lián)分析，最后候選了12個與光合作用高度相關(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)鍵酶。

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