巨菌草幼葉及根轉錄組功能基因測序及分析

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

（1.福建農林大學國家菌草工程技術研究中心，福建福州350002；2.福建農林大學動物科學學院，蜂學學院，福建 福州350002）

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

近年來，分子生物學發展極為迅猛，這就為揭示植物生長發育潛在調控機制提供了便利。通過采用RNA 高通量深度測序技術（RNA-seq），在分子水平上挖掘植物的功能基因和轉錄本變化，以及植物發育的生理生化代謝，進而闡釋植物生長發育的響應過程［4］。與以往測試方法相比，RNA-seq 具有價格低廉，靈敏度強，通量高等優點［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］等重要 牧草中，均 有應用RNA?seq 技術，進行植物生長發育、抗性等研究的相關報道。

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

1 材料與方法

1.1 試驗材料

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

1.2 RNA 提取及轉錄組測序

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

1.3 De novo 組裝及質量控制

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

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

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

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

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

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

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

表1 巨菌草幼葉、根8 個DEGs 和內參基因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 數據處理與分析

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

2 結果與分析

2.1 巨菌草轉錄組測序與De nove 組裝

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

2.2 Unigene 注釋結果

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

表2 巨菌草轉錄組組裝結果Table 2 Summary of transcriptome assembly for Giant Jun?cao（bp）

圖1 基因注釋結果統計（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 個亞類，包括生物學過程，細胞組成和分子功能（圖2）。其中生物學過程中以細胞進程（cellular process）和代謝進程（metabolic process）富集unige?ne 個數較多；細胞組成中則主要以細胞（cell）和細胞成分（cell part）富集較多的unigene；結合（binding）和催化活性（catalytic activity）在分子功能中富集到的unigene 個數較多。巨菌草生長相關的GO 富集分析顯示，與發育生長調節（regulation of developmental growth）相關的unigene 個數最多，為108 個，其次與涉及形態發生的生長發育（developmental growth involved in morphogenesis）相關的unigene 有67 個，與細胞生長發育（developmental cell growth）相關的unigene 有52 個（表3）。表明巨菌草調控生物量的基因主要分布于以上幾個生物學過程中。

2.4 KOG 功能分析

圖2 GO 功能分布統計Fig.2 Functional distribution of GO annotation

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

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

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

2.5 差異表達基因分析

篩選巨菌草幼葉與根差異表達基因（圖4），共鑒定出5735 個差異表達基因，其中上調基因有3435 個，下調基因2300 個，分別占總DEGs 的59.90% 和40.10%，得到的差異基因將進行后續GO 功能，KEGG 代謝通路及轉錄因子分析。

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

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

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

表4 巨菌草幼葉、根上調差異基因GO 富集結果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 個數分別為52，25，396，35 和109 個（表6）。

2.7 差異表達基因轉錄因子分析

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

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

續表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 個基因在幼葉和根中的表達程度不同，但表達趨勢與高通量測序結果基本一致，表明測序結果真實可信。

3 討論

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

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

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

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

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

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

轉錄調控是植物生長發育的一個重要環節，是通過對其轉錄靶基因的臨時和空間調控來實現的［22?23］。本研究顯示，bHLH 和WRKY 是巨菌草幼葉、根DEGs 最主要的兩大轉錄因子家族。bHLH 和WRKY 在光信號、激素信號轉導、創傷和干旱脅迫反應中發揮重要作用，是植物中最常見的TFs［24?25］。已有研究發現，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 個基因屬于該TFs 家族［32］。不同物種TFs 個數不同，這與物種所處環境，生理狀況，發育階段有很大的關系，對轉錄因子的研究也可以為植物生長發育調控機制提供寶貴基因。

4 結論

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