劉 暢，劉思章，于靖輝，陳 屏，雷 軍，王 義，張美萍

? 藥材與資源?

基于權重基因共表達網(wǎng)絡分析挖掘人參皂苷Rh1生物合成相關基因

劉 暢，劉思章，于靖輝，陳 屏，雷 軍，王 義，張美萍*

吉林農(nóng)業(yè)大學生命科學學院，吉林省人參基因資源開發(fā)與利用工程研究中心，吉林 長春 130118

通過權重基因共表達網(wǎng)絡分析（weighted gene co-expression network analysis，WGCNA）挖掘人參皂苷Rh1生物合成相關基因。以種植于吉林省人參主產(chǎn)區(qū)304個農(nóng)家品種的轉錄組測序數(shù)據(jù)為基礎，以人參皂苷Rh1含量為表型進行差異表達基因的挖掘。通過WGCNA獲得6個與人參皂苷Rh1含量密切相關的基因共表達模塊，根據(jù)關聯(lián)分析結果確定Green模塊為與人參皂苷Rh1含量顯著相關的關鍵模塊。功能注釋結果顯示Green模塊內(nèi)的基因富集到mRNA結合、蛋白修飾等多個相關途徑。通過構建基因互作網(wǎng)絡篩選獲得7個核心基因，功能預測表明這些基因可能在人參皂苷Rh1的生物合成途徑中起著重要作用。對7個候選基因進行茉莉酸甲酯（methyl jasmonate，MeJA）體外調(diào)控分析。對人參不定根進行MeJA誘導處理，隨MeJA的誘導，、、、的表達量與人參皂苷Rh1的含量同時呈現(xiàn)顯著變化，但只有的表達趨勢與人參皂苷Rh1在MeJA誘導下的表達趨勢相同，表明這條候選基因可能與人參皂苷Rh1的合成密切相關。為人參皂苷Rh1生物合成途徑的解析提供了理論基礎。

人參；人參皂苷Rh1；差異表達基因；權重基因共表達網(wǎng)絡分析；茉莉酸甲酯

人參C. A. Meyer是一種藥用價值極高的五加科人參屬多年生草本植物，廣泛應用于中藥領域。其根、莖、葉和花等部位都含有豐富的生物活性成分，如人參皂苷、多糖、維生素、氨基酸等。人參皂苷作為最重要的活性成分之一，具有提高免疫力、抗氧化、改善中樞神經(jīng)系統(tǒng)紊亂等多種藥理作用[1]。目前已從人參中分離鑒定出多種不同的人參皂苷，如人參皂苷Rh1、Rc、Rd、Rf等，這些不同種類的人參皂苷具有不同的藥理活性和生物利用度。隨著生物信息學技術的發(fā)展，利用基因組學、轉錄組學、代謝組學等手段對人參遺傳信息進行解析已經(jīng)成為人參皂苷合成途徑研究的重要手段。這些研究不僅揭示了人參皂苷合成途徑的分子機制，還為人參的遺傳育種提供了重要的理論基礎。

人參皂苷Rh1是一種四環(huán)三萜達瑪烷型皂苷單體，按結構分類屬于原三醇型人參皂苷，是一種具有廣泛藥理活性的稀有皂苷[2-3]。已有研究表明，人參皂苷Rh1可以抑制腫瘤細胞的擴散，對神經(jīng)系統(tǒng)和免疫系統(tǒng)也有積極的影響，具有改善哮喘[4]、抗癌[5]、預防肝損傷[6]多種藥理作用。人參皂苷Rh1的生物合成途徑包括3個部分：上游人參皂苷母核苷元前體物的合成、中游人參皂苷母核苷元碳骨架延長以及下游母核苷元形成及后修飾。其中，上游人參皂苷母核苷元前體物的合成是人參皂苷Rh1合成途徑中的關鍵步驟之一。這些前體物質(zhì)可以通過甲羥戊酸（mevalonic acid，MVA）和甲基赤蘚糖醇磷酸化（methylerythritol phosphate pathway，MEP）途徑合成，其中異戊二烯焦磷酸（isopentenyl pyrophosphate，IPP）和二甲基烯丙基焦磷酸（dimethylallyl pyrophosphate，DMAPP）是合成人參皂苷母核苷元骨架的重要前體物。中游人參皂苷母核苷元碳骨架延長是人參皂苷Rh1合成途徑的第二個關鍵步驟。IPP和DMAPP在多種萜類合酶的催化下，依次形成牻牛兒基焦磷酸（geranyl pyrophosphate）、法尼基焦磷酸（farnesyl pyrophosphate）和角鯊烯（squalene），完成了人參皂苷母核苷元骨架的延長。角鯊烯通過角鯊烯環(huán)氧酶形成不同人參皂苷母核苷元的共同前體物質(zhì)2,3-氧化角鯊烯，為人參皂苷Rh1的合成提供了重要的前體物質(zhì)。下游母核苷元形成及后修飾是人參皂苷Rh1合成途徑的最后一個步驟。這一步驟主要包括一系列的細胞色素氧化酶（P450）、糖基轉移酶等酶的修飾。通過這些酶的作用，人參皂苷Rh1的結構得到了進一步的修飾和調(diào)整，最終形成了人參皂苷Rh1單體[7-9]。

轉錄組測序是目前廣泛應用于生命科學研究中的一種高通量技術。它可以幫助研究人員準確鑒定和獲得性狀功能基因。然而，由于基因的表達與性狀間存在復雜的關聯(lián)，因此我們需要發(fā)現(xiàn)一種更為精準的方法來挖掘控制目標性狀的候選基因[10-11]。差異表達基因（differentially expressed genes，DEGs）是指當基因在RNA水平處于不同環(huán)境、壓力、和時間等條件時，表達有顯著性差異的基因[12]。近年來差異表達基因相關分析在藥用植物研究領域中廣泛應用，對控制相關生物學基因性狀的候選基因鑒定以及植物品種選育具有重要意義[13-17]。權重基因共表達網(wǎng)絡分析（weighted gene co-expression network analysis，WGCNA）作為一種新的基因互作網(wǎng)絡建模方法，能夠將基因分為不同的模塊，并與表型相關聯(lián)，從而挖掘出關鍵模塊中的基因進行深入研究[18-19]。WGCNA已經(jīng)被廣泛應用于藥用植物領域[20-22]。然而，在人參研究領域中，尚未有研究使用WGCNA的分析方法挖掘與人參皂苷生物合成相關的基因。人參皂苷是人參中最主要的有效成分之一，對于揭示其生物合成機制具有重要意義。因此，使用WGCNA方法，挖掘與人參皂苷合成相關的基因，將有助于深入了解人參皂苷的生物合成機制，提高人參的遺傳育種水平。

本研究以吉林省人參主產(chǎn)區(qū)304個農(nóng)家品種的轉錄組測序數(shù)據(jù)為基礎，并挖掘了與人參皂苷Rh1差異表達相關的基因，構建了WGCNA共表達網(wǎng)絡。通過將共表達網(wǎng)絡中的基因劃分為不同的基因模塊，再通過與人參皂苷Rh1含量的相關性分析進一步確定關鍵模塊。篩選出人參皂苷Rh1的生物合成相關中的核心基因，該核心基因在人參皂苷Rh1的生物合成途徑中具有重要作用，本研究還對這些核心基因的表達模式進行了分析，為深入了解人參皂苷Rh1的生物合成途徑提供了參考。同時為實現(xiàn)利用植物細胞工廠定向化生產(chǎn)人參皂苷Rh1奠定了基礎。為其他藥用植物重要性狀相關基因挖掘提供參考。

1 儀器與材料

1.1 儀器

Waters 2695型高效液相色譜儀，美國Waters公司；電子分析天平，上海精密科學儀器有限公司；超聲清洗機，寧波新芝生物科技有限公司；37 ℃恒溫烘箱，Thermo公司；恒溫搖床，上海智誠分析儀器制造有限公司；超低溫冰箱，Thermo公司；水浴鍋，上海一恒科技有限公司；旋轉蒸發(fā)器，上海亞榮生化儀器廠；熒光定量PCR儀，賽默飛世爾科技有限公司。

1.2 材料

對照品人參皂苷Rh1（批號為B21061），質(zhì)量分數(shù)≥98%，上海源葉生物科技有限公司；甲醇（批號M813895）購自上海麥克林生化科技股份有限公司；色譜甲醇批號為MS1922-801，色譜乙腈批號為AS1122-801，美國天地有限公司；屈臣氏蒸餾水，北京屈臣氏蒸餾水有限公司；茉莉酸甲酯（methyl jasmonate，MeJA）批號為M46630，上海吉至生化科技有限公司。

本研究所用304份人參材料為吉林省人參基因資源開發(fā)與利用工程研究中心提供，來源于種植在吉林省人參主產(chǎn)區(qū)的四年生人參主根，經(jīng)吉林農(nóng)業(yè)大學生命科學學院王義教授鑒定為五加科人參屬植物人參C. A. Meyer；本研究所用數(shù)據(jù)庫為吉林省人參基因資源開發(fā)與利用工程研究中心根據(jù)種植于吉林省人參主產(chǎn)區(qū)的304個農(nóng)家品種（S1～S304）4年生人參樣本根部轉錄組測序數(shù)據(jù)所建立的轉錄組序列及表達量數(shù)據(jù)庫（NCBI: SRR1313640～SRR132505）[23]，包括304個農(nóng)家品種4年生人參樣本根部人參皂苷Rh1含量。

2 方法

2.1 人參皂苷Rh1生物合成差異基因的挖掘

2.1.1 樣本數(shù)據(jù)的預處理 將數(shù)據(jù)庫中的304個樣本按照人參皂苷Rh1含量進行排序，分別取人參皂苷Rh1含量較高的前25個樣本及人參皂苷Rh1含量較低的后25個樣本，設置為高含量組及低含量組，繪制布魯克斯圖，探討樣本間差異性。使用Perl語言分別調(diào)取2組樣本的轉錄組序列及表達量，用于差異基因挖掘。

2.1.2 差異基因的挖掘 調(diào)用R語言中的DESeq2包對上述高低2組樣本中的數(shù)據(jù)樣本進行差異基因分析，將≤0.05的基因視為差異基因。利用TBtools繪制差異基因的火山圖及熱圖，并通過主成分分析（principal component analysis，PCA）進一步分析驗證高低兩組樣本間的差異性。

2.2 權重共表達網(wǎng)絡的構建

利用權重共表達網(wǎng)絡分析基因與人參皂苷Rh1含量之間的關系，以“2.1.2”項中挖掘獲得的差異基因為輸入數(shù)據(jù)，調(diào)用R語言中的WGCNA包構建共表達網(wǎng)絡。利用Pearson相關矩陣和網(wǎng)絡拓撲分析分別計算基因相關性和軟閾值能力。選擇合適的軟閾值對原始矩陣進行轉化，得到鄰接矩陣，然后將鄰接矩陣轉化為拓撲重疊矩陣，并運用動態(tài)切割法進行基因聚類及劃分模塊。設定模塊最少基因數(shù)為5，相似模塊合并閾值為0.85。繪制樣本表達模式熱圖，并通過模塊特征值來展示模塊基因在各個樣本中的表達模式，根據(jù)模塊特征向量基因分析確定與人參皂苷Rh1含量顯著相關的特異性模塊，選擇相應的模塊進行深入研究。

2.3 關鍵模塊的功能富集分析

為進一步分析關鍵模塊的功能，使用Blast2GO及京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）在線平臺對關鍵模塊內(nèi)的基因進行基因本體論（GeneOntology，GO）功能注釋及KEGG通路分析。經(jīng)過多重檢驗校正后，以≤0.05為閾值，滿足此條件的定義為顯著富集的GO條目及KEGG 通路。

2.4 基因互作網(wǎng)絡構建及核心基因的挖掘

為獲得關鍵模塊中的核心基因，為獲得關鍵模塊中的核心基因，利用Cytoscape3.9.0軟件中的cytoHubba插件，選擇MCC算法，對關鍵模塊內(nèi)的基因進行可視化處理，每個節(jié)點代表1個基因，連線代表基因間的關系。根據(jù)互作網(wǎng)絡快速篩選出Weight值及模塊連通性都較高的基因，作為關鍵模塊核心基因。利用NCBI數(shù)據(jù)庫（https://www. ncbi.nlm.nih.gov/）獲取這些核心基因的相關信息，并在韓國人參基因組[24]中查找其同源基因的注釋結果，進一步預測其功能。

2.5 MeJA誘導人參皂苷Rh1

將2 g在固體B5培養(yǎng)基上生長的人參不定根接種于液體B5培養(yǎng)基中，在22 ℃、110 r/min的搖床中黑暗培養(yǎng)23 d，在對數(shù)生長期后，向培養(yǎng)基中添加200 μmol/L MeJA進行誘導處理，誘導時間設置為6、12、24、36、48、60、72、84、96、108、120 h，在每個時間點收獲材料，同時以MeJA處理0 h的材料作為對照組收獲，對上述每組材料進行3次生物學重復。取樣時，將0.5 g新鮮樣品進行液氮速凍，保存在?80 ℃冰箱中，將剩余樣品在37 ℃恒溫烘箱烘干至恒定質(zhì)量后提取人參皂苷。利用索氏提取法進行皂苷提取。

2.6 人參皂苷Rh1的測定

2.6.1 對照品溶液的制備 分別稱取人參皂苷Rh1對照品適量，加入甲醇定容，經(jīng)0.22 μm膜濾過，制成含人參皂苷Rh1質(zhì)量濃度為0.075 mg/mL的對照品溶液。

2.6.2 供試品溶液的制備 取“2.5”項所述方法的人參皂苷提取液，用0.22 μm有機濾膜濾過，得供試品溶液。

2.6.3 色譜條件[25]高效液相色譜系統(tǒng)是Waterse2695，色譜柱為Waters C18柱。流動相組為屈臣氏蒸餾水（A）和色譜乙腈（B），梯度洗脫，進樣量為10 μL，柱溫30 ℃，流動相體積流量為1.0 mL/min。檢測波長為203 nm。梯度洗脫：0～40 min，18%～21% B；40～42 min，21%～26% B；42～46 min，26%～32% B；46～66 min，32%～33.5% B；66～71 min，33.5%～38% B；71～86 min，38%～65% B；86～91 min，65% B；91～96 min，65%～85% B；96～103 min，85% B；103～105 min，85%～18% B；105～120 min，82% B。

2.6.4 線性關系考察 精密吸取上述混合對照品溶液，按照“2.6.3”項色譜條件進行測定，以色譜峰峰面積為縱坐標（），以人參皂苷Rh1對照品質(zhì)量濃度為橫坐標（），繪制標準曲線，進行線性回歸，得到回歸方程，相關系數(shù)（2）和線性范圍分別為人參皂苷Rh1＝347 766.69－6 332.14，2＝0.999 65，線性范圍6.0～0.375 μg/mL。

2.6.5 穩(wěn)定性試驗 精密吸取適量的人參皂苷Rh110 μL對照品溶液，依照“2.6.3”項色譜條件分別在0、4、8、12、16、20、24 h進樣，測得人參皂苷Rh1的峰面積的RSD為0.80%，表明儀器穩(wěn)定性良好。

2.6.6 精密度試驗 精密吸取適量的人參皂苷Rh110 μL對照品溶液，按照“2.6.3”項色譜條件進行分析，連續(xù)進樣6次，測得人參皂苷Rh1的峰面積的RSD為0.89%，表明儀器精密度良好。

2.6.7 重復性試驗 精密吸取人參皂苷Rh1對照品溶液6份，按照按照“2.6.3”項色譜條件進行分析，測定人參皂苷Rh1的峰面積，通過線性回歸方程計算其質(zhì)量分數(shù)RSD值為0.96%，結果表明該方法重復性較好。

2.6.8 加樣回收率的測定 精密吸取已知質(zhì)量分數(shù)的人參皂苷Rh1對照品溶液，按照“2.6.2”項和“2.6.3”項方法制備供試液和進樣，計算得到人參皂苷Rh1的加樣回收率為99.63%，RSD值為1.14%，數(shù)據(jù)表明測試結果的準確度良好。

2.6.9 人參皂苷Rh1的測定 根據(jù)“2.6.3”項色譜條件對MeJA處理的人參不定根中的人參皂苷Rh1進行含量測定。

2.7 RNA提取、cDNA文庫構建及數(shù)據(jù)組裝

使用Trizol試劑盒提取人參不定根的總RNA，按照Jiang等[26]的方法構建RNA-seq文庫。文庫經(jīng)過鑒定、定量和復用，使用HiSeq X Ten（Illumina，Inc.，San Diego，美國）進行測序，使用Trinity （Version 2.14.0）軟件，以人參unigene的轉錄組為參考，從過濾后的讀數(shù)中組裝單個基因轉錄本。使用RSEM（Version 1.3.3）軟件對單個轉錄本的表達量和基因的整體表達量進行了量化。轉錄本的表達量以每百萬轉錄本（TPM）為單位，用于進一步分析。PgRh候選基因轉錄本的表達量是從每個時間點取樣的不定根的每個生物重復中提取的。

2.8 人參皂苷合成候選基因在MeJA誘導下的功能驗證

為了驗證MeJA處理對人參皂苷Rh1生物合成的影響，通過檢驗，將6～120 h每個時間點的3個生物重復的人參皂苷Rh1含量與0 h時間點的人參皂苷Rh1含量進行比較。使用TRIpure Reagent Total RNA Extraction Reagent從每個時間點采樣的3個生物重復中分離出總RNAs。PgRh候選基因轉錄本的表達量是從每個時間點取樣的不定根的每個生物重復中提取的。通過檢驗，比較每個時間點MeJA處理的樣本和0 h時間點對照樣本的基因表達，確認MeJA處理對PgRh候選基因表達的影響。為了進一步確認PgRh候選基因在人參皂苷Rh1生物合成中的作用，分析其單個轉錄本的表達，并對每個時間點（包括0 h）的單個轉錄本的表達與每個時間點的人參皂苷Rh1含量進行了Pearson相關分析。如果一個PgRh候選基因的轉錄本表達量與人參皂苷Rh1含量在≤0.05的雙尾顯著性水平下有相關性，則認為拼接到該轉錄本的基因參與人參皂苷Rh1的生物合成。

3 結果與分析

3.1 人參皂苷Rh1生物合成差異基因的挖掘

3.1.1 樣本數(shù)據(jù)的預處理 首先對304個農(nóng)家品種4年生根的人參皂苷Rh1含量分布情況進行了統(tǒng)計。如圖1所示，304個農(nóng)家品種4年生根中人參皂苷Rh1皂苷含量呈正態(tài)分布，存在較大差異。將304個農(nóng)家品種4年生根中的人參皂苷Rh1含量進行分組，每組25個農(nóng)家品種，含量高的為高組，含量低的為低組。圖2顯示了這2組人參皂苷Rh1皂苷含量顯著不同（≤0.05），可用于進一步進行差異表達基因分析。

3.1.2 差異基因的挖掘 將“3.1.1”項得到的差異基因進行了火山圖（圖3）和熱圖的繪制（圖4）。在703個差異基因中，有578個基因上調(diào)，而125個基因下調(diào)。對這些差異基因進一步進行了PCA分析。PCA結果（圖5）表明高低含量2組之間存在明顯差異，這為進一步的分析提供了基礎。

圖1 304份樣本人參皂苷Rh1含量的頻率分布直方圖

圖2 人參皂苷Rh1含量差異分組

圖3 人參皂苷Rh1差異基因火山圖

圖4 人參皂苷Rh1差異基因表達量熱圖

圖5 人參皂苷Rh1差異基因PCA分析

3.1.3 差異基因的功能富集分析 對“3.1.2”項獲得的差異基因進行了GO功能注釋和KEGG信號通路分析。進一步總結了差異基因的GO功能注釋和KEGG信號通路分析結果。根據(jù)GO功能注釋的結果（圖6），發(fā)現(xiàn)核心模塊基因的功能主要涉及核酸結合（GO: 0003676）、RNA結合（GO: 0003723）、mRNA結合（GO: 0003729）、細胞液（GO: 0005829）、核質(zhì)（GO: 0005654）、染色質(zhì)組織（GO: 0006325）、細胞蛋白修飾過程（GO: 0006464）、質(zhì)體（GO: 0009536）、轉移酶活性（GO: 0016740）、蛋白修飾過程（GO: 0036211）和催化活性（GO: 0140096）等多個GO term。KEGG富集結果顯示（圖7），差異基因的代謝通路主要集中在碳水化合物代謝、轉運、信號轉導、運輸和分解代謝等途徑。

圖6 核心模塊的GO富集分析

圖7 核心模塊KEGG富集功能分析

3.2 權重共表達網(wǎng)絡的構建

3.2.1 軟閾值的確定 根據(jù)加權共表達網(wǎng)絡的要求，利用無尺度拓撲結構準則來確定軟閾值β。運用R語言中WGCNA包中的pickSoftThreshold函數(shù)計算軟閾值。數(shù)值范圍設置為1～20，計算相關系數(shù)的平方和基因的平均連接度，以確定最佳的軟閾值。根據(jù)本研究的選擇標準，選擇擬合曲線第一次超過0.85時的power值來篩選共表達模塊，當power＝12時，無標度拓撲擬合指數(shù)曲線在達到較高值后開始變平（圖8）。因此，選擇12作為后續(xù)分析的power值。

圖8 軟閾值的確定與平均連通度圖

3.2.2 構建加權共表達網(wǎng)絡和基因模塊的確定 在網(wǎng)絡構建和模塊劃分的過程中，選擇了β值為12為最佳軟閾值，并根據(jù)基因拓撲重疊矩陣的相異度進行聚類。通過聚類分析，得到了基因拓撲重疊聚類樹和相關性熱圖。根據(jù)圖9的結果，將與人參皂苷Rh1相關的703個差異基因分為了6個模塊。在這6個模塊中，turquoise模塊擁有最多的基因數(shù)量，共有333條基因，而red模塊則只包含6條基因，是基因數(shù)量最少的模塊。一般WGCNA中grey模塊被視為不能分配給任何其他模塊的基因集合。

3.2.3 模塊與人參皂苷Rh1的相關性系數(shù)計算及關鍵模塊確定 通過計算基因模塊與人參皂苷Rh1的相關性系數(shù)，確定了關鍵模塊。根據(jù)圖10的相關性熱圖，發(fā)現(xiàn)在6個基因模塊中，有4個基因模塊與人參皂苷Rh1呈現(xiàn)正相關，而另外2個基因模塊與人參皂苷Rh1呈現(xiàn)負相關。在這些基因模塊中，Green模塊的相關系數(shù)達到了0.74，值為0.008。因此，選擇Green模塊中的基因進行研究。

圖9 6個基因模塊中基因數(shù)量分布

圖10 6個基因模塊與人參皂苷Rh1含量關聯(lián)熱圖

3.2.4 基因互作網(wǎng)絡與核心基因挖掘 可視化得到的核心基因進一步在NCBI數(shù)據(jù)庫與韓國基因組數(shù)據(jù)庫進行比對（表1）。用Cytoscape軟件對基因互作網(wǎng)絡進行可視化，根據(jù)weight值與基因之間的連通性建立基因互作網(wǎng)絡（圖11）。

3.3 MeJA誘導處理后人參皂苷Rh1含量變化

本實驗在MeJA誘導處理后的不同時間內(nèi)，通過高效液相色譜檢測人參不定根中人參皂苷Rh1皂苷含量。根據(jù)圖12結果，與對照組（0 h）相比，MeJA處理48 h后人參皂苷Rh1含量顯著變化，表明人參不定根中的人參皂苷Rh1對MeJA有響應作用。MeJA顯著促進了人參皂苷Rh1的生物合成，隨著MeJA處理時間的增加，人參皂苷Rh1含量顯著增加。在處理84 h時，人參皂苷Rh1含量達到最高值1.09 mg/g，相較于對照組，人參皂苷Rh1含量增加了3.19倍，說明MeJA對人參不定根中人參皂苷Rh1的合成起到了促進作用。

在MeJA處理的不同時間點（48、72、84、96、108、120 h）下，人參皂苷Rh1的含量與對照組相比均發(fā)生了變化。人參皂苷Rh1在MeJA處理下的48 h時與對照組相比發(fā)生顯著差異，在84 h時人參皂苷Rh1含量達到頂峰。由圖13可知，在MeJA處理72 h時，人參皂苷Rh1候選基因、和的表達量與對照組也發(fā)生了顯著變化。同樣地，在MeJA處理108 h時，人參皂苷Rh1候選基因的表達量與對照組也發(fā)生了顯著變化。但是在MeJA處理48 h和84 h時，、和的表達情況與人參皂苷Rh1含量沒有相同的趨勢，以上結果說明只有候選基因在人參皂苷Rh1合成過程中可能發(fā)揮重要作用。

表1 Green模塊中核心基因的功能注釋

圖11 核心模塊中核心基因互作網(wǎng)絡分析

與對照組（0 h）比較：*P＜0.05

圖13 人參皂苷Rh1合成候選基因在MeJA處理下的表達分析

4 討論

由于生物系統(tǒng)的復雜性，許多關鍵基因在植物次級代謝物的生成過程中尚未明確，這些基因廣泛分布在植物染色體中，因此挖掘與表征它們具有一定的挑戰(zhàn)性。人參中的主要活性成分，人參皂苷的合成路徑及其機制，一直是人參研究的關鍵話題。

目前已經(jīng)對基因家族[27]、基因家族[28]、基因家族[29]、基因家族[30]、P基因家族[31]以及基因家族[32]等在人參皂苷生物合成過程中起作用的多個基因家族進行了全面的分析和研究，這極大地加強了對人參皂苷合成機制的理解，并為其提供了寶貴的基因資源。

隨著高通量測序技術的不斷突破與創(chuàng)新和成本的不斷降低，多樣本的轉錄組測序在系統(tǒng)研究生命科學問題中已被廣泛應用。然而，如何從這些龐大數(shù)據(jù)中挖掘具有生物學意義的信息已成為轉錄組分析的核心問題[33-34]。在這種背景下，涌現(xiàn)出了大量的生物信息分析工具和方法。WGCNA作為一種系統(tǒng)生物學方法，可以闡述基因在不同樣本間關聯(lián)模式。WGCNA將大量的基因分組為若干模塊，這些模塊內(nèi)的基因表達模式具有相同的表達模式。通過與目標性狀進行相關性分析確定關鍵模塊，并根據(jù)連通性篩選模塊內(nèi)的核心基因，以此可以快速定位與目標表型顯著關聯(lián)的基因。對共表達基因模塊的鑒定、與表型性狀的關聯(lián)以及關鍵樞紐基因的定位具有重要意義[35]。

Yao等[36]對重組自交系的WDD01514（E1）、ZYD00463（E2）及其重組自交系的后代(E23和E171)在大豆發(fā)育的2個關鍵時期進行RNA-seq、WGCNA等分析鑒定出6個與種子含油量和種子重量相關的模塊后，并結合后續(xù)生物信息學分析、q-PCR等驗證6個可能參與種子形成和含油量積累的關鍵基因，為提高大豆的產(chǎn)量和種子油量提供理論支持。Zhu等[37]利用RNA-seq和差異基因分析，在水稻轉錄組的28 432個表達基因中鑒定出457個響應鹽脅迫的核心基因，進一步利用WGCNA分析、GO功能注釋、KEGG分析以及qRT-PCR鑒定得到編碼、、、、、和轉錄因子等不同家族蛋白的重要樞紐基因，推動了培養(yǎng)水稻耐鹽新品種的研究。Li等[38]對C2P5A和C2P5B 2種不同的棉花品種在發(fā)育的3個不同階段（花粉母細胞、四分體和單核階段）進行RNA-seq、WGCNA分析、差異基因分析等，挖掘得到3個與細胞質(zhì)雄性不育（CMS）高度相關的模塊和7個與生長發(fā)育相關的基因，為棉花育種和生產(chǎn)提供依據(jù)。

次生代謝物的生物合成與調(diào)控構建了一個多層次的網(wǎng)絡，這需要多種關鍵誘導子（如MeJA和水楊酸）和信號通道的合作。已有證據(jù)確認，MeJA是一種強大的誘導子，能夠在多種藥用植物中刺激高價值藥用次生代謝物的生產(chǎn)[39]。對于人參培養(yǎng)細胞和人參不定根，MeJA發(fā)揮了顯著作用，增加了人參皂苷含量，并且能夠上調(diào)相關的代謝酶基因表達[40-41]。Kim等[42]發(fā)現(xiàn)，MeJA能誘發(fā)、等基因的表達，同時也刺激了皂苷的積累。因此，在本研究中選擇使用MeJA誘導處理人參不定根，并進行了深入分析。

目前，人參皂苷生物合成關鍵酶基因，如[43-45]、[46]、[47]、[48-50]、[51]，均響應MeJA調(diào)控，皂苷含量和基因表達量都發(fā)生了變化，以此證明這些基因參與人參皂苷的生物合成。

本研究結果顯示，可能的人參皂苷Rh1皂苷候選基因、、、在MeJA誘導下表達量均發(fā)生顯著變化，但只有的表達趨勢與人參皂苷Rh1在MeJA誘導下的表達趨勢相同，表明這條候選基因可能與人參皂苷Rh1的合成密切相關。

總的來說，通過WGCNA分析方法，挖掘出了與人參皂苷Rh1含量高度相關的關鍵基因模塊，成功地鑒定了與人參皂苷Rh1生物合成途徑相關的關鍵候選基因。這些結果為進一步研究人參皂苷Rh1的生物合成機制提供了重要的線索，也為人參的遺傳改良和代謝工程提供了可能的目標。

通過對304個人參樣本中人參皂苷Rh1含量進行分組，利用差異基因分析的方法成功獲得703個差異基因?；?04份材料中703個差異基因的表達量與人參皂苷Rh1的含量，通過WGCNA分析共獲得6個基因表達模塊。通過模塊與表型的關聯(lián)分析結果顯示Green模塊與人參皂苷Rh1含量的關聯(lián)性最高。對Green模塊進行GO富集分析和KEGG富集分析結果顯示Green模塊中的基因可以富集到細胞代謝、RNA結合和細胞質(zhì)修飾等GO term；同時參與碳水化合物代謝、轉運、信號轉導、運輸和分解代謝等途徑。因此，Green模塊中的基因參與植物代謝產(chǎn)物的合成。通過Cytoscape軟件對Green模塊內(nèi)的基因進行了可視化處理，并篩選出了7個核心基因。通過查詢韓國人參基因組注釋。結果顯示7條核心基因注釋到了人參皂苷合成基因和生長素響應因子等生物功能。

采用MeJA誘導處理人參不定根后，發(fā)現(xiàn)、、、4條合成候選基因的表達量變化趨勢只有與人參皂苷Rh1在MeJA誘導處理下的表達量變化趨勢相同，因此這條基因可能在人參皂苷Rh1皂苷合成過程中發(fā)揮重要作用。

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

[1] Ratan Z A, Haidere M F, Hong Y H,. Pharmacological potential of ginseng and its major component ginsenosides [J]., 2021, 45(2): 199-210.

[2] 孫小單, 王天鳴, 李慧, 等. 人參皂苷Rh2抑制人非小細胞肺癌細胞增殖的機制研究[J]. 中草藥, 2022, 53(2): 441-448.

[3] 畢云楓, 姜珊, 鄭明珠, 等. 稀有原人參三醇型皂苷的生物轉化研究進展 [J]. 中草藥, 2017, 48(19): 4120-4125.

[4] Jin Y J, Tangchang W, Kwon O S,. Ginsenoside Rh1ameliorates the asthma and allergic inflammation via inhibiting Akt, MAPK, and NF-κB signaling pathwaysand[J]., 2023, 321: 121607.

[5] Yoon J H, Choi Y J, Lee S G. Ginsenoside Rh1suppresses matrix metalloproteinase-1 expression through inhibition of activator protein-1 and mitogen-activated protein kinase signaling pathway in human hepatocellular carcinoma cells [J]., 2012, 679(1/2/3): 24-33.

[6] Su W Y, Fan M L, Li Y,. 20()-ginsenoside Rh1alleviates T2DM induced liver injury via the Akt/FOXO1 pathway [J]., 2022, 20(9): 669-678.

[7] Yang J L, Hu Z F, Zhang T T,. Progress on the studies of the key enzymes of ginsenoside biosynthesis [J]., 2018, 23(3): 589.

[8] Hou M Q, Wang R F, Zhao S J,. Ginsenosides ingenus and their biosynthesis [J]., 2021, 11(7): 1813-1834.

[9] Yin J X, Zhang D H, Zhuang J J,. Study on the correlation between gene expression and enzyme activity of seven key enzymes and ginsenoside content in ginseng in over time in Ji’an, China [J]., 2017, 18(12): 2682.

[10] Velculescu V E, Zhang L, Zhou W,. Characterization of the yeast transcriptome [J]., 1997, 88(2): 243-251.

[11] 李欣, 李小俊, 陳曉麗, 等. 轉錄組數(shù)據(jù)分析與功能基因挖掘 [J]. 畜牧獸醫(yī)學報, 2019, 50(3): 474-484.

[12] 章魯瑤, 張洋, 張書筠, 等. 類風濕關節(jié)炎免疫相關差異表達基因分析 [J]. 牡丹江醫(yī)學院學報, 2023, 44(1): 44-50.

[13] Chaudhary S, Sharma P C. DeepSAGE based differential gene expression analysis under cold and freeze stress in seabuckthorn (L.) [J]., 2015, 10(3): e0121982.

[14] 吳萍, 王曉宇, 李青苗, 等. 基于高通量轉錄組測序的白芷差異表達基因分析 [J]. 西南農(nóng)業(yè)學報, 2020, 33(2): 233-240.

[15] 肖麗君, 陳謙明, 吳玉琪, 等. 基于轉錄組測序分析芍藥苷誘導下口腔黏膜上皮細胞差異表達基因 [J]. 口腔醫(yī)學研究, 2023, 39(6): 527-533.

[16] Chang Y J, Wang M Z, Li J,. Transcriptomic analysis reveals potential genes involved in tanshinone biosynthesis in[J]., 2019, 9(1): 14929.

[17] Lei H Y, Niu T Z, Song H F,. Comparative transcriptome profiling reveals differentially expressed genes involved in flavonoid biosynthesis between biennial and triennial[J]., 2021, 161: 113217.

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

[19] Liu W, Li L, Ye H,. Weighted gene co-expression network analysis in biomedicine research [J]., 2017, 33(11): 1791-1801.

[20] Liu D F, Wang Z J, Zhong L,. Targets and potential mechanism ofin treatment of primary hepatocellular carcinoma based on bioinformatics analysis [J]., 2022, 2022: 8762717.

[21] Zeng G P, Li Z, Zhao Z. Analysis of weighted gene co-expression network of triterpenoid-related transcriptome characteristics from different strains of[J]., 2021, 11(1): 18207.

[22] 鄔東東, 楊啟帆, 方柔柔, 等. 基于網(wǎng)絡藥理學聯(lián)合加權基因共表達網(wǎng)絡分析黃精治療代謝性相關脂肪性肝病的作用機制 [J]. 現(xiàn)代藥物與臨床, 2023, 38(7): 1606-1614.

[23] Liu S Z, Jiang Y, Wang Y F,. Genetic and molecular dissection of ginseng (Mey.) germplasm using high-density genic SNP markers, secondary metabolites, and gene expressions [J]., 2023, 14: 1165349.

[24] Kim N H, Jayakodi M, Lee S C,. Genome and evolution of the shade-requiring medicinal herb[J]., 2018, 16(11): 1904-1917.

[25] 于靖輝. 人參內(nèi)生菌廣泛靶向代謝組學分析及其對人參不定根的影響 [D]. 長春: 吉林農(nóng)業(yè)大學, 2023.

[26] Jiang Y, Liu S Z, Li L,. Transcriptome and phenotype integrated analysis identifies genes controlling ginsenoside Rb1biosynthesis and reveals their interactions in the process in[J]., 2022, 23(22): 14016.

[27] Liu Q, Sun C Y, Han J Z,. Identification, characterization and functional differentiation of thegene family and its roles in response to cold stress in ginseng,C.A. Meyer [J]., 2020, 15(6): e0234423.

[28] Chen J, Zhou Y H, Zhang Q,. Structural variation, functional differentiation and expression characteristics of the AP2/gene family and its response to cold stress and methyl jasmonate inC.A. Meyer [J]., 2020, 15(3): e0226055.

[29] Zhu L, Zhao M Z, Chen M Y,. Thegene family and its response to saline stress in Jilin ginseng,C.A. Meyer [J]., 2020, 295(4): 877-890.

[30] Li H J, Chen J, Zhao Q,. Basic leucine zipper () transcription factor genes and their responses to drought stress in ginseng,C.A. Meyer [J]., 2021, 22(1): 316.

[31] Wang Y, Li X Y, Lin Y P,. Structural variation, functional differentiation, and activity correlation of the cytochrome P450 gene superfamily revealed in ginseng [J]., 11: 170106.

[32] Wang N, Wang K Y, Li S K,. Transcriptome-wide identification, evolutionary analysis, and GA stress response of thegene family inC. A. Meyer [J]., 2020, 9(2): 190.

[33] Tian J D, Ma K S, Saaem I. Advancing high-throughput gene synthesis technology [J]., 2009, 5(7): 714-722.

[34] 王艷花, 劉景森, 李加納. 整合GWAS和WGCNA篩選鑒定甘藍型油菜生物產(chǎn)量候選基因 [J]. 作物學報, 2021, 47(8): 1491-1510.

[35] 龍佳麗, 鄒奕, 邳植, 等. 轉錄組測序揭示激素介導的信號通路參與甜菜低溫應答 [J]. 中國農(nóng)學通報, 2021, 37(31): 5-14.

[36] Yao Y J, Xiong E H, Qu X L,. WGCNA and transcriptome profiling reveal hub genes for key development stage seed size/oil content between wild and cultivated soybean [J]., 2023, 24(1): 494.

[37] Zhu M D, Xie H J, Wei X J,. WGCNA analysis of salt-responsive core transcriptome identifies novel hub genes in rice [J]., 2019, 10(9): 719.

[38] Li Y Q, Qin T F, Wei C Y,. Using transcriptome analysis to screen for key genes and pathways related to cytoplasmic male sterility in cotton (L.) [J]., 2019, 20(20): 5120.

[39] Ho T T, Murthy H N, Park S Y. Methyl jasmonate induced oxidative stress and accumulation of secondary metabolites in plant cell and organ cultures [J]., 2020, 21(3): 716.

[40] Lu J, Li J X, Wang S H,. Advances in ginsenoside biosynthesis and metabolic regulation [J]., 2018, 65(4): 514-522.

[41] Cao H Z, Nuruzzaman M, Xiu H,. Transcriptome analysis of methyl jasmonate-elicitedadventitious roots to discover putative ginsenoside biosynthesis and transport genes [J]., 2015, 16(2): 3035-3057.

[42] Kim O T, Bang K H, Kim Y C,. Upregulation of ginsenoside and gene expression related to triterpene biosynthesis in ginseng hairy root cultures elicited by methyl jasmonate [J]., 2009, 98(1): 25-33.

[43] Kim Y J, Zhang D B, Yang D C. Biosynthesis and biotechnological production of ginsenosides [J]., 2015, 33(6 Pt 1): 717-735.

[44] Lee M H, Jeong J H, Seo J W,. Enhanced triterpene and phytosterol biosynthesis inoverexpressing squalene synthase gene [J]., 2004, 45(8): 976-984.

[45] Kim T D, Han J Y, Huh G H,. Expression and functional characterization of three squalene synthase genes associated with saponin biosynthesis in[J]., 2011, 52(1): 125-137.

[46] Kim Y K, Kim Y B, Uddin M R,. Enhanced triterpene accumulation inhairy roots overexpressing mevalonate-5-pyrophosphate decarboxylase and farnesyl pyrophosphate synthase [J]., 2014, 3(10): 773-779.

[47] Han J Y, Kwon Y S, Yang D C,. Expression and RNA interference-induced silencing of the dammarenediol synthase gene in[J]., 2006, 47(12): 1653-1662.

[48] Chen S, Luo H, Li Y,. 454 EST analysis detects genes putatively involved in ginsenoside biosynthesis in[J]., 2011, 30(9): 1593-1601.

[49] Yan X, Fan Y, Wei W,. Production of bioactive ginsenoside compound K in metabolically engineered yeast [J]., 2014, 24(6): 770-773.

[50] Wang P P, Wei Y J, Fan Y,. Production of bioactive ginsenosides Rh2and Rg3by metabolically engineered yeasts [J]., 2015, 29: 97-105.

[51] Han J Y, Hwang H S, Choi S W,. Cytochrome P450 CYP716A53v2 catalyzes the formation of protopanaxatriol from protopanaxadiol during ginsenoside biosynthesis in[J]., 2012, 53(9): 1535-1545.

Mining of genes related to biosynthesis of ginsenoside Rh1 based on WGCNA in

LIU Chang, LIU Sizhang, YU Jinghui, CHEN Ping, LEI Jun, WANG Yi, ZHANG Meiping

Jilin Province Engineering Research Center of Ginseng Gene Resources Development and Utilization,College of Life Science, Jilin Agriculturial University, Changchun, 130118, China

The aim of this study is to unearth genes associated with the biosynthesis of ginsenoside Rh1employ through weighted gene co-expression network analysis (WGCNA).Utilizing transcriptome sequencing data from 304 distinct domesticated varieties of ginseng cultivated in the primary ginseng-producing region of Jilin Province, China, we conducted a differentially expressed gene mining with ginsenoside Rh1content as the phenotypic parameter.A total ofsix co-expression modules of genes intimately correlated with ginsenoside Rh1content was gained by employing WGCNA. Subsequent correlation analyses highlighted the Green module as a pivotal module significantly associated with ginsenoside Rh1content. Functional annotation results revealed that genes within the Green module were enriched in several pathways, including mRNA binding and protein modification. A total of seven core genes were selected by constructing gene interaction network, and functional prediction indicated that these genes may play crucial roles in the biosynthesis pathway of ginsenoside Rh1.modulation analysis of methyl jasmonate (MeJA) was performed on these seven candidate genes. Upon MeJA induction in ginseng adventitious roots, the expression levels of,,, andand ginsenoside Rh1content exhibited significant simultaneous changes. However, only the expression trend of compwas the same as that of ginsenoside Rh1induced by MeJA, indicating that this candidate gene may be closely related to ginsenoside Rh1synthesis.This study lays the theoretical foundation for unraveling the biosynthesis pathway of ginsenoside Rh1.

C. A. Mey.; ginsenosides Rh1; differentially expressed gene; WGCNA; methyl jasmonate

R286.12

A

0253 - 2670(2024)08 - 2723 - 11

10.7501/j.issn.0253-2670.2024.08.021

2023-10-05

吉林省科技廳項目（20180101027JC，20200403011SF）；國家自然科學基金青年基金項目（31401445）

劉 暢，碩士研究生，研究方向為人參基因組學。Tel: 15148182822 E-mail: 2454278361@qq.com

通信作者：張美萍，博士生導師，研究方向為人參功能基因組學、人參分子設計育種。Tel: 13504315977 E-mail: meiping.zhang@jlau.edu.cn

[責任編輯 時圣明]