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

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

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

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

1 材料和方法

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

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

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

1.2 構建加權基因共表達網(wǎng)絡

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

1.3 多糖合成通路相關性

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

2 結果

2.1 組裝與注釋

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

表1 多糖合成通路

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

2.2 加權基因共表達網(wǎng)絡構建

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

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

WGCNA結果如圖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模塊作為本研究的相關特異性模塊,并對模塊中的基因進行過濾篩選,要求基因與分組信息的關聯(lián)性同時滿足Membership值大于0.8且Gene Significance值大于0.4,最后得到1 608個核心基因.此外,本研究也選取了Blue模塊(|r|=0.55,P=0.04),但與Turquoise模塊相比,Blue模塊中基因與分組信息的關聯(lián)性較低,因此過濾標準設定為Membership值大于0.6且Gene Significance值大于0.4,根據(jù)實際數(shù)據(jù)共篩選出核心基因1 753個.

圖3 WGCNA結果Fig. 3 Result of WGCNA

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

2.3 差異基因及相關性分析

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

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

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

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

3 討論

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

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

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