巨菌草幼叶及根转录组功能基因测序及分析

周晶，陈思齐，史文娇，阳伏林，林辉，林占熺

（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］。巨菌草巨大地上生物量及发达根系正是其被不断开发利用的主要原因，因此揭示影响巨菌草生物量的主要因子就显得尤为重要。

本试验利用转录组测序方法，将巨菌草以根为代表的地下生物量和以叶为代表的地上生物量进行比较，找到二者差异表达基因及其基因富集情况。基于De 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 最主要的转录因子家族。本研究极大地丰富了巨菌草不同组织中转录组信息，为今后其分子生物学研究提供了宝贵数据基础，同时也为其品种改良提供了参考。