張軍霞，郭巧生，朱再標(biāo)，徐碧霞

? 藥材與資源 ?

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析挖掘老鴉瓣芽莖發(fā)育關(guān)鍵基因

張軍霞，郭巧生，朱再標(biāo)*，徐碧霞

南京農(nóng)業(yè)大學(xué) 中藥材研究所，江蘇 南京 210095

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

老鴉瓣；轉(zhuǎn)錄組；加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析；芽莖；富集分析；核心基因

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

老鴉瓣芽莖兼具根、根狀莖或匍匐莖的部分特點(diǎn)，作為重要的無(wú)性生殖器官，對(duì)老鴉瓣快速繁育有重要作用，也可作為研究老鴉瓣物種進(jìn)化和生態(tài)適應(yīng)等方面的重要材料。課題組前期研究了影響老鴉瓣芽莖形成的環(huán)境條件[5]，芽莖發(fā)育過(guò)程的糖代謝物質(zhì)和內(nèi)源激素含量變化[6]，通過(guò)對(duì)老鴉瓣芽莖轉(zhuǎn)錄組測(cè)序，篩選到大量差異表達(dá)基因，可能參與了細(xì)胞生長(zhǎng)、激素調(diào)節(jié)、代謝等過(guò)程并篩取部分基因進(jìn)行了驗(yàn)證[7]。通過(guò)對(duì)芽莖發(fā)生的miRNA測(cè)序并分析，發(fā)現(xiàn)5個(gè)miRNA可能在芽莖發(fā)育過(guò)程具有重要作用[8]。

隨著生命科學(xué)的發(fā)展，更多系統(tǒng)生物學(xué)的研究方法應(yīng)用于植物器官發(fā)育的研究，如加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（weighted correlation network analysis，WGCNA）[9]。WGCNA被廣泛應(yīng)用于醫(yī)學(xué)和植物學(xué)研究領(lǐng)域，核心基因與其他基因連接最緊密，更容易發(fā)現(xiàn)與性狀相關(guān)的生物學(xué)意義[10]，用于鑒定候選生物標(biāo)記基因或治療靶點(diǎn)[11-12]。

本研究將WGCNA法與轉(zhuǎn)錄組測(cè)序結(jié)合，構(gòu)建老鴉瓣6個(gè)組織的加權(quán)基因共表達(dá)網(wǎng)絡(luò)并劃分模塊，探索老鴉瓣芽莖不同發(fā)育階段基因表達(dá)模式并分析特異性模塊，以期發(fā)掘芽莖發(fā)生的核心基因，為進(jìn)一步研究調(diào)控老鴉瓣芽莖發(fā)育的分子機(jī)制提供依據(jù)。

1 材料

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

2 方法

2.1 WGCNA的構(gòu)建流程

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

2.2 特異性模塊的篩選

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

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

通過(guò)轉(zhuǎn)錄組數(shù)據(jù)獲得基因的GO編號(hào)，利用基迪奧云平臺(tái)網(wǎng)頁(yè)https://www.omicshare.com/tools 在線分析模塊的GO富集和KEGG通路富集。

2.4 核心基因的篩選與功能預(yù)測(cè)

通過(guò)分析結(jié)果中模塊kME中最高的150個(gè)基因的“點(diǎn)”“線”關(guān)系文件，繪制基因的Cytoscape網(wǎng)絡(luò)關(guān)系圖，選定位于網(wǎng)絡(luò)關(guān)系圖中顏色較深的基因?yàn)楹诵幕?。?jīng)NCBI的blast對(duì)篩選得到的基因進(jìn)行序列比對(duì)，根據(jù)結(jié)果，初步確定其所屬基因家族，預(yù)測(cè)其功能，進(jìn)一步確定與老鴉瓣芽莖階段生長(zhǎng)發(fā)育相關(guān)的候選基因。

3 結(jié)果與分析

3.1 WGCNA的構(gòu)建

通過(guò)WGCNA方法構(gòu)建老鴉瓣加權(quán)基因共表達(dá)網(wǎng)絡(luò)，篩選出具有較高表達(dá)量的39 931個(gè)Unigene參與網(wǎng)絡(luò)劃分，并進(jìn)行相似模塊合并，最終獲得了15個(gè)共表達(dá)模塊（圖1），每個(gè)模塊用不同的聚類樹和不同顏色表示，其中，darkturquoise（暗綠松石色）模塊中基因數(shù)目最多，達(dá)16 101個(gè)，palevioletred（蒼紫羅藍(lán)色）模塊基因數(shù)目最少，為63個(gè)。圖2分為3部分，一是基于拓?fù)渲丿B（topological overlap dissimilarity measure，TOM）計(jì)算基因間相關(guān)聯(lián)程度繪制的基因系統(tǒng)聚類樹，表現(xiàn)每個(gè)基因的聚類關(guān)系；二是對(duì)應(yīng)第1部分結(jié)果，使用動(dòng)態(tài)剪枝算法將基因劃分模塊的結(jié)果；三是對(duì)第2部分結(jié)果進(jìn)行優(yōu)化合并，獲得的最終模塊劃分結(jié)果。

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

3.2 組織特異性模塊的分析

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

3.3 芽莖的組織特異性表達(dá)

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

圖2 老鴉瓣6個(gè)組織特異性模塊

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

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

3.4 模塊核心基因的篩選與功能預(yù)測(cè)

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

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

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

4 討論

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

A-T1 MEplum1模塊基因網(wǎng)絡(luò) B-T2 MEdarkturquoise模塊基因網(wǎng)絡(luò) C-T3 MElightcyan模塊基因網(wǎng)絡(luò)

箭頭代表篩選的4個(gè)基因

WGCNA能特異地篩選出與目標(biāo)性狀相關(guān)的基因，通過(guò)進(jìn)行模塊化分類，得到相關(guān)度高的共表達(dá)模塊，已經(jīng)被證明是一種高效的數(shù)據(jù)挖掘手段[15]。Hollender等[16]利用WGCNA法發(fā)現(xiàn)了草莓中與花托組織特異性有關(guān)的模塊，并挖掘出7個(gè)相關(guān)核心基因。王思瑤等[14]通過(guò)WGCNA法構(gòu)建網(wǎng)絡(luò)，對(duì)特異性模塊進(jìn)行分析，發(fā)掘出36個(gè)潛在的與偃松種子萌發(fā)相關(guān)的核心基因。鞠正等[17]利用番茄組織的RNA-Seq通過(guò)WGCNA法構(gòu)建網(wǎng)絡(luò)分析，在網(wǎng)絡(luò)和組織特異性模塊發(fā)現(xiàn)了與果實(shí)成熟相關(guān)的轉(zhuǎn)錄因子。

本研究利用老鴉瓣芽莖與其他組織轉(zhuǎn)錄組數(shù)據(jù)，通過(guò)WGCNA法構(gòu)建了15個(gè)基因共表達(dá)模塊，分析3個(gè)芽莖發(fā)育相關(guān)的組織模塊，分別為T1 MEplum1、T2 MEdarkturquoise、T3 MElightcyan模塊，同時(shí)對(duì)模塊進(jìn)行GO、KEGG富集分析并繪制模塊內(nèi)基因的網(wǎng)絡(luò)分析圖，于T2 MEdarkturquoise模塊篩選出4個(gè)核心基因，其可能在芽莖發(fā)育過(guò)程起到重要作用。

篩選出的4個(gè)核心基因中，、比對(duì)到核糖體蛋白。核糖體蛋白多存在于快速增殖、分泌或功能旺盛的細(xì)胞中，并參與植物發(fā)育過(guò)程[18]。除了與核糖體RNA組裝形成核糖體參與細(xì)胞內(nèi)蛋白質(zhì)合成外，還與DNA復(fù)制、RNA加工、細(xì)胞增殖、生長(zhǎng)發(fā)育等調(diào)控有關(guān)[19]。比對(duì)到糖基轉(zhuǎn)移酶，以家族形式存在，在次生代謝產(chǎn)物的糖基化修飾、內(nèi)源或外源物質(zhì)的解毒、機(jī)體防御反應(yīng)、植物激素調(diào)節(jié)等方面發(fā)揮著重要作用[20]。擬南芥中新型糖基轉(zhuǎn)移酶基因能對(duì)生長(zhǎng)素的主要前體吲哚丙酮酸（IPyA）進(jìn)行特異糖基化修飾，對(duì)活性生長(zhǎng)素本身卻沒有活性。通過(guò)分析激素水平，發(fā)現(xiàn)UGT76F1催化的IPyA糖基化能夠調(diào)控IAA生物合成的代謝，參與生長(zhǎng)素穩(wěn)態(tài)調(diào)節(jié)[21]。同時(shí)，糖基轉(zhuǎn)移酶家族在糖類代謝過(guò)程中參與糖類轉(zhuǎn)移，其中經(jīng)催化活化的木糖基轉(zhuǎn)移酶是多糖生物合成中的關(guān)鍵酶[22]。比對(duì)到細(xì)胞色素P450，是一種氧化酶家族，廣泛參與植物代謝過(guò)程，具有催化活性，主要涉及植物激素、生物堿、萜類等代謝產(chǎn)物的合成[23]。糖基轉(zhuǎn)移酶基因?yàn)榻痂F索的三萜皂苷類化合物合成修飾環(huán)節(jié)的關(guān)鍵基因，與P450單加氧酶基因?qū)θ圃碥疹惢衔锏墓餐负恕?香樹素進(jìn)行修飾，最終形成三萜皂苷[24]。

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

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

本研究重點(diǎn)關(guān)注了與老鴉瓣芽莖發(fā)育相關(guān)模塊，通過(guò)WGCNA法的模塊組織相關(guān)性分析挖掘出潛在的與芽莖發(fā)育相關(guān)的4個(gè)核心基因，為進(jìn)一步探究核心基因與芽莖發(fā)育的關(guān)系奠定了堅(jiān)實(shí)基礎(chǔ)，后續(xù)將對(duì)這些基因的功能進(jìn)行驗(yàn)證。

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

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

[2] 國(guó)家中醫(yī)藥管理局《中華本草》編委會(huì). 中華本草: 藏藥卷 [M]. 上海: 上海科學(xué)技術(shù)出版社, 2002: 176-177.

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

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

[5] 繆媛媛. 老鴉瓣芽莖發(fā)育特征及形成機(jī)制研究 [D]. 南京: 南京農(nóng)業(yè)大學(xué), 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] 劉源, 隋正紅, 劉昊昕, 等. 加權(quán)基因共表達(dá)網(wǎng)絡(luò)探究鏈狀亞歷山大藻爆發(fā)性生長(zhǎng)的分子機(jī)制 [J]. 中國(guó)海洋大學(xué)學(xué)報(bào): 自然科學(xué)版, 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] 王思瑤, 叢日征, 閆曉娜, 等. 利用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析 (WGCNA) 的方法挖掘偃松種子萌發(fā)過(guò)程關(guān)鍵基因 [J]. 溫帶林業(yè)研究, 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] 鞠正, 曹東艷, 梁巖, 等. 利用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析 (WGCNA) 的方法挖掘番茄果實(shí)成熟相關(guān)的轉(zhuǎn)錄因子 [J]. 中國(guó)食品學(xué)報(bào), 2018, 18(6): 240-248.

[18] 廖昌敏, 張?jiān)侅? 劉小紅, 等. 水杉40S核糖體蛋白S8基因的克隆及生物信息學(xué)分析 [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] 馬風(fēng)偉, 鄧青芳, 陳海江, 等. 植物UDP-糖基轉(zhuǎn)移酶結(jié)構(gòu)及活性研究進(jìn)展 [J]. 貴陽(yáng)學(xué)院學(xué)報(bào): 自然科學(xué)版, 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] 朱靈英, 郭娟, 張愛麗, 等. 參與植物三萜生物合成的細(xì)胞色素P450酶研究進(jìn)展 [J]. 中草藥, 2019, 50(22): 5597-5610.

[24] 江舟. 金鐵鎖糖基轉(zhuǎn)移酶的克隆與生物信息學(xué)分析 [D]. 昆明: 云南中醫(yī)學(xué)院, 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

國(guó)家自然科學(xué)基金資助項(xiàng)目（81773834）；中央高?；究蒲袠I(yè)務(wù)費(fèi)項(xiàng)目（KYZZ2022004）

張軍霞，碩士研究生，研究方向?yàn)樗幱弥玻▌?dòng)）物種植（養(yǎng)殖）理論與技術(shù)。E-mail: 2017104127@njau.edu.cn

朱再標(biāo)，教授，研究方向?yàn)樗幱弥玻▌?dòng)）物種植（養(yǎng)殖）理論與技術(shù)。Tel: (025)84395980 E-mail: zhuzaibiao@njau.edu.cn

[責(zé)任編輯 時(shí)圣明]