周晶，陳思齊，史文嬌，陽伏林，林輝，林占熺

（1.福建農(nóng)林大學(xué)國家菌草工程技術(shù)研究中心，福建福州350002；2.福建農(nóng)林大學(xué)動物科學(xué)學(xué)院，蜂學(xué)學(xué)院，福建 福州350002）

巨菌草（Pennisetum giganteum）屬禾本科狼尾草屬多年生C4植物，是一種適宜在我國熱帶、亞熱帶、溫帶地區(qū)生長和人工栽培的高產(chǎn)優(yōu)質(zhì)飼草［1］。巨菌草最重要的兩大特色，一是其生物量巨大，在適宜生長條件下，年鮮草產(chǎn)量可達200 t·hm?2以上，是當今世界生物量最大的飼草植物；二是其根系發(fā)達，尤以須根為主，對土壤適應(yīng)能力強，是環(huán)境脆弱地區(qū)生態(tài)治理的理想草種［2］。目前已在砒砂巖地（內(nèi)蒙古鄂爾多斯），灘涂鹽堿地（福建平潭），洪積扇（青海貴德）等地成功試種，且生物量仍顯著高于當?shù)貎?yōu)勢草種。然而由于巨菌草為野生栽培品種，其種質(zhì)資源遺傳背景模糊，對其分子生物學(xué)理論研究相對較薄弱，這也是限制巨菌草品種改良的主要因素［3］。

近年來，分子生物學(xué)發(fā)展極為迅猛，這就為揭示植物生長發(fā)育潛在調(diào)控機制提供了便利。通過采用RNA 高通量深度測序技術(shù)（RNA-seq），在分子水平上挖掘植物的功能基因和轉(zhuǎn)錄本變化，以及植物發(fā)育的生理生化代謝，進而闡釋植物生長發(fā)育的響應(yīng)過程［4］。與以往測試方法相比，RNA-seq 具有價格低廉，靈敏度強，通量高等優(yōu)點［5］，目前在小麥（Triticum aestivum）［6］、玉米（Zea mays）［7］、水稻（Oryza sativa）［8］、高粱（Sorghum bicolor）［9］、粟（Setaria italica）［10］等糧食作物，及紫花苜蓿（Medicago sativa）［11］，多年生黑麥草（Lolium perenne）［12］，霸王（Zy?gophyllum xanthoxylum）［13］，羊草（Leymus chinensis）［14］，草地早熟禾（Poa pratensis）［15］等重要 牧草中，均 有應(yīng)用RNA?seq 技術(shù)，進行植物生長發(fā)育、抗性等研究的相關(guān)報道。

針對現(xiàn)有研究中巨菌草轉(zhuǎn)錄組信息缺乏的情況，本研究采用Illumina 高通量測序技術(shù)對巨菌草幼葉和根進行轉(zhuǎn)錄組測序，結(jié)合生物信息分析對獲得的單基因簇（unigenes）和差異表達基因（differentially expressed genes，DEGs）開展功能注釋、代謝通路、轉(zhuǎn)錄因子等方面的分析研究。本研究極大地豐富了巨菌草的基因資源，為進一步開展巨菌草功能基因挖掘，分子育種等方面研究提供了數(shù)據(jù)和理論依據(jù)。

1 材料與方法

1.1 試驗材料

巨菌草于2018 年7 月采自福建農(nóng)林大學(xué)國家菌草工程技術(shù)研究中心菌草培育基地。以栽培時間一致，同時具有7 片葉的幼嫩植株作為試驗材料，分別取其幼葉和根，每個材料選取3 棵健康單株作為3 個重復(fù)，分別剪取葉片和根后立即經(jīng)液氮速凍，保存于?80 ℃冰箱備用，用于RNA 提取。

1.2 RNA 提取及轉(zhuǎn)錄組測序

分別取巨菌草幼葉、根各1 g 材料，根據(jù)Trizol 試劑盒說明書，按步驟提取總RNA。利用NanoDrop 1000 分光光度計（Thermo scientific 公司，美國）和Agilent 2100 分析儀（Agilent Technologies 公司，美國）對提取的總RNA純度和完整性進行檢測。之后將各樣品RNA 等量混合，用以構(gòu)建轉(zhuǎn)錄組文庫。對構(gòu)建好的文庫采用Illumina HiSeq 4000 進行雙端測序（paired-end）分析，讀長為150 bp。

1.3 De novo 組裝及質(zhì)量控制

首先將預(yù)處理測序得到的原始數(shù)據(jù)（raw reads），通過移除接頭、poly-N 和低質(zhì)量片段，獲得高質(zhì)量的干凈數(shù)據(jù)（clean data）。其次計算clean data 的GC 含量、Q20、Q30 以及序列重復(fù)水平。最后利用Trinity 方法對clean reads 進行無參轉(zhuǎn)錄組拼接，構(gòu)建轉(zhuǎn)錄本（transcripts）和單基因簇（unigene）序列［16］。

1.4 Unigene 序列注釋、分類及功能注釋

利用BLAST 軟件將巨菌草unigene 序列與7 個最常見的公共數(shù)據(jù)庫進行比對，并根據(jù)基因的相似性進行功能注釋［17］。采用的數(shù)據(jù)庫有非冗余蛋白質(zhì)序列數(shù)據(jù)庫（non-redundant protein database，Nr）（https：//www.ncbi.nlm.nih.gov/）；核 酸 序 列 數(shù) 據(jù) 庫（non-redundant nucleotide sequences，Nt）（http：//ftp.ncbi.nlm.nih.gov/BLAST/db/）；蛋白質(zhì)序列數(shù)據(jù)庫（swiss prot protein database，Swiss-Prot）（https：//web.expasy.org/docs/Swis?sprot）；真 核 直 系 同 源 基 因 數(shù) 據(jù) 庫（eukaryotic orthologous groups，KOG）（http：//www.ncbi.nlm.nih.gov/COG/）；基因本體數(shù)據(jù)庫（gene ontology，GO）http：//www.gene ontology.org/）；蛋白質(zhì)家族數(shù)據(jù)庫（protein fami?lies database，Pfam）（http：//Pfam.sanger.ac.uk/）；京都基因和基因組途徑數(shù)據(jù)庫百科全書（kyoto encyclopedia of genes and genomes，KEGG）（http：//www.genome.jp/KEGG/）。

1.5 差異表達基因鑒定及分析

利用每百萬reads 中來自某基因每千堿基長度的reads 數(shù)（reads per kilobases per million mapped reads，RP?KM）方法計算所有鑒定到的unigene 的表達量。使用TMM 方法進行差異基因分析，同時結(jié)合錯誤發(fā)現(xiàn)率（false discovery rate，F(xiàn)DR）（FDR<0.05）和基因表達量倍數(shù)變化（|log2 fold change|>1）這兩個指標，獲得巨菌草幼葉、根DEGs。對獲得的DEGs 進行GO 功能顯著性富集和KEGG 途徑分析，從而獲得DEGs 顯著相關(guān)的生物學(xué)功能和代謝通路［18?19］。同時將所有鑒定到的DEGs 進行轉(zhuǎn)錄因子分析，將DEGs 比對到植物轉(zhuǎn)錄因子數(shù)據(jù)庫4.0（plant transcription factor database，PlantTFDB 4.0）（http：//planttfdb.cbi.pku.edu.cn/），所設(shè)閾值為1×e?5，最終獲得巨菌草葉、根表達轉(zhuǎn)錄因子分析結(jié)果。

1.6 差異表達基因qRT?PCR 驗證

為了驗證高通量數(shù)據(jù)，隨機選取8 個DEGs 進行qRT?PCR 分析。分別提取巨菌草葉片與根的RNA，采用FastKing gDNA Dispelling RT SuperMix 試劑盒反轉(zhuǎn)錄成cDNA，設(shè)計引物（表1），并以PgACT為內(nèi)參基因，利用CFX Connect 熒光定量PCR 儀（BIO?RAD 公司，美國）進行qRT?PCR，每個樣品測試3 次，總反應(yīng)體系為20 μL，利用公式2???Ct計算其相對表達量。

表1 巨菌草幼葉、根8 個DEGs 和內(nèi)參基因PgACT 的qRT-PCR 引物Table 1 Primers of eight DEGs and reference gene PgACT used for quantitative real-time PCR in leaves and roots of Giant Juncao

1.7 數(shù)據(jù)處理與分析

利用Microsoft Office Excel 2016 軟件進行DEGs 相對表達量水平的統(tǒng)計分析，并使用SigmaPlot 14.0 軟件和R 語言作圖。

2 結(jié)果與分析

2.1 巨菌草轉(zhuǎn)錄組測序與De nove 組裝

通過采用Illumina Hiseq 4000 高通量測序技術(shù)，對巨菌草幼葉、根組織進行轉(zhuǎn)錄組測序，共有144670426 條raw reads 和141886800 條clean reads，總長度分別為21.68 和21.28 Gb。對所有樣品的轉(zhuǎn)錄組clean reads 進行無參組裝，共得到210806 條轉(zhuǎn)錄本和150336 條unigene。其中最短transcripts 為201 bp，最長為14723 bp，平均為760 bp；最短unigene 為201 bp，最長為14723 bp，平均unigene 為642 bp。巨 菌 草transcripts 和unigene 長度分布如表2 所示。

2.2 Unigene 注釋結(jié)果

為了預(yù)測組裝到的unigene 的功能，將其與7 個主要生物學(xué)數(shù)據(jù)庫（Nr，Nt，KEGG，Swiss-Prot，Pfam，GO 和KOG）進行比對。7 個數(shù)據(jù)庫均可注釋到的unigene 有9671 個，至少注釋到一個數(shù)據(jù)庫的unigene有88765 個，其余各個數(shù)據(jù)庫分別注釋到的unigene 有76043，50882，24322，44054，54032，54440 和31705 個（圖1a）。近緣物種相似序列匹配分析結(jié)果顯示，巨菌草與粟相似度最高，達41.53%；其次為高粱，相似度為6.49%，表明巨菌草與狗尾草屬粟有較高的近緣關(guān)系（圖1b）。

表2 巨菌草轉(zhuǎn)錄組組裝結(jié)果Table 2 Summary of transcriptome assembly for Giant Jun?cao（bp）

圖1 基因注釋結(jié)果統(tǒng)計（a）及Nr 注釋的物種分布（b）Fig.1 Unigene information annotated in different databases（a）and distributed into different species from Nr database（b）

2.3 GO 功能分析

將注釋到的unigene 進行GO 分析，顯示unigene 主要被分為3 個大類和41 個亞類，包括生物學(xué)過程，細胞組成和分子功能（圖2）。其中生物學(xué)過程中以細胞進程（cellular process）和代謝進程（metabolic process）富集unige?ne 個數(shù)較多；細胞組成中則主要以細胞（cell）和細胞成分（cell part）富集較多的unigene；結(jié)合（binding）和催化活性（catalytic activity）在分子功能中富集到的unigene 個數(shù)較多。巨菌草生長相關(guān)的GO 富集分析顯示，與發(fā)育生長調(diào)節(jié)（regulation of developmental growth）相關(guān)的unigene 個數(shù)最多，為108 個，其次與涉及形態(tài)發(fā)生的生長發(fā)育（developmental growth involved in morphogenesis）相關(guān)的unigene 有67 個，與細胞生長發(fā)育（developmental cell growth）相關(guān)的unigene 有52 個（表3）。表明巨菌草調(diào)控生物量的基因主要分布于以上幾個生物學(xué)過程中。

2.4 KOG 功能分析

圖2 GO 功能分布統(tǒng)計Fig.2 Functional distribution of GO annotation

表3 與巨菌草生長相關(guān)的GO 富集分析Table 3 GO enrichment analysis related to growth annotation in Giant Juncao

將巨菌草unigene 與KOG 數(shù)據(jù)庫進行比對，并對比對到的unigene 進行分類統(tǒng)計和功能預(yù)測。研究結(jié)果表明，巨菌草unigene 按其功能大致可以分為26 類，并對每一類unigene 個數(shù)進行統(tǒng)計（圖3）。其中涉及unigene 數(shù)目最多的是一般功能預(yù)測（general function prediction only），為7352 個；其次是翻譯后修飾，蛋白質(zhì)轉(zhuǎn)換，分子伴侶（posttranslational modification，protein turnover，chaperones），有3685 個unigene。涉 及 影 響巨 菌 草 生 物量 的unigene 主要分布于能量產(chǎn)生和轉(zhuǎn)換（energy production and conversion），碳水化合物運輸和代謝（carbohydrate transport and metabolism），次生代謝產(chǎn)物的生物合成、運輸和分解代謝（secondary metabolites biosynthesis，trans?port and catabolism）和脂質(zhì)轉(zhuǎn)運與代謝（lipid transport and metabolism），其 個 數(shù) 分 別 為2050，1889，1668 和1405 個。

圖3 COG/KOG 功能分布統(tǒng)計Fig.3 COG/KOG function classification of unigenes involved in drought tolerance in the Giant Juncao

2.5 差異表達基因分析

篩選巨菌草幼葉與根差異表達基因（圖4），共鑒定出5735 個差異表達基因，其中上調(diào)基因有3435 個，下調(diào)基因2300 個，分別占總DEGs 的59.90% 和40.10%，得到的差異基因?qū)⑦M行后續(xù)GO 功能，KEGG 代謝通路及轉(zhuǎn)錄因子分析。

2.6 差異表達基因GO 功能和KEGG 通路富集分析

圖4 巨菌草幼葉、根差異基因表達分析火山圖Fig.4 Volcanic map of DEGs analysis between leaves and roots of Giant Juncao

分別將上調(diào)和下調(diào)的差異基因進行GO 功能富集分析，將FDR<0.05 作為顯著富集選項，共鑒定到上調(diào)基因GO 顯著富集結(jié)果20 條，下調(diào)基因顯著富集結(jié)果53 條。其中上調(diào)的DEGs 顯著富集的前5 個功能主要包括有銅離子結(jié)合（copper ion binding），化學(xué)刺激反應(yīng)（response to chemical stimulus），抗氧化活性（an?tioxidant activity），應(yīng)對水（response to water），過氧化物酶活性（peroxidase activity）等功能上，其DEGs 個數(shù)分別為35，73，42，25 和38。DEGs 富集較多的功能主要在電子傳遞（electron transport），碳水化合物代謝過程（carbohydrate metabolic process），離子結(jié)合（ion binding），陽離子綁定（cation binding），金屬離子結(jié)合（metal ion binding），每個功能都有超過150 個DEGs 富集（表4）。下調(diào)的DEGs 主要集中在生物學(xué)過程和分子功能方面，顯著富集的前5 項主要包括有光合作用、光收獲（photosynthesis，light harvesting），光合作用、光反應(yīng)（photosynthe?sis，light reaction），電子傳遞（electron transport），光合作用（photosynthesis）和氧化還原酶活性（oxidoreductase ac?tivity）。DEGs 富集個數(shù)較多的功能是代謝過程（metabolic process）和催化活性（catalytic activity），分別為883 和767（表5）。

表4 巨菌草幼葉、根上調(diào)差異基因GO 富集結(jié)果Table 4 Giant Juncao GO enrichment of leaves and roots upregulated DEGs（FDR<0.05）

對巨菌草幼葉、根DEGs 進行KEGG 富集分析，共鑒定到120 條代謝通路。將FDR<0.05 作為顯著富集指標，共篩選出光合生物的固碳作用（carbon fixation in photosynthetic organisms），類胡蘿卜素生物合成（carotenoid biosynthesis），次生代謝物的生物合成（biosynthesis of secondary metabolites），磷酸戊糖途徑（pentose phosphate pathway）和碳代謝（carbon metabolism）5 條代謝途徑，且每條代謝途徑包含的DEGs 個數(shù)分別為52，25，396，35 和109 個（表6）。

2.7 差異表達基因轉(zhuǎn)錄因子分析

對巨菌草幼葉、根DEGs 進行轉(zhuǎn)錄因子分析，結(jié)果顯示，共鑒定到3100 個DEGs，隸屬于56 個轉(zhuǎn)錄因子家族，其中bHLH 轉(zhuǎn)錄因子所含的DEGs 數(shù)量最高，為316 個；其次是WRKY 轉(zhuǎn)錄因子，有220 個DEGs 注釋到其中；剩余 轉(zhuǎn) 錄 因 子 中，DEGs 個 數(shù) 超 過100 的 分 別 為MYB-related，B3，NAC，bZIP，F(xiàn)AR1，ERF，GRAS，C3H 和C2H2（圖5）。

表5 巨菌草幼葉、根下調(diào)差異基因GO 富集結(jié)果Table 5 Giant Juncao GO enrichment of leaves and roots downregulated DEGs（FDR<0.05）

續(xù)表5 Continued Table 5

表6 巨菌草幼葉、根差異基因KEGG 富集分析Table 6 Giant Juncao KEGG enrichment of leaves and roots DEGs（FDR＜0.05）

2.8 qRT?PCR 熒光定量驗證

對隨機選取的8 個DEGs 進行qRT?PCR 驗證，由圖6 可知，8 個基因在幼葉和根中的表達程度不同，但表達趨勢與高通量測序結(jié)果基本一致，表明測序結(jié)果真實可信。

3 討論

圖5 巨菌草幼葉、根差異表達基因轉(zhuǎn)錄因子Fig.5 Transcription factor family of leaves and roots differential expression genes in Giant Juncao

巨菌草由最初的以替代林木作為培養(yǎng)基栽培食藥用菌，逐步發(fā)展成為具有可供家畜采食、修復(fù)脆弱生態(tài)環(huán)境等特性的綜合型草種，目前已在我國三十多個省被廣泛種植［3］。巨菌草巨大地上生物量及發(fā)達根系正是其被不斷開發(fā)利用的主要原因，因此揭示影響巨菌草生物量的主要因子就顯得尤為重要。

本試驗利用轉(zhuǎn)錄組測序方法，將巨菌草以根為代表的地下生物量和以葉為代表的地上生物量進行比較，找到二者差異表達基因及其基因富集情況?；贒e novo組裝，共獲得150336 條unigenes，平均長度為642 bp。該unigene 個數(shù)要小于草地早熟禾［15］和羊草unigene 個數(shù)（254331 和180770）［20］，但是多于海濱雀稗（Paspalum vagi?natum）unigene 個數(shù)（117619）［21］，這可能是由于不同物種自身攜帶的遺傳物質(zhì)堿基序列不同，因此在轉(zhuǎn)錄組組裝過程中會出現(xiàn)較大差異。本研究中，由于巨菌草尚未完成全基因組測序，因此只有88765 個unigene 被注釋到至少一個數(shù)據(jù)庫中，說明還有很多未知功能的基因可能是控制巨菌草巨大生物量的主要原因，對其進一步挖掘還有待于全基因組序列的測序結(jié)果。通過對巨菌草進行Nr 數(shù)據(jù)庫物種分布比對，認為巨菌草與狗尾草屬粟的親緣關(guān)系最近，相似度高達41.53%。粟的全基因組序列測序已經(jīng)完成，這也為本研究提供了一種思路，通過利用近緣種粟來挖掘與巨菌草生物量相關(guān)的調(diào)控基因?qū)⒊蔀橐环N可能。

圖6 隨機挑選巨菌草8 個DEGs 的相對表達量水平Fig.6 Relative expression levels of eight randomly selected DEGs of Giant Juncao

巨菌草幼葉、根差異表達基因共有5735 個，其中包含3435 個上調(diào)DEGs 和2300 個下調(diào)DEGs。上調(diào)的DEGs顯著富集在GO 的電子傳遞（202 個unigenes），碳水化合物代謝（191 個unigenes），金屬離子結(jié)合（152 個unige?nes），陽離子綁定（152 個unigenes）和離子結(jié)合（152 個unigenes）功能中；下調(diào)的DEGs 主要富集在光合作用、光收獲（21 個unigenes）、光合作用、光反應(yīng)（21 個unigenes）等過程中。本研究中有大量DEGs 在離子、電子傳遞結(jié)合過程中起主要調(diào)控作用，推測可能是由植物根呼吸作用引起的。根的呼吸過程與對礦物質(zhì)元素的吸收及分化密不可分，因此導(dǎo)致了大量帶電荷物質(zhì)的轉(zhuǎn)運。這也反映了巨菌草發(fā)達根系是由復(fù)雜的生理生化過程調(diào)控的，要了解其原因也不是通過單個基因可以完成的。巨菌草葉在光合作用和光反應(yīng)方面較根具有更顯著的差異表達水平，反映了葉通過光合作用在能量積累、物質(zhì)貯存等方面具有更強的能力。通常認為，巨菌草是C4植物，其在對光能利用方面要較C3植物更具優(yōu)勢，生物量也因此更高。

對差異基因的KEGG 富集分析顯示，葉與根在光合生物固碳作用，類胡蘿卜素生物合成，次生代謝物的生物合成，磷酸戊糖途徑和碳代謝途徑中具有顯著的差異性。以往研究表明，不同處理條件下，植物差異表達基因富集通路通常不相同。對羊草轉(zhuǎn)錄組結(jié)果分析顯示，經(jīng)鹽堿脅迫處理過的羊草，其差異表達基因主要富集在脅迫忍受功能，信號轉(zhuǎn)導(dǎo)，能量生產(chǎn)和轉(zhuǎn)化，無機離子傳輸?shù)韧罚?4］。在干旱脅迫處理過的草地早熟禾中，蛋白激酶，蛋白磷酸酶，碳代謝和植物激素信號等通路，是差異表達基因的主要富集結(jié)果［15］。而多年生黑麥草經(jīng)冷凍脅迫之后，大量的差異表達基因則富集在信號轉(zhuǎn)導(dǎo)，ABA 刺激響應(yīng)，代謝過程等通路［12］。由此可見，不同處理下，植物轉(zhuǎn)錄組經(jīng)信息分析后得到的差異表達基因通路差別較大，這也體現(xiàn)了轉(zhuǎn)錄組對不同生長條件下植物生長響應(yīng)的差別反應(yīng)。本研究中，差異基因的富集通路主要與光合作用相關(guān)，表明巨菌草在光合作用方面具有明顯優(yōu)勢，這也是其生物量巨大的一個主要因素。

轉(zhuǎn)錄調(diào)控是植物生長發(fā)育的一個重要環(huán)節(jié)，是通過對其轉(zhuǎn)錄靶基因的臨時和空間調(diào)控來實現(xiàn)的［22?23］。本研究顯示，bHLH 和WRKY 是巨菌草幼葉、根DEGs 最主要的兩大轉(zhuǎn)錄因子家族。bHLH 和WRKY 在光信號、激素信號轉(zhuǎn)導(dǎo)、創(chuàng)傷和干旱脅迫反應(yīng)中發(fā)揮重要作用，是植物中最常見的TFs［24?25］。已有研究發(fā)現(xiàn)，bHLH 家族在水稻、玉米和小麥［26］基因組中分別有183、231 和571 個成員，在擬南芥（Arabidopsis thaliana）［27］、白菜（Brassica rapassp.）［28］、二穗短柄草（Brachypodium distachyon）［29］和蘋果（Malus domestica）［30］中則分別有162、230、146 和175 個成員。對WRKY 的研究顯示，水稻中有103 個基因［31］，玉米和高粱中則分別為116 和68 個基因?qū)儆谠揟Fs 家族［32］。不同物種TFs 個數(shù)不同，這與物種所處環(huán)境，生理狀況，發(fā)育階段有很大的關(guān)系，對轉(zhuǎn)錄因子的研究也可以為植物生長發(fā)育調(diào)控機制提供寶貴基因。

4 結(jié)論

本研究以正常生長的巨菌草為材料，以其地上部分葉和地下部分根為比較進行轉(zhuǎn)錄組分析，共獲得150336條轉(zhuǎn)錄unigenes，物種分布鑒定到巨菌草與粟的親緣關(guān)系更近。地下和地上生物量比較共鑒定到3435 個上調(diào)差異基因和2300 個下調(diào)差異基因，GO 分析的上調(diào)DEGs 主要富集在電子傳遞，碳水化合物代謝過程，離子結(jié)合等功能；下調(diào)DEGs 主要富集在代謝過程和催化活性等過程。KEGG pathway 分析則顯示參與光合生物的固碳作用代謝通路差異基因富集程度最高；參與次生代謝物的生物合成代謝通路的DEGs 數(shù)目最多。bHLH 和WRKY 是巨菌草幼葉、根DEGs 最主要的轉(zhuǎn)錄因子家族。本研究極大地豐富了巨菌草不同組織中轉(zhuǎn)錄組信息，為今后其分子生物學(xué)研究提供了寶貴數(shù)據(jù)基礎(chǔ)，同時也為其品種改良提供了參考。