玉米穗发芽突变体vp-like8 的转录组分析

郭宇航，郑军，李润植*

（1.山西农业大学农学院，山西太谷030801；2.中国农业科学院作物科学研究所，北京100081）

玉米是重要的粮食和饲料作物，目前已成为我国第一大作物。穗发芽（Preharvest sprouting），又称胎萌（vivipary），是指种子收获前在母体上发芽的现象。穗发芽的种子部分营养和贮藏物质被分解消耗，导致其品质和播种质量下降［1］。我国玉米、水稻、小麦在收获期经常遭遇阴雨连绵的天气，直接导致穗发芽，生产造成重大损失。因此，研究穗发芽相关基因及调控网络对培育耐穗发芽作物品种，保障粮食生产具有重要意义。

目前，在玉米中已经发现了十多个穗发芽突变体。这些突变体根据Robertson［2］提出的标准分为两类：第一类突变体包括vp1/vp4、vp6、vp8、vp10/vp13、vp14和vp15，这些突变体的胚乳和幼苗颜色不受影响；第二类突变体包括vp2、vp5、vp7/ps1、vp9、vp12/lw2、y9、w3和rea1，这些突变体的类胡萝卜素和叶绿素合成受到影响。目前已知，这些穗发芽突变体中有9 个已克隆到基因，包括vp1［3］、vp8［4］、vp10/vp13［5］、vp14［6］、vp15［7］、vp5［8］、vp7/ps1［9］、y9［10］、vp9［11］。VP1 对玉米种子成熟至关重要［12］。vp1突变体对脱落酸（ABA）敏感性降低，不能正常成熟，并引起种子发芽。玉米穗发芽基因VP1和拟南芥的ABI3是同源基因，其编码蛋白是具有多结构域的转录因子，行使激活或者抑制功能取决于启动子的上下游环境［12～15］。VP1 的 3 种蛋白结构域 B1，B2，B3 在不同植物中是高度保守的［12，16～21］。VP1 通过 B3 结构域结合到Sph/RY 元件激活花青素调节基因C1基因的表达［22］。玉米vp1突变体无法激活C1，进而影响糊粉层着色［3］。VP1 还能抑制糊粉层细胞中发芽相关基因如淀粉酶的表达，表明VP1 也能行使转录抑制的功能［14，23，24］。

第二代测序技术（Next-generation sequencing，NGS），又称高通量测序技术（High-throughput sequencing），具有高通量、高准确性、低成本、耗时短的特点，广泛应用于各物种的基因组测序与转录组分析。例如利用转录组测序（RNA-Seq）数据，在全基因组的水平上发掘对玉米籽粒物质积累过程具有重要影响的功能基因［25］。

本研究利用实验室新发现的Vp1的等位突变体vp-like8［26］的杂合 体与自交 系 Mo17 和 Zheng58构建的F2分离群体，利用生物信息学比较分析穗发芽和野生型的转录组，寻找差异表达基因（differentially expressed genes，DEGs），通过基因本体（Gene Ontology，GO）分析和 KEGG 代谢通路富集分析，为揭示玉米vp-like8突变体穗发芽相关基因的表达网络提供理论依据。

1 材料和方法

1.1 试验材料

在玉米收获时，玉米穗发芽的籽粒已干枯致死，因此穗发芽突变体（viviparous）vp-like8只能以杂合体的形式进行保存。利用vp-like8杂合体分别杂交自交系Mo17 和Zheng58，然后再自交，获得具有穗发芽表型分离的F2果穗作为后续转录组测序的样品。

1.2 总RNA 提取及测序

在授粉后30 d，可在果穗上看到明显的穗发芽表型。此时，随机挑选有穗发芽表型分离的F2果穗，发芽籽粒和正常籽粒各取40 粒，剥去种皮，分别混池Mo17 背景的突变体（Mo17_mut）和野生型（Mo17_wt） ， Zheng58 背 景 的 突 变 体（Zheng58_mut）和野生型（Zheng58_wt）等 4 份实验样品，由于利用了不同背景的材料，故每种材料不再单独做生物学重复。液氮研磨样品，利用植物总RNA 提取试剂盒（北京天根生化科技有限公司，Cat#DP432）提取籽粒总RNA，在北京贝瑞和康生物技术有限公司利用HiSeq2000 仪器完成转录组测序。

1.3 差异表达基因的鉴定

测序得到的原始数据，经过去除接头序列和低质量测序片段之后利用软件Tophat2 v2.1.0［27］比对到玉米B73 参考基因组上（B73.AGPv3.27，https：//www.maizegdb.org），所 用参数“-i 5 -I 60000 --no-novel-juncs --microexon-search -r 100 --mate-std-dev 75”，其它参数保持默认，得到mapped data。 利 用 samtools［28］对 bam 文件进行过滤，设置参数“-q 50”。之后利用 HTSeq v0.9.1［29］的htseq-count 脚本计算出唯一比对到参考基因组上的 reads 数，参数为“-f bam -q -a 10 -m union --nonunique none -s no -r pos”。随机挑选基因，使用可视化软件 Integrative Genomics Viewer（IGV）［30］对HTSeq-count 统计计算的read counts 进行验证。最后使用R 语言DESeq2［31］包进行差异表达分析。计算每个基因在突变体与野生型中表达差异倍数（Fold Change，FC），以（|log2FC|＞1）并且差异显著性padj（p 值经过多重检验校正后的值）＜0.01作为判断基因表达差异是否显著的阈值。

1.4 qRT-PCR 验证

随机挑选7 个差异表达基因（G R M Z M 2 G 133398；GRMZM2G382534；GRMZM2G082487；GRMZM2G123107；GRMZM2G126261；GRMZM 2G017290；GRMZM5G835903）进行验证。将 1.2中保留下的RNA 样品，利用反转录试剂盒（北京天根生化科技有限公司，Code# AH311-02）反转录为 cDNA。利用 Primer3［32］设计特异性引物，以GAPDH 为内参（表 1），利用 ABI 7300 实时定量PCR 仪进行测定。qRT-PCR 体系为Dye Ⅰ 0.4 μL ，Q-mix 10 μL，ddH2O 7.8 μL，F/R 引物各 0.4 μL，cDNA 1 μL，qRT-PCR 反应条件为 95 ℃预变性 2 min，95 ℃变性 10 s，58 ℃退火 20 s，72 ℃延伸30 s，设置 40 个循环，在 72 ℃延伸 30 s 阶段采集荧光，并添加溶解曲线，结果用2-ΔΔCt计算穗发芽籽粒与正常籽粒中基因表达差异。试验包含2 次生物学重复，每个生物学重复进行3 次技术重复。

1.5 差异表达基因的注释

基因本体（Gene Ontology，GO）分析是根据GO 功能显著性富集来分析差异表达基因的功能，本研究利用AgriGo v2.0 在线分析网站（http：//systemsbiology.cau.edu.cn/agriGOv2/index.php）进行分析，每条条目至少有3 个基因，并且FDR＜0.05 作为筛选标准。

京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）数据库可以系统分析基因的代谢通路，本研究利用KOBAS3.0［33，34］在线分析网站（http：//kobas.cbi.pku.edu.cn/index.php）进行分析，P 值经 Benjamini and Hochberg 方法校正后，筛选Corrected P-Value＜0.01的通路作为显著富集的通路。

表1 qRT-PCR 所用引物Table 1 Primers used in qRT-PCR

2 结果与分析

2.1 序列比对结果

在本研究中，Mo17 背景的突变体（Mo17_mut）和野生型（Mo17_wt），Zheng58 背景的突变体（Zheng58_mut）和野生型（Zheng58_wt）等4 份样品的测序原始数据，经过去除接头序列和低质量测序片段，之后通过Tophat2 比对到玉米B73 参考基因组（B73.AGPv3.27），比对率均达80%以上（表2），可以用于进一步分析。对于每份数据的比对结果，通过HTSeq-count 对其片段计数（read counts），以基因组注释为基础，计算每个基因区域中比对上的片段数，以此作为每个基因的原始表达数据。以GRMZM2G000743 基因为例，利用IGV 对HTSeq-count 统计计算的read counts 进行验证，在 Mo17_mut 材料中的 read counts 为 12，IGV 结果显示有 12 个成对的 reads 比对到该基因区域（图1），与计算得到的结果一致，验证了HTSeq-count 统计计算read counts 结果的准确性。

表2 Tophat 比对结果Table 2 The tophat mapped results

2.2 差异表达基因分析

分别以Mo17 背景的野生型和Zheng58 背景的野生型数据为对照，计算每个基因在突变体与野生型中表达量的比值，即表达差异倍数。根据DESeq2 分 析结果，以|log2FC|＞1 并且 padj＜0.01 为条件进行筛选（图2A），共筛选出27 50 个DEGs。其中，1 812 个基因表达上调，占总DEGs 的65.89%左右，938 个基因表达下调，占总DEGs 的34.11%左右（图2B）。Vp1基因（GeneID：GRMZM2G1333 98）确实存在于差异基因中，并且在突变体中表达下调（表3）。

2.3 qRT-PCR 验证

本研究从发现的DEGs 中随机选取7 个基因进行荧光定量PCR 验证，通过与差异表达分析结果比较，虽然差异基因在上调和下调倍数上存在差异，但是总体变化趋势是一致的，验证了RNA-seq 分析结果的可靠性（图3）。

2.4 差异表达基因GO 富集分析

图1 利用IGV 验证read counts 结果Fig.1 Using IGV to validate results of read counts

图2 差异表达基因火山图和数量Fig.2 Vocalno plot and number of DEGs

通过GO 功能分类将差异表达基因进行功能分类注释，1 791 个差异基因得到注释，分布于699 个GO 分类条目。在生物学过程（biological process）、细胞组分（cellular component）和分子功能（molecular function）中所占的比例分别为 49.6%（347）、11.9%（83）、38.5%（269）（图4）。进一步对差异基因的生物进程类别富集分析发现（图5），碳水化合物代谢过程（carbohydrates metabolic process）、细胞氮化合物代谢过程（cellular nitrogen compound metabolic process）、脂质代谢过程（lipid metabolic process）、光合作用（photosynthesis）和对刺激响应（response to stress）显著富集。穗发芽需要大量能量，涉及多种代谢过程，我们发现在碳水化合物代谢过程中，α-淀粉水解酶基因和其他糖基水解酶基因在突变体中表达上调；在细胞氮化合物代谢过程中，多种与氮代谢相关的基因表达上调；在脂质代谢过程中，酰基辅酶A 氧化酶基因表达上调，并且与ABA 合成相关的β-胡萝卜素羟化酶基因表达下调；在光合作用中，多种光合作用相关的基因表达上调；在胁迫响应过程中，多个过氧化物酶基因表达上调（表3）。这些显著富集的GO 类别的相关基因共同作用，为玉米籽粒发芽提供能量和保护，与玉米穗发芽是紧密相关的。

图3 qRT-PCR 验证随机挑选的DEGsFig.3 Validation of the random selected DEGs by qRT-PCR

表3 部分差异表达基因信息Table 3 Part of differentially expressed genes’informatation

续表

2.5 差异基因代谢途径分析

图4 差异表达基因的GO 分类Fig.4 Go annotation of differentially expressed genes

图5 差异表达基因的生物进程（有向无环图）Fig.5 Biological process of differentially expressed genes（Directed acyclic graph）

利用KEGG 数据库对差异表达基因所参与的代谢途径进行显著性富集分析，结果显示，这些差异基因分布于110 个代谢通路中，分别为代谢通路（Metabolic pathways）、次生代谢产物的生物合成（Biosynthesis of secondary metabolites）、苯丙素生物合成（Phenylpropanoid biosynthesis）、光合作用（Photosynthesis）、植物激素信号转导（Plant hormone signal transduction）、碳代谢（Carbon metabolism）、淀粉和蔗糖代谢（Starch and sucrose metabolism）等。其中代谢通路和次级代谢产物生物合成通路富集到的基因是最多的，分别有314 和211 个基因被富集（表4）。植物内源激素与植物生长发育密切相关，有48 个基因富集到植物激素信号传导通路（图6），图中红框为差异表达基因在通路中位置，涉及生长素（Auxin）、细胞分裂素（Cytokinin，CTK）、赤霉素（Gibberellin，GA）、脱落酸（Abscisic acid，ABA）、乙烯（Ethylene，ETH）等，其中ABA 直接影响种子休眠，我们发现蛋白磷酸酶2C基因（PP2C）、丝氨酸/苏氨酸蛋白激酶基因（SnRK2）、ABA 响应元件结合因子（ABF）在该通路中受到影响。这些受影响的代谢通路中的基因为研究玉米穗发芽相关表达网络提供了丰富的基因资源。

表4 显著富集的KEGG 代谢通路Table 4 Highly enriched KEGG terms

3 讨论与结论

玉米穗发芽是一个复杂的生物学过程，影响因素较多，往往各因素相互影响。在本研究中，借助本实验室新发现的玉米穗发芽vp-like8突变体［26］进行研究。赵燕［35］对玉米穗发芽差异基因分析发现光合作用通路最显著富集，我们发现与光合作用相关的基因大部分都表达上调（表3），说明光合作用与穗发芽确实存在紧密联系。花青素调节基因C1与糊粉层着色相关，本研究中，该基因表达下调（表 3），与前人研究结果一致［3］。VP1 能够抑制糊粉层中和发芽相关的α-淀粉酶的表达［14，23］，本研究中我们发现在淀粉和蔗糖代谢途径中，α-淀粉酶基因表达上调（表3）。这些结果印证了我们分析的可靠性，同时为利用转录组测序技术研究玉米穗发芽机制提供了支持。

研究表明，种子内源激素 ABA［36，37］对种子发芽有重要作用。Suzuki 等人研究发现［38］，VP1 能影响ABA 信号传导中的bZIP基因和蛋白磷酸酶2C基因（PP2C）的表达。本研究中，同样发现这些基因在ABA 信号传导通路中是差异表达的，并且还发现丝氨酸/苏氨酸蛋白激酶基因（SnRK2）也在ABA 信号传导中受到影响（表3）。同时，还发现生长素转运载体，生长素响应蛋白相关基因的表达也有显著差异（表3）。综上，这些差异基因表达水平的变化为穗发芽相关调控网络的研究提供了参考。

图6 植物内源激素信号传导Fig.6 Plant hormone signal transduction

其它作物穗发芽相关研究也较为丰富。廖泳祥等［39］通过比较不同的水稻材料种子，发现在不耐穗发芽的水稻材料的种子中α-淀粉酶含量明显高于耐穗发芽的水稻种子，类似结果在小麦中也见报道［40］。孙果忠等［41］人对不同生理状态的小麦离体胚萌发过程研究发现，ABA 可以抑制未成熟胚与休眠成熟胚中α-淀粉酶-1 的表达，但在不休眠成熟胚中则诱导α-淀粉酶-1 表达。此外，也有部分学者认为，种皮颜色、子粒和穗部其它物理因素与穗发芽存在关联。

本研究利用转录组测序技术系统分析了玉米vp-like8突变体和野生型的籽粒转录组。通过比较突变体和野生型间的表达水平，发现了2 750 个差异表达基因，其中1 812 个基因在突变体中上调表达，938 个下调表达。通过功能注释分析发现，这些差异基因主要涉及代谢通路、光合作用、植物激素信号传导等途径，表明这些途径在穗发芽过程中具有一定作用。本文通过转录组测序的方法初步研究了vp-like8基因的表达差异，为深入了解Vp1基因功能和玉米穗发芽的机制及基因调控网络提供了数据信息。