基于转录组测序揭示玉米抗倒伏相关基因和代谢通路

刘松涛 田再民 刘子刚 高志佳 张静 贺东刚 黄智鸿 兰鑫

（1 河北北方学院/河北省农产品食品质量安全分析检测重点实验室，075000，河北张家口；2 河北巡天农业科技有限公司，075000，河北张家口；3 河北兆育种业集团有限公司，050000，河北石家庄）

玉米是我国重要的粮食作物之一，具有高产、适应性广和营养物质丰富等特点，因此玉米高产稳产在农业生产中占有重要地位[1]。玉米生长发育过程中经常发生倒伏现象，造成产量大幅下降，据统计[2-3]，倒伏可造成玉米减产15%～20%，导致我国每年玉米产量损失约100 万t。玉米倒伏主要包括根倒、茎倒和茎折断3 种形式。导致玉米倒伏的因素有气候因素、种植密度、土壤质地以及玉米株高、茎秆柔韧性和根系发达程度等自身遗传因素[4-6]。

茎秆力学特征是影响玉米茎秆倒伏的重要因素之一，代表玉米抗倒伏能力的茎秆力学特征有茎秆拉力、抗推力、穿刺强度和压碎强度，且与玉米抗倒伏能力呈显著正相关[7]。茎秆显微结构是影响玉米茎秆抗倒伏性能的另一个重要因素，维管束数目越多、细胞表皮越厚的玉米茎秆强度越好，抗倒性越强。冯素伟等[8]对茎秆显微结构进行分析发现，茎秆显微结构中大维管束数量越多，玉米的抗倒能力越强。前人对玉米倒伏的研究主要集中在茎秆力学特征、形态学特征以及相应的生理生化指标方面，在分子水平对玉米抗倒伏相关基因的挖掘的研究较少。因此，本研究以高抗倒的京农科728、中抗倒的金农738 和低抗倒的先玉335 为材料，在玉米抽雄期取第3 节间茎秆进行转录学测序，挖掘与玉米抗倒性相关的基因，为培育高产和高抗倒伏的玉米品种提供理论依据。

1 材料与方法

1.1 试验设计

试验于2019 年在河北省张家口市沙岭子镇农业科学院（115°05′E，40°6′N）进行，土壤质地为壤质土。供试玉米品种为京农科728（JNK728，高抗倒性，H）、金农738（JN738，中抗倒性，M）和先玉335（XY335，低抗倒性，L），试验采取完全随机区组设计，3 次重复，行距55cm，株距30cm，种植密度67 500 株/hm2，每个品种12 行，行长10m。在玉米抽雄期，每个小区取3 株长势一致的玉米，取第3 节间茎秆液氮速冻后于-80℃保存待测。

1.2 玉米茎秆样本总RNA 的提取、文库构建及Illumina 测序

使用康为全能型植物R N A 提取试剂盒（OminiPlant RNA Kit）提取9 个玉米茎秆样本的总RNA。使用Qubit 2.0 荧光计测定RNA 的浓度，用Agilent 2100 生物分析仪和1%的琼脂糖凝胶电泳检测RNA 的完整性。使用Illumina 公司的TruSeqTM RNA 试剂盒（San Diego）构建测序所需cDNA 文库，用Agilent 2100 生物分析仪检测文库质量。由武汉迈特维尔生物科技公司完成建库测序。

1.3 数据处理、比对及表达量评估

使用Illumina HiSeq 高通量测序平台对cDNA文库进行测序，对生成的原始数据进行处理，除去低质量读段、含接头或引物的读段，利用HISAT2将高质量读段与参考基因组（B73 RefGen_v4）进行序列比对，获取在参考基因组或基因上的位置信息以及测序样品特有的序列特征信息。所有下游分析均基于与参考数据完美匹配或只有一处不匹配的读段。采用FPKM（fragments per kilobase per million mapped reads）计算转录本或基因的表达水平。

1.4 差异表达基因的鉴定与功能注释

为了鉴定与玉米茎折相关基因，用DESeq2 分别对3 个品种组间的差异表达进行分析，基因表达差异倍数在2 倍以上、且FDR＜0.05 确定为差异表达基因（differentially expressed genes，DEGs）。为了进一步鉴定DEGs 的功能，使用Blast2GO（ https://www.blast2go.com/ ） 和 KOBAS 2.1.1(http://kobas.cbi.pku.edu.cn/download.php)进行GO和代谢通路富集分析（FDR ＜0.05）；使用Venny2.1.0 构建Venn 图。

1.5 实时荧光定量PCR（qRT-PCR）验证

为了验证RNA 测序结果的准确性，选择10个DEGs 进行qRT-PCR 验证。采用Primer 5 对DEGs 进行荧光定量引物设计。内参基因为玉米GAPDH，引物信息见表1。使用Bio-Rad iQ5 荧光定量PCR 仪（Bio-RAD，美国），选用2×Fast Super EvaGreen ®qPCR Mastermix（EverbrightInc.，美国）进行荧光定量分析。反应体系为cDNA 模板1µL、前后引物各1µL、ddH2O 7µL、2×AugeGreen TMMaster Mix 10µL。用2-ΔΔCT法[9]计算基因的相对表达量。

表1 荧光定量引物序列Table 1 The primer sequence of RT-qPCR

2 结果与分析

2.1 转录组测序结果统计与质量评估

用HiSeq6000 平台对9 个玉米茎秆样本进行双端测序后共获得59.26Gb 高质量读段，与玉米参考基因组比对后发现，84.84%～87.28%的高质量读段比对到参考基因的唯一位置，3.04%～3.58%的有效数据比对到参考基因的多个位置，Q30 和GC 含量分别高于93.90%和53.16%（表2）。其次，用主成分分析（PCA）评估样本之间相关性，结果（图1）表明，生物学重复之间高度相关，每个材料间明确分离。以上结果证实转录组测序数据是可靠的，可用于进一步分析。

图1 转录组测序样本PCA 分析Fig.1 Principal component analysis of samples used for transcriptome sequencing

表2 转录组测序数据统计Table 2 The summary of RNA-Seq data

2.2 差异表达基因鉴定

为了挖掘玉米抗倒伏相关基因，对3 个抗倒性不同的玉米材料进行品种间基因差异表达分析。由表3 可知，先玉335 与京农科728 比较（L-vs-H）鉴定到差异基因最多，其中4411 个DEGs 上调表达，3368 个DEGs 下调表达。先玉335 与京农738 比较（L-vs-M）鉴定到5373 个DEGs，包括2896 个DEGs 下调表达，2477 个DEGs上调表达。京农738 与京农科728 比较（M-vs-H）鉴定到4905 个DEGs，2831 个DEGs 下调表达，2074 个DEGs 上调表达。3 个比较组共鉴定到10 093 个DEGs，其中2095、769、606 个DEGs分别在L-vs-H、L-vs-M、M-vs-H 分组特异表达，1341 个DEGs 在3 个分组中共同表达（图2）。

图2 分组鉴定到的DEGs Venn 图分析Fig.2 Venn diagram analysis of DEGs observed in comparison groups

表3 差异表达基因统计Table 3 Statistics of differentially expressed genes

2.3 DEGs GO 功能富集分析

为了探究抗倒性不同玉米品种差异表达基因的生物学功能，对鉴定到的DEGs 进行GO 功能富集分析。由图3 可知，生物学过程分类中，L-vs-H分组显著富集到植物型细胞壁的生物发生、光合作用、木质素合成过程，L-vs-M 分组显著富集到植物型次生细胞壁的生物发生，M-vs-H 分组显著富集到次生代谢过程、细胞壁发生、苯丙素的代谢过程。分子功能分类中，L-vs-H、M-vs-H 分组显著富集到血红素结合、单加氧酶活性、碳水化合物结合；L-vs-M 分组显著富集到叶绿素结合、四吡咯结合、糖基转移酶活性。细胞成分分类中，光合体系、光合体系Ⅰ、光合体系Ⅱ在L-vs-H、L-vs-M 分组中显著富集，非原质体在M-vs-H 分组中显著富集。因此，3 个玉米品种抗倒性不同可能是因为富集到相同GO 条目的DEGs 数量不同。

图3 DEGs 的GO 富集分析Fig.3 Gene ontology(GO)enrichment analysis of DEGs

2.4 DEGs 代谢通路富集分析

为了近一步了解DEGs 的功能，对3 个分组鉴定到的DEGs 进行KEGG 代谢通路富集分析。结果（图4）表明，L-vs-H 和M-vs-H 分组均富集到6 条代谢通路，L-vs-M 分组富集到5 条代谢通路。其中，苯丙素生物合成、次生代谢产物的生物合成、类黄酮生物合成在L-vs-H 和M-vs-H 分组显著富集。L-vs-M 分组中，光合作用―天线蛋白、植物―病原体相互作用、次生代谢产物的生物合成显著富集。

图4 DEGs 代谢通路富集分析Fig.4 Metabolic pathway enrichment analysis of DEGs

2.5 3 个玉米品种茎秆的显微结构差异

3 个玉米品种的显微镜可视范围内茎秆维管束数目、单个维管束面积、茎秆表皮细胞厚度如图5所示。先玉335、金农738 和京农科728 显微镜可视范围内茎秆维管束数目分别为32、39、44 个（图5a 和图6a）；先玉335 单个维管束面积最小，京农科728 单个维管束面积最大（图5b 和图6b）；茎秆表皮细胞厚度的顺序依次为京农科728＞金农738＞先玉335（图5c 和图6c）。由此可知，京农科728 的高抗倒性可能与较少的维管束数目、较大的单个维管束面积以及较厚的茎秆表皮细胞厚度有关。

图5 3 个玉米品种茎秆的显微结构Fig.5 Microstructure of the stems of the three maize varieties

图6 3 个玉米品种茎秆的显微结构Fig.6 Microstructures of the stalk of the three maize varieties

2.6 DEGs 的荧光定量验证

为了验证转录组测序的准确性，通过qRT-PCR检测了10 个已报道与茎秆抗倒相关或注释到与茎秆抗倒相关途径的DEGs 表达量，结果（图7）表明，所有10 个基因的表达水平与RNA-seq 测序的结果一致，且相关系数达到0.9171，因此qRT-PCR的分析结果证实了转录组测序结果的可靠性。

图7 10 个DEGs 的RNA-seq 测序结果的qRT-PCR 验证Fig.7 The qRT-PCR validation of the RNA-seq data for the ten DEGs

3 讨论

玉米在我国粮食安全战略中占有重要地位，随着玉米种植密度增加，茎秆增高，茎秆倒伏问题日益严重，成为限制玉米产量和品质的因素之一。作物倒伏指由各种外界因素引发的植株茎秆由自然直立状态到永久错位的现象，是农作物生产中存在的一个普遍性问题。导致作物倒伏的原因主要有遗传、栽培管理和自然环境3 个方面，其中由作物基因型决定的遗传因素是内因，也是决定作物自身是否抗倒的最直接、最根本的因素[10]。植物生长处于风、雨等自然条件的影响下，当作物的茎秆弹性不足以使弯曲的植株恢复时容易发生大面积倒伏。玉米茎秆单位面积内维管束数目和表皮细胞厚度决定了玉米茎秆的抗压强度，本研究通过比较分析不同玉米品种的维管束数目和表皮细胞厚度，结果表明维管束数目越多，表皮细胞越厚，玉米茎秆抗折强度越高。杨硕等[11]对玉米茎秆显微结构研究表明，茎秆维管束数目、大维管束数目和小维管束数目在不同环境下均与茎秆抗倒伏能力呈正相关。同样，对水稻的研究[12]发现，抗倒伏能力越强的水稻，其茎壁厚，大、小维管束数目多。

茎秆抗倒伏是由多个基因控制的复杂的数量性状，目前对于玉米茎秆抗倒伏相关基因与代谢通路的研究较少。因此，本研究对3 个抗倒性不同玉米品种进行转录学测序，并对DEGs 进行了GO 和KEGG 富集分析，报道了玉米抗倒伏的关键差异表达基因和代谢通路。木质素是一类重要的大分子有机物质，由类苯丙酸途径的单体衍生形成，在维持细胞壁结构完整性方面起着关键性作用，茎秆的木质素含量越高作物的抗倒性越强[13]。木质素合成途径中，CoA 连接酶（4CL）引导光合代谢产物向木质素合成途径的流动，处于代谢分支点的位置[14]。研究[10]表明，当4CL 活性比原来提升1 倍时，木质素含量也会相应提高1/4。本研究中苯丙素生物合成、木聚糖生物合成过程、细胞壁生物发生和植物型次生细胞壁的生物发生在3 个分组中显著富集，对大豆茎倒伏的研究[10]表明DEGs 也与这些功能相关。已报道的与茎倒伏相关的苯丙烷类代谢途径中的关键酶苯丙氨酸解氨酶、肉桂酸4-羟化酶和连接植物苯丙烷复合途径和木质素特异生物合成途径的4CL 显著差异表达。推测这些显著富集的通路和差异基因可能与玉米抗倒性相关，其功能还需进一步验证。

4 结论

对3 个抗倒性不同的玉米品种抽雄期茎秆进行转录组学分析，3 个比较组共鉴定到10 093 个DEGs，其中2095、769、606 个DEGs 分别在L-vs-H、L-vs-M 和M-vs-H 分组特异表达，1341 个DEGs在3 个分组中共同表达。GO 富集分析表明，3 个玉米品种抗倒性不同可能是因为富集到相同GO 条目的DEGs 数量不同。代谢通路富集分析表明，L-vs-H、M-vs-H 分组显著富集到苯丙素生物合成、次生代谢产物的生物合成、类黄酮生物合成途径，而L-vs-M 分组显著富集到光合作用―天线蛋白、植物―病原体相互作用途径，注释到这些途径的基因为与玉米抗倒伏相关的关键候选基因。茎秆的显微结构分析表明，京农科728 的抗倒性强可能与茎秆维管束数目多、单个维管束面积大、茎秆表皮细胞厚有关。此外，qRT-PCR 验证了转录组测序的结果。本研究确定的关键基因和代谢途径可作为未来定向克隆和下游分析研究的宝贵遗传资源或选择目标。