利用WGCNA挖掘黃精屬植物多糖合成通路關(guān)鍵基因*

汪 蓉,孟 然,陳雪梅,孟 盈

(吉首大學(xué)生物資源與環(huán)境科學(xué)學(xué)院,湖南 吉首 416000)

黃精屬(PolygonatumMill.)植物是多年生草本植物,隸屬于天門冬科洛林亞科黃精族[1],該屬在我國廣泛分布,可藥食兩用,常見于南北方不同的中醫(yī)藥派系及相關(guān)食譜中[2-4].在最新版《中國藥典》中,中藥黃精具體是指黃精(PolygonatumsibiricumDelar. ex Redoute)、多花黃精(PolygonatumcyrtonemaHua)和滇黃精(PolygonatumkingianumCollett &Hemsl.)3個種[4].其中黃精主要分布在我國北方,滇黃精主要分布于我國西南地區(qū),多花黃精廣泛分布在國內(nèi)多個地區(qū)且已有人工培育居群[5-6].現(xiàn)代藥理學(xué)及市場應(yīng)用中,黃精及相關(guān)提取物主要用于抗腫瘤、抗氧化、抗衰老等領(lǐng)域,并有較好的療效.黃精的主要活性成分有多糖、皂苷、黃酮等化學(xué)成分,其中糖類較多,以黃精多糖和低聚糖為主[7-10].另外,黃精屬中的玉竹(Polygonatumodoratum(Mill.) Druce)也是重要的中藥資源.最新研究表明,玉竹的主要活性成分具有降血糖、抗炎、調(diào)節(jié)免疫和抗腫瘤作用[11].黃精屬植物具有重要的醫(yī)藥價值,當(dāng)前研究方向主要集中在其藥理活性、提取工藝和人工培育技術(shù)等方面[12-14].在后基因組時代,黃精屬植物仍然缺乏高質(zhì)量的全基因組數(shù)據(jù)作為參考.目前,在有效成分的調(diào)控機(jī)制方面,部分研究人員基于轉(zhuǎn)錄組數(shù)據(jù)及相關(guān)算法,發(fā)現(xiàn)了黃精屬植物中多糖有效成分的多個調(diào)控基因,這為黃精屬植物種質(zhì)資源的保護(hù)與合理開發(fā)利用奠定了研究基礎(chǔ)[14-16].隨著測序技術(shù)的提升,這個方向的研究也在不停地推進(jìn)與深化.

生物信息分析技術(shù)的更新讓轉(zhuǎn)錄組學(xué)的分析方案不再局限于差異的查找與定位,越來越多的研究強調(diào)通過整合多個因素解決更加具體的科學(xué)問題[17-19],以便在基因數(shù)據(jù)中找到一些高度關(guān)聯(lián)的基因表達(dá)網(wǎng)絡(luò)來解釋復(fù)雜的生物學(xué)現(xiàn)象.近年來,加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(Weighted Gene Co-expression Network Analysis,WGCNA)被廣泛應(yīng)用于多個研究領(lǐng)域,其可通過無監(jiān)督的機(jī)器學(xué)習(xí)內(nèi)核,將不同樣本的表達(dá)量數(shù)據(jù)進(jìn)行表達(dá)網(wǎng)絡(luò)計算,并根據(jù)基因表達(dá)量的相關(guān)性劃分模塊,分析這些相對獨立的模塊與分組之間的關(guān)系來定位相關(guān)靶基因.本研究基于黃精屬中多花黃精和玉竹的根狀莖及葉片的轉(zhuǎn)錄組數(shù)據(jù),采用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析方法,結(jié)合KEGG數(shù)據(jù)庫中多糖合成相關(guān)通路信息,計算核心基因在多花黃精、玉竹的葉片及根狀莖之間表達(dá)模式的相關(guān)性,以獲取影響黃精屬植物多糖類有效成分合成的關(guān)鍵基因.

1 材料和方法

1.1 數(shù)據(jù)獲取與分析

本研究以黃精屬下多花黃精和玉竹2個物種根狀莖和葉片的轉(zhuǎn)錄組數(shù)據(jù)開展后續(xù)分析.在NCBI中篩選SRP166090和SRP170979隊列的相關(guān)樣本,共計14個.其中SRP166090是多花黃精的轉(zhuǎn)錄組測序結(jié)果,基于Illumina HiSeq 4000平臺完成建庫測序,選取其中的根狀莖(SRR8074960)和葉片(SRR8074961)樣品參與分析[20];SRP170979是玉竹的轉(zhuǎn)錄組測序結(jié)果,同樣基于Illumina HiSeq 4000平臺完成轉(zhuǎn)錄組建庫測序,篩減其中的非目標(biāo)器官數(shù)據(jù),共保留12個樣品[21].使用Aspera軟件下載所有的fastq文件,并通過FastQC軟件對原始數(shù)據(jù)進(jìn)行質(zhì)量評估,然后使用Fastp軟件過濾掉低質(zhì)量的reads和接頭序列.

使用Trinity軟件完成轉(zhuǎn)錄組組裝,然后執(zhí)行CD-HIT-EST腳本過濾組裝結(jié)果中的冗余轉(zhuǎn)錄本.基于liliopsida_odb10數(shù)據(jù)庫,使用BUSCO軟件完成組裝結(jié)果的質(zhì)量檢測,并通過BLAST2GO軟件對組裝結(jié)果進(jìn)行注釋.用align_and_estimate_abundance.pl軟件調(diào)用RSEM分析不同樣品的表達(dá)定量,同時用FPKM值來衡量不同測序深度造成的表達(dá)量差異.

1.2 構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)

使用WGCNA 1.69(R 4.1)軟件包進(jìn)行加權(quán)基因共表達(dá)網(wǎng)絡(luò)構(gòu)建,以均一化后的樣本基因表達(dá)矩陣及樣品來源組分信息作為輸入數(shù)據(jù),利用表達(dá)量最高的8 000個基因構(gòu)建共表達(dá)網(wǎng)絡(luò).為進(jìn)行WGCNA計算中的軟閾值篩選,通過pickSoftThreshold函數(shù)來比較不同軟閾值參數(shù)下相關(guān)性的差異,選擇“Power=12”計算最終的無尺差異鄰接矩陣.基于該矩陣數(shù)據(jù)進(jìn)一步計算拓?fù)溆胁町惖木仃?利用動態(tài)剪切算法將基因進(jìn)行聚類分析.通過比較不同模塊與分組信息之間的相關(guān)性,挑選出最為相關(guān)的模塊,利用Module Membership及Gene Significance參數(shù)過濾模塊中的基因.

1.3 多糖合成通路相關(guān)性

通過KEGG數(shù)據(jù)庫的注釋數(shù)據(jù),篩選WGCNA分析結(jié)果中的核心基因,提取其中與多糖合成通路相關(guān)的基因.使用R中的limma軟件包對所有基因進(jìn)行篩選,挑選分組間存在顯著表達(dá)差異的基因.差異基因的挑選標(biāo)準(zhǔn)設(shè)定為P_value<0.05,并且logFC閾值服從整體數(shù)據(jù)分布的95%分位數(shù).將差異基因的表達(dá)量矩陣按照分組信息拆分,使用R中的cor ()函數(shù)計算表達(dá)矩陣的相關(guān)性,執(zhí)行g(shù)gpubr軟件包中的ggscatter函數(shù)來查驗這些基因在根狀莖數(shù)據(jù)與葉片數(shù)據(jù)間的皮爾森相關(guān)性系數(shù).

2 結(jié)果

2.1 組裝與注釋

黃精屬植物作為傳統(tǒng)中藥,具有重要的藥用價值,但是由于其基因組龐大,全基因組測序數(shù)據(jù)組裝困難,目前仍然缺乏高質(zhì)量的全基因組數(shù)據(jù).在本研究中,筆者將所有樣品的原始數(shù)據(jù)合并后使用Trinity軟件進(jìn)行無參組裝.組裝得到的結(jié)果共計207 Mb,其中共獲得了394 149條轉(zhuǎn)錄本,GC占比43.25%,N50值664 bp.BUSCO分析共匹配到數(shù)據(jù)庫中3 236個基因,其中單拷貝基因1 853個(占比57.26%),多拷貝基因969個(占比29.94%),還有277個僅部分比對上的基因,以及137個未比對到數(shù)據(jù)庫的基因.筆者在多個數(shù)據(jù)庫中進(jìn)行比對注釋,大量轉(zhuǎn)錄本被比對到多糖合成相關(guān)通路的關(guān)鍵基因序列上,如糖酵解/糖異生(ko00010)、淀粉和蔗糖代謝(ko00500)、氨基糖和核苷酸糖代謝(ko00520)等通路.表1列出了多糖合成的相關(guān)通路及相應(yīng)的基因數(shù)量.此外,本研究還進(jìn)行了NR,eggNOG,KEGG,GO數(shù)據(jù)庫的注釋,注釋結(jié)果如圖1所示.

表1 多糖合成通路

圖1 非冗余組裝結(jié)果注釋信息Fig. 1 Annotation Information of Unigenes

2.2 加權(quán)基因共表達(dá)網(wǎng)絡(luò)構(gòu)建

筆者在表達(dá)量矩陣中選擇高表達(dá)的前8 000基因進(jìn)行分析,有效避免低表達(dá)基因?qū)Ψ治鲞^程的影響.軟閾值計算結(jié)果分布如圖2所示.由圖2可知,當(dāng)Power值為12時,R2數(shù)值大于0.85,且平均連接度小于100,分布較為合理,因此選擇該軟閾值進(jìn)行后續(xù)分析.

圖2 軟閾值計算結(jié)果分布Fig. 2 Soft Threshold Calculation Result Distribution

WGCNA結(jié)果如圖3所示.在動態(tài)剪切算法分析中,共獲得9個基因模塊.其中:Turquoise模塊包含的基因數(shù)量最多,共3 386個;其次為Blue模塊,包含1 689個基因;Pink模塊最少,僅包含54個基因;Grey模塊功能為回收未能分配的基因,包含880個基因;其余模塊中所包含的基因數(shù)量為100～900個.在9個基因模塊中,選擇同時滿足|r|>0.5,P<0.05的Turquoise模塊作為本研究的相關(guān)特異性模塊,并對模塊中的基因進(jìn)行過濾篩選,要求基因與分組信息的關(guān)聯(lián)性同時滿足Membership值大于0.8且Gene Significance值大于0.4,最后得到1 608個核心基因.此外,本研究也選取了Blue模塊(|r|=0.55,P=0.04),但與Turquoise模塊相比,Blue模塊中基因與分組信息的關(guān)聯(lián)性較低,因此過濾標(biāo)準(zhǔn)設(shè)定為Membership值大于0.6且Gene Significance值大于0.4,根據(jù)實際數(shù)據(jù)共篩選出核心基因1 753個.

圖3 WGCNA結(jié)果Fig. 3 Result of WGCNA

圖4 葉片與根狀莖之間的差異基因Fig. 4 Differential Expression Genes Between Leaves and Rhizomes

2.3 差異基因及相關(guān)性分析

為了進(jìn)一步驗證多糖合成通路相關(guān)基因在根狀莖與葉片中表達(dá)模式的相關(guān)性,筆者挑選出在葉片和根狀莖間表達(dá)模式高度相關(guān)的基因.葉片與根狀莖之間的差異基因如圖4所示.由圖4可知,本研究在3 236個基因中共找到88個上調(diào)基因和14個下調(diào)基因.筆者將這102個基因的序列放在NCBI數(shù)據(jù)庫中檢索,確定了基因名稱和基因?qū)?yīng)的功能.表2列出了多糖合成通路的顯著關(guān)聯(lián)基因.從表2可以看出,多糖相關(guān)合成通路中的2個關(guān)鍵基因是果糖二磷酸醛縮酶(Fructose-Bisphosphate Aldolase,FBA)和肌醇半乳糖苷合成酶(Galactinol Synthase,GOLS),它們在根狀莖與葉片之間顯著差異表達(dá).

表2 多糖合成通路顯著關(guān)聯(lián)基因

多糖合成基因在根狀莖與葉片中的相關(guān)性分析結(jié)果如圖5所示.根據(jù)基因表達(dá)矩陣,筆者計算了樣品間的皮爾森相關(guān)性系數(shù).皮爾森相關(guān)性熱圖見圖5(a).為了更好地呈現(xiàn)2組數(shù)據(jù)間關(guān)鍵基因表達(dá)模式的相關(guān)性,筆者計算了每個基因在不同分組中的相關(guān)性系數(shù)平均值,并繪制線性回歸圖(圖5(b)).在95%的置信區(qū)間內(nèi),這些基因在根狀莖和葉片之間的表達(dá)量具有很高的相關(guān)性(R=0.88,P=8.5×10-5).

圖5 多糖合成基因在根狀莖與葉片中的相關(guān)性分析Fig. 5 Correlation Analysis of Polysaccharide Synthesis Genes Between Rhizomes and Leaves

3 討論

黃精屬植物是多年生草本植物,入藥僅使用地下的根狀莖部分,其有效成分的檢驗需要挖掘根狀莖進(jìn)一步提取分析,因而無法進(jìn)行廣泛有效檢測[22-24].當(dāng)前已經(jīng)廣泛使用的轉(zhuǎn)錄組測序方案中提出的基因表達(dá)量與表型之間的調(diào)控關(guān)系,為解決該困境提供了行之有效的方法.本研究通過獲取黃精屬下多花黃精和玉竹2個物種多個樣本的轉(zhuǎn)錄組數(shù)據(jù),利用WGCNA方法,在Blue和Turquoise模塊中篩選出5 075個與分組相關(guān)的基因,另外還查找到與FBA基因和GOLS基因相對應(yīng)的多個轉(zhuǎn)錄本,并根據(jù)這些轉(zhuǎn)錄本在根狀莖與葉片的表達(dá)模式之間的皮爾森系數(shù),確定了2個基因的表達(dá)模式在器官間具有顯著關(guān)聯(lián)性(R=0.88,P=8.5×10-5),為后續(xù)基于葉片轉(zhuǎn)錄組數(shù)據(jù)預(yù)測根狀莖中多糖成分含量提供了研究基礎(chǔ).

FBA基因在多個物種中均有深入研究,并且經(jīng)過大量的實驗驗證,確認(rèn)該基因及相關(guān)蛋白是糖酵解、糖異生、磷酸戊糖通路中的重要原件,其在光合作用中也是卡爾文循環(huán)的重要參與者[25].FBA基因常常被富集到糖酵解/糖異生(ko00010)、戊糖磷酸途徑(ko00030)和半乳糖代謝(ko00052)等與多糖合成直接相關(guān)的通路中,這對于判斷黃精中多糖相關(guān)有效成分含量具有指示性作用.此外,該基因也被報道在棉花、玉米等物種不同器官中參與了多糖類化合物的合成調(diào)控[26-27].本研究結(jié)果顯示,半乳糖代謝通路(ko00052)中的GOLS基因差異表達(dá),并與葉片和根狀莖之間的表達(dá)模式高度相關(guān).目前針對GOLS基因的研究主要圍繞該蛋白處于非生物因素下對植物生長的影響方面,在不同的影響因素下,GOLS基因的應(yīng)答機(jī)制存在明顯的異質(zhì)性[28-29].這說明植株生境受到影響時,GOLS基因能夠反映植株的生長狀態(tài),為進(jìn)一步確定有效成分含量提供判斷依據(jù).綜上,FBA基因和GOLS基因在根狀莖與葉片間的表達(dá)模式為黃精屬植物多糖有效成分的觀測提供了分子指示物.

隨著多組學(xué)時代的到來,植物研究領(lǐng)域中數(shù)據(jù)量大爆炸的時代也已到來,如何更好地應(yīng)用大數(shù)據(jù)來解決科學(xué)問題,是當(dāng)前最大的科研需求.本研究采用加權(quán)共表達(dá)網(wǎng)絡(luò)分析方法整合了多個轉(zhuǎn)錄組數(shù)據(jù),從全新角度為黃精屬植物有效成分的觀測提供了分子證據(jù).盡管在當(dāng)前分析中仍然存在缺乏黃精屬植物的高質(zhì)量基因組數(shù)據(jù),以及部分現(xiàn)有數(shù)據(jù)無法使用的問題,但轉(zhuǎn)錄組數(shù)據(jù)分析方案的多樣性可基本彌補這些缺陷.雖然FBA基因參與多糖類化合物的合成調(diào)控在多個物種中均有驗證,但是在黃精屬植物中的實際調(diào)控過程不明,后續(xù)仍需要通過分子生物學(xué)實驗來進(jìn)一步驗證.