加權基因共表達網絡分析挖掘老鴉瓣芽莖發育關鍵基因

張軍霞，郭巧生，朱再標，徐碧霞

? 藥材與資源 ?

加權基因共表達網絡分析挖掘老鴉瓣芽莖發育關鍵基因

張軍霞，郭巧生，朱再標*，徐碧霞

南京農業大學 中藥材研究所，江蘇 南京 210095

利用老鴉瓣各部位及芽莖不同發育時期RNA-Seq及基因表達量數據，挖掘老鴉瓣芽莖發育過程關鍵基因。通過加權基因共表達網絡分析（weighted gene co-expression network analysis，WGCNA）法構建網絡，根據模塊功能富集分析及基因表達模式進行共表達模塊和核心基因的篩選。通過基因表達量相關性進一步將網絡劃分為15個模塊，將共表達模塊與老鴉瓣芽莖3個發育時期相關聯，鑒定到與芽莖發生高度相關的3個模塊，即T1 MEplum1模塊、T2 MEdarkturquoise模塊和T3 MElightcyan模塊。對3個模塊內的基因進行動態基因本體（gene ontology，GO）和京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）功能富集分析并繪制網絡關系圖，挖掘出4個與芽莖發育過程相關核心基因（、、、），涉及蛋白質合成、次生代謝產物的糖基化修飾、植物激素調節等。挖掘出的3個共表達模塊和4個核心基因有助于進一步闡明老鴉瓣芽莖發育機制。

老鴉瓣；轉錄組；加權基因共表達網絡分析；芽莖；富集分析；核心基因

老鴉瓣(Miq.) Honda為百合科老鴉瓣屬多年生草本植物[1]。其去除膜質皮和披絨毛的干燥鱗莖為藥材光慈菇，味甘、性寒，有小毒，有清熱解毒、散結消腫功效，用于治療咽喉腫痛、瘰疬、乳腺癌等多種癌癥[2-3]。由于近年來光慈菇藥用價值的不斷開發，對老鴉瓣的需求日益增長。但目前老鴉瓣以野生資源為主，人工栽培技術還在逐步完善中。老鴉瓣人工栽培中一個重要限制因子是繁殖系數。前期研究發現在1月底2月初，生長年限較短的老鴉瓣從鱗莖盤的側面形成1～4條地下莖狀結構（暫命名“芽莖”），而芽莖是影響老鴉瓣繁殖系數的主要因素[4]。

老鴉瓣芽莖兼具根、根狀莖或匍匐莖的部分特點，作為重要的無性生殖器官，對老鴉瓣快速繁育有重要作用，也可作為研究老鴉瓣物種進化和生態適應等方面的重要材料。課題組前期研究了影響老鴉瓣芽莖形成的環境條件[5]，芽莖發育過程的糖代謝物質和內源激素含量變化[6]，通過對老鴉瓣芽莖轉錄組測序，篩選到大量差異表達基因，可能參與了細胞生長、激素調節、代謝等過程并篩取部分基因進行了驗證[7]。通過對芽莖發生的miRNA測序并分析，發現5個miRNA可能在芽莖發育過程具有重要作用[8]。

隨著生命科學的發展，更多系統生物學的研究方法應用于植物器官發育的研究，如加權基因共表達網絡分析（weighted correlation network analysis，WGCNA）[9]。WGCNA被廣泛應用于醫學和植物學研究領域，核心基因與其他基因連接最緊密，更容易發現與性狀相關的生物學意義[10]，用于鑒定候選生物標記基因或治療靶點[11-12]。

本研究將WGCNA法與轉錄組測序結合，構建老鴉瓣6個組織的加權基因共表達網絡并劃分模塊，探索老鴉瓣芽莖不同發育階段基因表達模式并分析特異性模塊，以期發掘芽莖發生的核心基因，為進一步研究調控老鴉瓣芽莖發育的分子機制提供依據。

1 材料

本實驗數據來源于課題組前期對老鴉瓣不同器官及芽莖不同發育時期的轉錄組測序結果[7-8]。老鴉瓣于2018年10月栽種于南京農業大學試驗大棚，實驗材料經南京農業大學郭巧生教授鑒定為老鴉瓣(Miq.) Honda。選取其葉（L）、花冠（F）、鱗莖（B）及芽莖發育的前、中、后期（T1、T2、T3）6個時期樣品，每個重復3次，共18個樣品。

2 方法

2.1 WGCNA的構建流程

基于WGCNA軟件包和R包進行樣品聚類以及數據缺失值檢測，并基于基因表達數據變異程度過濾，選擇絕對中位差（median absolute deviation，MAD）前75%且MAD＞0.1的基因進行后續分析。通過使用WGCNA軟件包的選擇軟閾值的功能進行網絡拓撲分析[12]，選擇軟閾值β＝21（擬合指數＞0.9）進行鄰矩陣的構建，建立網絡，基于cutHeight＝0.25（相當于模塊相似的高于0.75）進行相似模塊合并，最獲得16個模塊（其中，grey模塊含未匹配到任何模塊的基因），選擇除grey模塊之外的15個模塊做分析。結合WGCNA網絡、模塊劃分結果與樣本性狀聯合分析模塊和樣本之間的關系。再通過模塊內連通性（module eigengene-based connectivity，kME）找尋每個模塊的核心基因，在共表達模塊內具有高網絡連接性的基因被定義為核心基因[13]。

2.2 特異性模塊的篩選

模塊與性狀相關性絕對值代表著模塊與性狀之間的相關性。某一性狀與某一模塊的相關性絕對值越接近1，性狀與模塊的基因功能相關性越大[14]。選擇絕對值大于0.6的模塊作為進一步分析的組織特異性模塊。

2.3 模塊基因本體（gene ontology，GO）和京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）功能富集分析

通過轉錄組數據獲得基因的GO編號，利用基迪奧云平臺網頁https://www.omicshare.com/tools 在線分析模塊的GO富集和KEGG通路富集。

2.4 核心基因的篩選與功能預測

通過分析結果中模塊kME中最高的150個基因的“點”“線”關系文件，繪制基因的Cytoscape網絡關系圖，選定位于網絡關系圖中顏色較深的基因為核心基因。經NCBI的blast對篩選得到的基因進行序列比對，根據結果，初步確定其所屬基因家族，預測其功能，進一步確定與老鴉瓣芽莖階段生長發育相關的候選基因。

3 結果與分析

3.1 WGCNA的構建

通過WGCNA方法構建老鴉瓣加權基因共表達網絡，篩選出具有較高表達量的39 931個Unigene參與網絡劃分，并進行相似模塊合并，最終獲得了15個共表達模塊（圖1），每個模塊用不同的聚類樹和不同顏色表示，其中，darkturquoise（暗綠松石色）模塊中基因數目最多，達16 101個，palevioletred（蒼紫羅藍色）模塊基因數目最少，為63個。圖2分為3部分，一是基于拓撲重疊（topological overlap dissimilarity measure，TOM）計算基因間相關聯程度繪制的基因系統聚類樹，表現每個基因的聚類關系；二是對應第1部分結果，使用動態剪枝算法將基因劃分模塊的結果；三是對第2部分結果進行優化合并，獲得的最終模塊劃分結果。

圖1 基因聚類樹和模塊劃分

3.2 組織特異性模塊的分析

通過計算模塊與形狀相關性絕對值，發現3個組織和芽莖3個發育階段共有15個組織特異性模塊（絕對值＞0.6），包括呈正相關和負相關（圖2）。其中多數在芽莖T3時期、花、鱗莖、葉片。與芽莖發育后期T3正相關的模塊有淺青色（MElightcyan，0.73）、薊色（MEthistle1，0.68）模塊；與鱗莖正相關的模塊有黃色（MEyelow，0.90）、淺黃色（MElightyellow，0.67）、棕色4（MEbrown 4，0.76）模塊；與葉片呈正相關的是棕色（MEbrown，0.98）模塊；與花呈正相關的是綠色（MEgreen，0.98）、洋紅色（MEmagenta，0.82）模塊。與T1呈負相關的是紫紅色（MEplum1，0.66）模塊；與葉片負相關的模塊有深橙色2（MEdarkorange 2，0.89）；與花負相關的模塊有夜藍色（MEmidnightblue，0.62）。模塊與組織相關性絕對值表明模塊與組織的關聯性，數值越大關聯度越高。從中篩選的組織特異性模塊參與各組織特有的生物學過程的可能性越大。為了驗證基因共表達網絡和組織特異性模塊兩者結果的可靠性，同時為了研究芽莖發育過程，對T1、T3的組織特異性模塊MEplum1、MElightcyan與T2相關系數高的模塊MEdarkturquoise做進一步分析。

3.3 芽莖的組織特異性表達

通過對芽莖發育T1、T2、T3時期的MEplum1、MEdarkturquoise、MElightcyan模塊做GO富集分析（圖3-A～C）。3個模塊相比較，在生物學過程（biological progress，BP）、細胞組分（cellur component，CC）、分子功能（molecular function，MF）3個方面均發揮主要作用。其中的生物過程包括了細胞過程、代謝過程、固著生物過程，分子功能包括結合、催化活性，細胞組分包括細胞、細胞組成、細胞器等均占有較大比重。MEdarkturquoise模塊基因與MEplum1模塊中的基因相比，生物學過程中的細胞成分組織或生物發生、生長發育過程、生物節律過程所占比率有所上升，MElightcyan模塊的定位、細胞成分組織或生物發生及生長發育過程的比重增加；在分子功能方面，MEdarkturquoise模塊的基因在抗氧化活性、核苷酸結合轉錄因子、信號轉換器所占比率增加，其中抗氧化活性提升的比例較高，MElightcyan模塊則是催化活性、轉運體活性的比率有所上升；在MEdarkturquoise模塊基因的細胞組分方面，胞外區的比例明顯增加，MElightcyan模塊基因的細胞膜、細胞膜組成這2項較MEplum1模塊基因的比例有增加。

圖2 老鴉瓣6個組織特異性模塊

A-T1 MEplum1模塊 B-T2 MEdarkturquoise模塊 C-T3 MElightcyan模塊

再以芽莖T1模塊內的基因為對照，著重對芽莖的T2模塊進行KEGG富集分析（圖4-A～B），結果發現，與MEplum1模塊相比，MEdarkturquoise模塊的基因在代謝過程的比例略有增加，主要是在氨基酸代謝、脂類代謝、能量代謝、次生代謝及生物合成、糖類生物合成與代謝；基因信息處理過程的翻譯等方面發揮著主要作用。而與MEdarkturquoise模塊相比，MEplum1模塊基因無明顯差異。

3.4 模塊核心基因的篩選與功能預測

為縮小篩選范圍，選取模塊kME較高的150個基因繪制網絡關系圖（圖5-A～C），根據網絡圖示信息，將T2時期kME較高同時權重值大的73個基因分布在5個類別下，分別含有45、11、7、3、7條基因，將權重值最高的第一個類別中的45個基因繪制熱圖分析（圖6），并通過BLASTP比對，發現15條序列能確切比對到信息，其中有6條比對到核糖體蛋白，找到2條已經試驗驗證過的核糖體蛋白序列、，同時基于本課題組前期對老鴉瓣芽莖的生理生化分析，找到注釋為糖基轉移酶的和注釋為細胞色素P450的的4條基因序列。

A-T1 MEplum1模塊 B-T2 MEdarkturquoise模塊

其中，T2 MEdarkturquoise模塊中的、比對到核糖體蛋白XP_013683027.1、KAE8711130.1，比對到細胞色素P450 704C1 OAY69779.1。為糖基轉移酶RWR91480.1。

4 討論

在前期研究中發現，通過轉錄組測序，用GO、COG、KEGG數據庫對芽莖形成發育的差異基因進行注釋，功能注釋顯示這些差異基因主要涉及物質及能量代謝、激素信號、細胞生長、轉錄調控等諸多生理生化過程[7]。蔗糖等可溶性糖及淀粉在眾多酶，包括蔗糖磷酸合酶（sucrose phosphate synthase，SPS）、淀粉酶（amylase，AMY）、蔗糖合酶（sucrose synthase，SS）、腺苷二磷酸葡萄糖焦磷酸化酶（adenosine diphosphoglucose pyrophosphorylase，AGPase）、可溶性淀粉合成酶（soluble starch synthase，SSS）、顆粒結合型淀粉合成酶（granule-bound starch synthase，GBSS）的共同作用下大量積累，為芽莖發生時的細胞分裂和生長提供物質和能量需要；赤霉素（gibberellins，GA）、玉米素核苷（zeatin riboside，ZR）、吲哚-3-乙酸（indole-3-acetic acid，IAA）及脫落酸（abscisic acid，ABA）協同作用，相互制約平衡，共同作用于芽莖的發生[6]。

A-T1 MEplum1模塊基因網絡 B-T2 MEdarkturquoise模塊基因網絡 C-T3 MElightcyan模塊基因網絡

箭頭代表篩選的4個基因

WGCNA能特異地篩選出與目標性狀相關的基因，通過進行模塊化分類，得到相關度高的共表達模塊，已經被證明是一種高效的數據挖掘手段[15]。Hollender等[16]利用WGCNA法發現了草莓中與花托組織特異性有關的模塊，并挖掘出7個相關核心基因。王思瑤等[14]通過WGCNA法構建網絡，對特異性模塊進行分析，發掘出36個潛在的與偃松種子萌發相關的核心基因。鞠正等[17]利用番茄組織的RNA-Seq通過WGCNA法構建網絡分析，在網絡和組織特異性模塊發現了與果實成熟相關的轉錄因子。

本研究利用老鴉瓣芽莖與其他組織轉錄組數據，通過WGCNA法構建了15個基因共表達模塊，分析3個芽莖發育相關的組織模塊，分別為T1 MEplum1、T2 MEdarkturquoise、T3 MElightcyan模塊，同時對模塊進行GO、KEGG富集分析并繪制模塊內基因的網絡分析圖，于T2 MEdarkturquoise模塊篩選出4個核心基因，其可能在芽莖發育過程起到重要作用。

篩選出的4個核心基因中，、比對到核糖體蛋白。核糖體蛋白多存在于快速增殖、分泌或功能旺盛的細胞中，并參與植物發育過程[18]。除了與核糖體RNA組裝形成核糖體參與細胞內蛋白質合成外，還與DNA復制、RNA加工、細胞增殖、生長發育等調控有關[19]。比對到糖基轉移酶，以家族形式存在，在次生代謝產物的糖基化修飾、內源或外源物質的解毒、機體防御反應、植物激素調節等方面發揮著重要作用[20]。擬南芥中新型糖基轉移酶基因能對生長素的主要前體吲哚丙酮酸（IPyA）進行特異糖基化修飾，對活性生長素本身卻沒有活性。通過分析激素水平，發現UGT76F1催化的IPyA糖基化能夠調控IAA生物合成的代謝，參與生長素穩態調節[21]。同時，糖基轉移酶家族在糖類代謝過程中參與糖類轉移，其中經催化活化的木糖基轉移酶是多糖生物合成中的關鍵酶[22]。比對到細胞色素P450，是一種氧化酶家族，廣泛參與植物代謝過程，具有催化活性，主要涉及植物激素、生物堿、萜類等代謝產物的合成[23]。糖基轉移酶基因為金鐵索的三萜皂苷類化合物合成修飾環節的關鍵基因，與P450單加氧酶基因對三萜皂苷類化合物的共同母核——β-香樹素進行修飾，最終形成三萜皂苷[24]。

Miao等[6]研究發現老鴉瓣的芽莖由發生T1時期至膨大T3時期到成為新鱗莖的過程中，淀粉和蛋白質的含量持續上升，可溶性糖被大量消耗，含量呈下降趨勢。同時，生長素的含量在T1時期達到頂峰，于T2時期維持較高水平，T3時期顯著下降。由此推測在芽莖發生過程中可能有大量的核糖體蛋白參與蛋白質的合成，且高含量生長素促進芽莖膨大，在此過程中核糖體蛋白的活性最高，通過能量合成大量芽莖發育膨大所需的蛋白質，為此提供了物質基礎。

老鴉瓣芽莖膨大后形成的鱗莖入藥，其中含有皂苷類、多糖類成分，推測、基因可能為老鴉瓣的次生代謝研究提供物質基礎。

本研究重點關注了與老鴉瓣芽莖發育相關模塊，通過WGCNA法的模塊組織相關性分析挖掘出潛在的與芽莖發育相關的4個核心基因，為進一步探究核心基因與芽莖發育的關系奠定了堅實基礎，后續將對這些基因的功能進行驗證。

利益沖突 所有作者均聲明不存在利益沖突

[1] 譚敦炎, 張震, 李新蓉, 等. 老鴉瓣屬(百合科)的恢復: 以形態性狀的分支分析為依據 [J]. 植物分類學報, 2005, 43(3): 262-270.

[2] 國家中醫藥管理局《中華本草》編委會. 中華本草: 藏藥卷 [M]. 上海: 上海科學技術出版社, 2002: 176-177.

[3] 陳彪, 焦淑萍, 尹榮, 等. 6種吉林抗癌中藥清除羥自由基及其抗DNA損傷體外實驗研究 [J]. 第三軍醫大學學報, 2004, 26(1): 88-89.

[4] 楊小花, 繆媛媛, 郭巧生, 等. 不同等級老鴉瓣生長繁殖特征研究 [J]. 中草藥, 2015, 46(24): 3746-3750.

[5] 繆媛媛. 老鴉瓣芽莖發育特征及形成機制研究 [D]. 南京: 南京農業大學, 2015.

[6] Miao Y Y, Zhu Z B, Guo Q S,. Dynamic changes in carbohydrate metabolism and endogenous hormones duringstolon development into a new bulb [J]., 2016, 59(2): 121-132.

[7] Miao Y Y, Zhu Z B, Guo Q S,. Transcriptome analysis of differentially expressed genes provides insight into stolon formation in[J]., 2016, 7: 409.

[8] Zhu Z B, Miao Y Y, Guo Q S,. Identification of miRNAs involved in stolon formation inby high-throughput sequencing [J]., 2016, 7: 852.

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

[10] 劉源, 隋正紅, 劉昊昕, 等. 加權基因共表達網絡探究鏈狀亞歷山大藻爆發性生長的分子機制 [J]. 中國海洋大學學報: 自然科學版, 2019, 49(9): 66-76.

[11] Yin L, Cai Z H, Zhu B A,. Identification of key pathways and genes in the dynamic progression of HCC based on WGCNA [J].(), 2018, 9(2): E92.

[12] Gao C, Ju Z, Li S,. Deciphering ascorbic acid regulatory pathways in ripening tomato fruit using a weighted gene correlation network analysis approach [J]., 2013, 55(11): 1080-1091.

[13] Fuller T F, Ghazalpour A, Aten J E,. Weighted gene coexpression network analysis strategies applied to mouse weight [J]., 2007, 18(6/7): 463-472.

[14] 王思瑤, 叢日征, 閆曉娜, 等. 利用加權基因共表達網絡分析 (WGCNA) 的方法挖掘偃松種子萌發過程關鍵基因 [J]. 溫帶林業研究, 2019, 2(1): 39-46.

[15] Zhao W, Langfelder P, Fuller T,. Weighted gene coexpression network analysis: State of the art [J]., 2010, 20(2): 281-300.

[16] Hollender C A, Kang C Y, Darwish O,. Floral transcriptomes in woodland strawberry uncover developing receptacle and anther gene networks [J]., 2014, 165(3): 1062-1075.

[17] 鞠正, 曹東艷, 梁巖, 等. 利用加權基因共表達網絡分析 (WGCNA) 的方法挖掘番茄果實成熟相關的轉錄因子 [J]. 中國食品學報, 2018, 18(6): 240-248.

[18] 廖昌敏, 張詠祀, 劉小紅, 等. 水杉40S核糖體蛋白S8基因的克隆及生物信息學分析 [J]. 分子植物育種, 2020, 18(21): 7008-7014.

[19] Labriet A, Lévesque é, Cecchin E,. Germline variability and tumor expression level of ribosomal protein gene RPL28 are associated with survival of metastatic colorectal cancer patients [J]., 2019, 9(1): 13008.

[20] 馬風偉, 鄧青芳, 陳海江, 等. 植物UDP-糖基轉移酶結構及活性研究進展 [J]. 貴陽學院學報: 自然科學版, 2019, 14(3): 73-82.

[21] Chen L, Huang X X, Zhao S M,. IPyA glucosylation mediates light and temperature signaling to regulate auxin-dependent hypocotyl elongation in[J]., 2020, 117(12): 6910-6917.

[22] Han W T, Fan X, Teng L H,. Identification, classification, and evolution of putative xylosyltransferases from algae [J]., 2019, 256(4): 1119-1132.

[23] 朱靈英, 郭娟, 張愛麗, 等. 參與植物三萜生物合成的細胞色素P450酶研究進展 [J]. 中草藥, 2019, 50(22): 5597-5610.

[24] 江舟. 金鐵鎖糖基轉移酶的克隆與生物信息學分析 [D]. 昆明: 云南中醫學院, 2016.

Key Genes involving in stolon development ofwith weighted gene co-expression network analysis

ZHANG Jun-xia, GUO Qiao-sheng, ZHU Zai-biao, XU Bi-xia

Institute of Chinese Medicinal Materials, Nanjing Agricultural University, Nanjing 210095, China

To detect the hub genes of the stolon development process inby using the RNA-Seq and gene expression data from various parts ofand bud stems at different developmental periods.The network was constructed by weighted gene co-expression network analysis (WGCNA) method, and co-expression modules and core genes were screened based on module functional enrichment analysis and gene expression patterns.Fifteen modules were divided through the correlation of gene expression, the co-expression module was associated with the three developmental stages ofstolon, and three modules were identified to be highly related to stolon development, namely T1 MEplum1 module, T2 MEdarkturquoise module and T3 MElightcyan module, respectively. The analysis of gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) function enrichment of the three modules and network relationship diagrams found that four hub genes related to the stolon development process, namely,,,, which were involved in protein synthesis, glycosylation modification of secondary metabolites, plant hormone regulation etc.The identification of three co-expression modules and four hub genes mined will be helpful for further elucidating the mechanism aboutstolon development.

(Miq.) Honda; transcriptome; weighted gene co-expression network analysis; stolon; enrichment analysis; hub genes

R282.12

A

0253 - 2670(2023)04 - 1228 - 08

10.7501/j.issn.0253-2670.2023.04.023

2022-08-20

國家自然科學基金資助項目（81773834）；中央高校基本科研業務費項目（KYZZ2022004）

張軍霞，碩士研究生，研究方向為藥用植（動）物種植（養殖）理論與技術。E-mail: 2017104127@njau.edu.cn

朱再標，教授，研究方向為藥用植（動）物種植（養殖）理論與技術。Tel: (025)84395980 E-mail: zhuzaibiao@njau.edu.cn

[責任編輯 時圣明]