高溫脅迫下不同玉米材料間的基因差異表達(dá)分析

曹丹，張紅梅，楊海鵬，趙歡，劉小紅

(1.西華師范大學(xué)生命科學(xué)學(xué)院，四川 南充 637000；2.山西省農(nóng)業(yè)科學(xué)院玉米研究所，山西 忻州 034000)

玉米(Zeamays)是世界上最重要的3大谷類作物之一，與水稻和小麥并列，在世界農(nóng)業(yè)中占有重要地位[1].高溫是限制全球玉米產(chǎn)量的最關(guān)鍵的非生物脅迫因素之一，世界上許多國家的玉米生產(chǎn)都遭受高溫危害[2-5].在中國，黃淮海地區(qū)是我國夏玉米集中產(chǎn)區(qū)，幾乎每年都會(huì)遭遇高溫危害.近年來，伴隨著全球氣候變暖，干旱和高溫協(xié)同脅迫已經(jīng)成為該地區(qū)玉米生產(chǎn)中頻繁遇到的一種自然災(zāi)害，且這一危害幾乎發(fā)生在玉米整個(gè)生育期，嚴(yán)重影響著玉米產(chǎn)量和品質(zhì)[6].

高溫脅迫可引起植物的生理、分子和生化變化，干擾細(xì)胞和整個(gè)植物的生長發(fā)育過程，從而對(duì)作物的生產(chǎn)產(chǎn)生不利影響[7].對(duì)玉米來說，除了在其他生長階段可能遭遇高溫脅迫外，在生殖生長及受精結(jié)實(shí)階段特別容易遇到高溫天氣，導(dǎo)致玉米減產(chǎn)[8].高溫對(duì)花粉(雄配子)發(fā)育的影響要大于花絲(雌配子)[9].在雌雄穗分化至籽粒成熟階段如果遭遇持續(xù)的高溫脅迫，可能造成敗育小花數(shù)量增加、花粉結(jié)構(gòu)破壞、雌穗結(jié)實(shí)率降低，最終導(dǎo)致減產(chǎn).研究發(fā)現(xiàn)，不同玉米材料耐高溫脅迫能力存在差異[10-11]，目前在國內(nèi)耐高溫的玉米材料較少，僅有鄭單958和中地88等較少的耐高溫玉米品種被推廣應(yīng)用.因此，以不同耐高溫脅迫能力的玉米花粉作為材料，研究在高溫環(huán)境下基因的差異表達(dá)具有重要意義.

基于此，本研究以耐高溫脅迫的中地88和高溫敏感的先玉335兩個(gè)玉米品種作為試驗(yàn)材料，在高溫脅迫環(huán)境下取其花粉作轉(zhuǎn)錄組測序分析，了解基因的差異表達(dá)情況，并進(jìn)一步對(duì)差異表達(dá)基因作功能注釋和代謝通路分析，以期為進(jìn)一步闡明玉米耐高溫脅迫的分子遺傳機(jī)制奠定理論基礎(chǔ)，為克隆應(yīng)用玉米耐高溫脅迫基因的分子育種提供參考.

1 材料與方法

1.1 試驗(yàn)材料

本試驗(yàn)選用了兩個(gè)玉米材料，一個(gè)是耐高溫脅迫的中地88，另一個(gè)是對(duì)高溫敏感的先玉335，兩個(gè)玉米材料均來自山西省農(nóng)業(yè)科學(xué)院玉米研究所.

1.2 試驗(yàn)方法

2019年將中地88和先玉335兩個(gè)玉米品種種植于山西省農(nóng)業(yè)科學(xué)院玉米研究所實(shí)驗(yàn)大棚，每個(gè)材料種植3行，行長3 m，每行10株.從播種到開始散粉前第5天，未作增溫處理，玉米生長環(huán)境包括溫度、濕度和光照等條件與田間自然生長的玉米材料一樣.在開始散粉前的第4天作增溫處理，但濕度和光照等條件與室外保持一致，在盛花期頭一天下午對(duì)雄穗進(jìn)行套袋，第2天12時(shí)30分時(shí)監(jiān)測到雄穗所處位置的溫度達(dá)到42 ℃，通過大棚通風(fēng)手段維持該溫度1 h，在13時(shí)30分，取兩個(gè)材料的花粉，各取4個(gè)重復(fù)，每個(gè)重復(fù)取約400 mg花粉裝進(jìn)2 mL離心管，立即液氮速凍，該過程在實(shí)驗(yàn)大棚內(nèi)完成.然后在實(shí)驗(yàn)室提取花粉總RNA，對(duì)符合要求的總RNA每個(gè)材料各選3份(作為3個(gè)生物學(xué)重復(fù))，中地88作為處理，用T1、T2和T3表示；先玉335作為對(duì)照，用CK1、CK2和CK3表示；再進(jìn)一步完成mRNA純化、反轉(zhuǎn)錄、建庫及測序等工作(該工作由上海派森諾生物科技股份有限公司完成).

1.3 數(shù)據(jù)統(tǒng)計(jì)分析

1.3.1 參考基因組選擇 試驗(yàn)結(jié)果分析選取玉米B73作為參考基因組，該參考物種在數(shù)據(jù)庫中注釋的基因數(shù)目及注釋比例為：NCBI_GI：16 607(41.94%)；KEGG：6 715(16.96%)；GO：30 113(76.06%)；EC：2 921(7.37%)；UniProtID：16 263(41.07%)；Ensembl：39 591(100%).

1.3.2 序列比對(duì)分析 對(duì)試驗(yàn)返回的結(jié)果，使用TopHaT2的升級(jí)版HISAT2(http://ccb.jhu.edu/software/hisaT2/index.shtml) 軟件分析，將過濾后的Reads比對(duì)到參考基因組上，并對(duì)獲得的基因在基因組上的區(qū)域分布進(jìn)行統(tǒng)計(jì).

1.3.3 基因表達(dá)量分析 使用HTSeq統(tǒng)計(jì)比對(duì)到每一個(gè)基因上的Read Count值作為基因的原始表達(dá)量.采用FPKM(Fragments Per Kilobase Million)對(duì)表達(dá)量進(jìn)行標(biāo)準(zhǔn)化，再繪制表達(dá)基因的FPKM密度分布的小提琴圖.

1.3.4 基因差異表達(dá)分析 根據(jù)基因在兩種玉米材料中的表達(dá)量，以先玉335材料為對(duì)照，以中地88材料為處理，分析上調(diào)基因、下調(diào)基因及無顯著差異表達(dá)基因的數(shù)量，并繪制差異表達(dá)基因的火山圖和在不同染色上的基因組圈圖.

對(duì)差異表達(dá)基因的生物學(xué)功能進(jìn)行GO(gene ontology)富集及功能注釋，對(duì)差異表達(dá)基因的GO富集分析結(jié)果，按照分子功能(molecular function，MF)、生物過程(biological process，BP) 和細(xì)胞組分(cellular component，CC) 進(jìn)行GO 分類，每個(gè)GO 分類中挑選P-value最小，即富集最顯著的前10個(gè)GO term條目進(jìn)行展示.此外，還對(duì)差異表達(dá)基因作KEGG(kyoto encyclopedia of genes and genomes)富集.根據(jù)差異表達(dá)基因的KEGG富集分析結(jié)果，挑選P-value最小即富集最顯著的前20個(gè)Pathway進(jìn)行展示.

2 結(jié)果與分析

2.1 Reads與參考基因組的比對(duì)

2.1.1 基本信息比對(duì) 使用TopHaT2的升級(jí)版HISAT2軟件將過濾后的 Reads 比對(duì)到參考基因組上.序列比對(duì)的基本信息統(tǒng)計(jì)結(jié)果見表1.對(duì)T1、T2、T3、CK1、CK2和CK3面言，比對(duì)上參考基因組的序列總數(shù)占比對(duì)的序列總數(shù)的百分比分別為92.50%、92.83%、92.24%、90.62%、91.25%和92.54%，平均值為91.96%；比對(duì)到多個(gè)位置的序列總數(shù)占比對(duì)上參考基因組的序列總數(shù)的百分比分別為5.88%、5.87%、6.29%、6.69%、6.47%和6.69%，平均值為6.34%；只比對(duì)到一個(gè)位置的序列總數(shù)占比對(duì)上參考基因組的序列總數(shù)分別為94.12%、94.13%、93.71%、93.31%、93.53%和93.31%，平均值為93.66%.

表1 Reads與參考基因組間的基本信息比對(duì)Table 1 Basic information comparison between reads and reference genome

2.1.2 比對(duì)區(qū)域分布 將比對(duì)到基因組上的Reads分布情況進(jìn)行統(tǒng)計(jì)，結(jié)果見表2.對(duì)T1、T2、T3、CK1、CK2和CK3面言，比對(duì)到基因區(qū)域的Reads總數(shù)占比對(duì)事件發(fā)生的總次數(shù)的百分比分別為94.24%、94.07%、93.53%、93.76%、94.10%和94.07%，平均值為93.96%；比對(duì)到基因間區(qū)的Reads總數(shù)占比對(duì)事件發(fā)生的總次數(shù)的百分比分別為5.76%、5.93%、6.47%、6.24%、5.90%和5.93%，平均值為6.04%；比對(duì)到外顯子區(qū)域的Reads總數(shù)占比對(duì)到基因區(qū)域的 Reads 總數(shù)的百分比分別為99.57%、99.59%、99.56%、99.45%、99.50%和99.45%，平均值為99.52%.

表2 比對(duì)到參考基因組上的Reads分布Table 2 Distribution of the Reads mapped to reference genome

2.2 基因表達(dá)的結(jié)果與分析

2.2.1 基因表達(dá)量 采用FPKM對(duì)表達(dá)量進(jìn)行標(biāo)準(zhǔn)化后，將表達(dá)量分為不同的區(qū)間，對(duì)各樣本在不同表達(dá)量區(qū)間內(nèi)基因的數(shù)目進(jìn)行統(tǒng)計(jì)，結(jié)果見表3.表達(dá)量值在0～0.01區(qū)間的基因數(shù)量最多，所有6個(gè)樣本都超過了20 000個(gè)，平均為21 228個(gè)，占總量的60.96%；隨表達(dá)量范圍值的提高，對(duì)應(yīng)基因數(shù)量逐漸減少，在>1000表達(dá)量區(qū)間所有6個(gè)樣本的基因數(shù)量都低于160個(gè)，平均為155.5個(gè)，僅占總量的0.45%.由此表得出，低表達(dá)的基因數(shù)量較多，而高表達(dá)的基因數(shù)量較少.

表3 不同表達(dá)量區(qū)間的基因數(shù)量統(tǒng)計(jì)Table 3 Statistics of the number for the genes in different expression quantity intervals

2.2.2 表達(dá)基因的FPKM密度分布 基因在各個(gè)樣品中的表達(dá)量利用小提琴圖進(jìn)行統(tǒng)計(jì)分析，結(jié)果如圖1所示.從圖1可以看出，這6個(gè)樣本具有相似的基因表達(dá)模式，即中等表達(dá)的基因占絕大多數(shù)，而低表達(dá)和高表達(dá)的基因僅占一小部分.

2.3 基因差異表達(dá)結(jié)果與分析

2.3.1 差異表達(dá)基因數(shù)量 中地88(Case:Traitment)相對(duì)于先玉335(Control：CK)材料而言，差異表達(dá)基因共計(jì)有1 822個(gè)，相比于Control，上調(diào)表達(dá)基因有857個(gè)，下調(diào)表達(dá)基因有965個(gè).另還有22 731個(gè)基因在這兩個(gè)玉米材料中表達(dá)差異不顯著.繪制的差異表達(dá)基因的火山圖見圖2.從圖2可看出，左側(cè)(藍(lán)色)為Case相比于Control下調(diào)基因，右側(cè)(紅色)為Case相比于Control上調(diào)基因，左右差異基因分布大致對(duì)稱.

2.3.2 差異表達(dá)基因的染色體分布 對(duì)在這兩個(gè)玉米材料中表達(dá)的基因在染色體上的分布作基因組圈圖，結(jié)果如圖3所示.從上調(diào)基因、下調(diào)基因及無顯著差異表達(dá)基因在整個(gè)玉米10條染色體上的分布來看，分布較均勻.

x為FPKM值；圖中盒型中間的橫線是中位數(shù)，盒型的上下邊緣為75%，上下限為90%，外部形狀為核密度估計(jì).x means FPKM value；The horizontal line in the middle of the box is the median，the upper and lower edges of the box are 75%，the upper and lower limits are 90%，and the outer shape is the kernel density estimation.圖1 FPKM密度分布Figure 1 Density distribution of FPKM

2.4 差異表達(dá)基因功能富集分析

2.4.1 GO富集分析 對(duì)差異表達(dá)基因進(jìn)行GO富集，結(jié)果顯示，與細(xì)胞組分(CC)相關(guān)的GO條目有491個(gè)，占總量的12.68%；與分子功能(MF)相關(guān)的GO條目有1 062個(gè)，占總量的27.42%；與生物過程(BP)相關(guān)的GO條目有2 320個(gè)，占總量的59.90%.對(duì)GO富集分析結(jié)果，每種分類挑選P-value最小，即富集最顯著的前10個(gè)GO條目進(jìn)行展示，結(jié)果見圖4.CC分類中，P-value值最小的基因的ID號(hào)為GO:0005887，與質(zhì)膜組成成分相關(guān)；MF分類中，P-value值最小的基因的ID號(hào)為GO:0004046，與氨?；D(zhuǎn)移酶活性相關(guān)；BP分類中，P-value值最小的基因的ID號(hào)為GO:0044281，與小分子代謝過程相關(guān).

x表示差異倍數(shù)；圖中兩條豎虛線為2倍表達(dá)差異閾值；橫虛線為P值的閾值(0.05).紅點(diǎn)表示該組上調(diào)基因，藍(lán)點(diǎn)表示該組下調(diào)基因，灰點(diǎn)表示非顯著差異表達(dá)基因.x means fold change；The two vertical dashed lines in the figure are 2 times of the expression difference threshold；the horizontal dashed line is the threshold value of P-value(0.05).The red，blue and gray dots represent the up-regulated，down-regulated and non-significantly differentially expressed genes，respectively.圖2 差異表達(dá)基因的火山圖Figure 2 Volcano map of differentially expressed genes

2.4.2 KEGG富集分析 參考KEGG通路系統(tǒng)的基因組信息數(shù)據(jù)庫，對(duì)檢測到的在這兩個(gè)玉米花粉材料中呈現(xiàn)差異表達(dá)的基因的代謝通路作富集分析.富集到的差異基因數(shù)目為349個(gè)，共涉及到102種通路，這些小通路可進(jìn)一步歸為6個(gè)大通路，包括細(xì)胞過程3個(gè)、環(huán)境信息處理4個(gè)、遺傳信息處理20個(gè)、人類疾病1個(gè)、新陳代謝72個(gè)和有機(jī)系統(tǒng)2個(gè).上調(diào)和下調(diào)表達(dá)基因分別有143和206個(gè)(注：一個(gè)基因可能同時(shí)富集于不同通路).挑選出的P-value最小即富集最顯著的前30個(gè)Pathway展示結(jié)果見圖5，其中1個(gè)與環(huán)境信息處理有關(guān)，4個(gè)與遺傳信息處理有關(guān)，25個(gè)與新陳代謝有關(guān).

最外圈是染色體條帶.紅色和綠色分別為上調(diào)和下調(diào)基因的log2x(Fold Change) 值的柱狀圖，灰色為無差異表達(dá)基因的log2x(Fold Change)值的散點(diǎn)圖.The outermost circle is the chromosomal banding.The red and green are the histogram of the log2x(Fold Change) value of up-regulated and down-regulated genes respectively，and the gray is the scatter graph of the log2x(Fold Change) value of the undifferentially expressed genes.圖3 基因組圈圖Figure 3 Genome circos

3 討論與結(jié)論

溫度是影響植物生長非常重要的因素[12-13]，高溫脅迫是最常見的非生物脅迫之一，它能顯著抑制植物生長發(fā)育[14].高溫脅迫一般會(huì)降低水分含量，影響植物光合作用和呼吸作用[15].研究發(fā)現(xiàn)，高溫脅迫已導(dǎo)致全球植物的產(chǎn)量下降和品質(zhì)降低[16-17].

1：質(zhì)膜組成成分;2：DNA指導(dǎo)的RNA聚合酶II，全酶；3：胞質(zhì)；4：核小體；5：胞質(zhì)大核糖體亞基；6：DNA包裝復(fù)合物；7：核DNA指導(dǎo)RNA聚合酶復(fù)合物；8：大核糖體亞基；9：MLL1/2復(fù)合物；10：MLL1復(fù)合物；11：氨基?；富钚?；12：蛋白異源二聚化活性；13：基本轉(zhuǎn)錄機(jī)制結(jié)合；14：基本RNA聚合酶II轉(zhuǎn)錄機(jī)制結(jié)合；15：肽基轉(zhuǎn)移酶活性；16：谷胱甘肽水解酶活性；17：輔因子結(jié)合；18：核苷跨膜轉(zhuǎn)運(yùn)體活性；19：磷酸吡哆醛結(jié)合；20：維生素B6結(jié)合；21：小分子代謝過程；22：赤霉素生物合成過程；23：應(yīng)答低氧含量；24：應(yīng)答氧水平；25：糖脂分解代謝過程；26：鞘糖脂分解代謝過程；27：萜類生物合成過程；28：磷代謝過程的正向調(diào)節(jié)過程；29：磷酸鹽代謝過程的正向調(diào)節(jié)過程；30：有機(jī)酸代謝過程.1：Integral component of plasma membrane；2：DNA-directed RNA polymerase II，homoenzyme；3：Cytosol；4：Nucleosome；5：Cytosolic large ribosomal subunit；6：DNA packaging complex；7：Nuclear DNA-directed RNA polymerase complex；8：Large ribosomal subunit；9：MLL1/2 complex；10：MLL1 complex；11：Aminoacylase activity；12：Protein heterodimerization activity；13：Basal transcription machinery binding；14：Basal RNA polymerase II transcription machinery binding；15：Peptidyltrasferase activity；16：Glutathione hydrolase activity；17：Cofactor binding；18：Nucleoside transmembrane transporter activity；19：Pyridoxal phosphate binding；20：Vitamin B6 binding；21：Small molecule metabolic process；22：Gibberellin biosynthetic process；23：Response to decreased oxygen levels；24：Response to oxygen levels；25：Glycolipid catabolic process；26：Glycosphingolipid catabolic process；27：Terpenoid biosynthetic process；28：Positive regulation of phosphorus metabolic process；29：Positive regulation of phosphate metabolic process；30：Organic acid metabolic process.圖4 差異表達(dá)基因的GO注釋和分類Figure 4 GO annotation and classification of differentially expressed genes

玉米是世界上最重要的農(nóng)作物之一，高溫脅迫同樣會(huì)顯著降低其產(chǎn)量和品質(zhì)[18].因此，研究玉米對(duì)高溫脅迫響應(yīng)的分子機(jī)制及選育耐高溫脅迫的玉米品種具有重要意義.

雖然已有關(guān)于玉米耐高溫脅迫的基因克隆或表達(dá)方面的研究，但是以前關(guān)于基因克隆的研究主要限于單個(gè)或少數(shù)幾個(gè)基因的克隆表達(dá)分析，如Wang等[19]從玉米中克隆了轉(zhuǎn)錄因子ZmWRKY106基因，將其轉(zhuǎn)入擬蘭芥，過量表達(dá)的結(jié)果顯示該基因會(huì)明顯提高對(duì)干旱和高溫脅迫的抗性；Ma等[20]研究發(fā)現(xiàn)ZmbZIP4基因的過表達(dá)可以調(diào)節(jié)ABA的生物合成及根的發(fā)育，從而提高玉米包括高溫脅迫等在內(nèi)的多種抗性；Li等[21]也發(fā)現(xiàn)玉米ZmbZIP60基因的表達(dá)會(huì)應(yīng)答于高溫脅迫.由于玉米對(duì)高溫脅迫的抗逆性屬多基因控制的數(shù)量性狀，因而對(duì)單個(gè)基因的克隆表達(dá)研究還不能完全滿足育種需要，有必要在基于轉(zhuǎn)錄組測序的在全基因組水平下研究耐高溫脅迫基因.

基于轉(zhuǎn)錄組測序?qū)τ衩啄透邷孛{迫的基因差異表達(dá)以前也有類似研究，且發(fā)現(xiàn)了大量與高溫脅迫相關(guān)的差異表達(dá)基因，但以前的研究取材部位主要限于玉米葉片[1，18，22-23].然而在包括中國在內(nèi)的許多國家，在玉米散粉期常遭遇高溫天氣，高溫脅迫導(dǎo)致花粉活力下降，從而影響玉米的授粉、結(jié)實(shí)，使其大幅減產(chǎn)，因此，選擇以花粉為材料，基于轉(zhuǎn)錄組高通量測序，研究不同玉米材料間的基因差異表達(dá)對(duì)于認(rèn)識(shí)玉米耐高溫脅迫的分子遺傳機(jī)制具有重要意義.

本試驗(yàn)在研究方法、材料選擇或取材部位等方面，不同于以前的研究[19-23]，本試驗(yàn)是以中地88(耐高溫脅迫)和先玉335(高溫敏感)兩個(gè)玉米品種的花粉為研究對(duì)象，研究其基因表達(dá)情況，共發(fā)現(xiàn)差異表達(dá)基因1 822個(gè)，相對(duì)于對(duì)照(先玉335)而言，中地88上調(diào)表達(dá)基因有857個(gè)，下調(diào)表達(dá)基因965個(gè).另外還發(fā)現(xiàn)有22 731個(gè)基因在兩個(gè)材料中表達(dá)差異不顯著.對(duì)差異表達(dá)基因進(jìn)行GO功能注釋，結(jié)果顯示大部分基因都與細(xì)胞組分、分子功能或生物過程功能相關(guān).KEGG功能富集結(jié)果顯示，這些差異表達(dá)基因可以歸為6個(gè)Level 1分類(包括102個(gè)Levle 2通路).綜合分析這些差異表達(dá)基因，發(fā)現(xiàn)有7個(gè)基因直接應(yīng)答于熱激效應(yīng)，其中2個(gè)上調(diào)表達(dá)，5個(gè)下調(diào)表達(dá)；有3個(gè)基因與熱激蛋白結(jié)合有關(guān)，其中上調(diào)表達(dá)基因1個(gè)，下調(diào)表達(dá)基因2個(gè).在這12個(gè)與玉米耐高溫脅迫相關(guān)的基因中，多數(shù)為在玉米中新發(fā)現(xiàn)的與耐高溫脅迫相關(guān)的基因.依據(jù)基因的ID號(hào)和玉米參考基因組序列信息，參考以前的研究[24-26]，進(jìn)一步對(duì)其基因全長序列、cDNA序列及表達(dá)蛋白的生物學(xué)特性等作進(jìn)一步研究，關(guān)于這些基因的分子結(jié)構(gòu)及遺傳功能研究正在進(jìn)行之中.本實(shí)驗(yàn)的研究結(jié)果為進(jìn)一步闡明玉米耐高溫脅迫的分子遺傳機(jī)制奠定了一些理論基礎(chǔ)，為在植物中克隆應(yīng)用這些耐高溫脅迫基因提供了參考.

1：植物MAPK信號(hào)通路；2：核糖體；3：堿基切除修復(fù)；4：硫傳遞系統(tǒng)；5：RNA運(yùn)輸；6：其他多糖降解；7：糖酵解/糖異生；8：α-亞麻酸代謝；9：淀粉和蔗糖代謝；10：苯并噁唑嗪酮類化合物生物合成；11：鞘脂代謝；12：丙酮酸代謝；13：丁酸代謝；14：脂肪酸降解；15：磷酸戊糖途徑；16：氰氨基酸代謝；17：抗壞血酸和醛糖酸代謝；18.苯丙氨酸代謝；19：咖啡因代謝；20：纈氨酸、亮氨酸和異亮氨酸降解；21：組氨酸代謝；22：苯丙氨酸、酪氨酸、色氨酸生物合成；23：果糖和甘露糖代謝；24：亞油酸代謝；25：?；撬崤c次?；撬岽x；26：煙酸和煙酰胺代謝；27：半乳糖代謝；28：硫胺代謝；29：有機(jī)含硒化合物代謝；30：泛醌和其他萜醌生物合成1：MAPK signaling pathway - plant；2：Ribosome；3：Base excision repair；4：Sulfur relay system；5：RNA transport；6：Other glycan degradation；7：Glycolysis / Gluconeogenesis；8：alpha-Linolenic acid metabolism；9：Starch and sucrose metabolism；10：Benzoxazinoid biosynthesis；11：Sphingolipid metabolism；12：Pyruvate metabolism；13：Butanoate metabolism；14：Fatty acid degradation；15：Pentose phosphate pathway；16：Cyanoamino acid metabolism；17：Ascorbate and aldarate metabolism；18：Phenylalanine metabolism；19：Caffeine metabolism；20：Valine，leucine and isoleucine degradation；21：Histidine metabolism；22：Phenylalanine，tyrosine and tryptophan biosynthesis；23：Fructose and mannose metabolism；24：Linoleic acid metabolism；25：Taurine and hypotaurine metabolism；26：Nicotinate and nicotinamide metabolism；27：Galactose metabolism；28：Thiamine metabolism；29：Selenocompound metabolism；30：Ubiquinone and other terpenoid-quinone biosynthesis圖5 KEGG通路富集結(jié)果柱狀圖Figure 5 Histogram of enrichment results for KEGG pathway