基于加權基因共表達網絡分析探究桑樹光合相關基因的調控模式

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

（湖北省農業科學院經濟作物研究所，武漢 430064）

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

1 材料與方法

1.1 數據材料

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

1.2 基因間相關系數的計算

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

1.3 軟閾值的選擇

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

1.4 節點之間的相關性計算

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

1.5 構建分層聚類樹

根據基因間的拓撲重疊計算基因之間的相異性系數，并以相異性系數為基礎構建基因間的系統聚類樹，在系統聚類樹中，表達模式相似的基因屬于一個樹分支，本分析通過動態剪切樹法將所有差異基因根據不同的表達模式分割為不同的模塊［1］。每一個模塊內的基因具有表達模式的連通性和相似性。

1.6 基因模塊與表型的關聯

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

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

對分析獲得的模塊內基因進行GO和KEGG富集分析，探索這些模塊內的基因參與哪些生物學功能和富集在哪些通路，可挖掘哪些基因與模塊內基因具有高度的連通性，進一步探索光合作用相關基因的調控網絡。

2 結果與分析

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

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

2.2 軟閾值確定

按照無尺度網絡的標準選擇合適的軟閾值（β），軟閾值（β）與無尺度網絡評價系數（R2）的關系以及軟閾值（β）與平均連通性的關系如圖1所示。構建基因加權共表達網絡，確定軟閾值時，一般選擇無尺度網絡評價系數（R2）在0.9左右對應的軟閾值，本研究中以無尺度網絡評價系數（0.85）對應的軟閾值（16）來構建基因共表達網絡。

圖1 軟閾值確定

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

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

表1 基因模塊和對應的基因數目

圖2 基因模塊聚類

2.4 基因模塊與表型值的關聯

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

圖3 基因模塊與性狀相關性

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

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

表3 基于模塊候選的光合作用相關基因

圖4 5個模塊KEGG富集通路

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

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

3 小結與討論

為了研究不同桑樹品種不同時間點光合作用相關基因的調控網絡，本研究利用加權基因共表達網絡分析（WGCNA）方法，結合基因模塊與表型關聯分析，發掘和篩選不同桑樹品種光合作用差異的關鍵基因。

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

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

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