高温胁迫下不同玉米材料间的基因差异表达分析

曹丹，张红梅，杨海鹏，赵欢，刘小红

(1.西华师范大学生命科学学院，四川 南充 637000；2.山西省农业科学院玉米研究所，山西 忻州 034000)

玉米(Zeamays)是世界上最重要的3大谷类作物之一，与水稻和小麦并列，在世界农业中占有重要地位[1].高温是限制全球玉米产量的最关键的非生物胁迫因素之一，世界上许多国家的玉米生产都遭受高温危害[2-5].在中国，黄淮海地区是我国夏玉米集中产区，几乎每年都会遭遇高温危害.近年来，伴随着全球气候变暖，干旱和高温协同胁迫已经成为该地区玉米生产中频繁遇到的一种自然灾害，且这一危害几乎发生在玉米整个生育期，严重影响着玉米产量和品质[6].

高温胁迫可引起植物的生理、分子和生化变化，干扰细胞和整个植物的生长发育过程，从而对作物的生产产生不利影响[7].对玉米来说，除了在其他生长阶段可能遭遇高温胁迫外，在生殖生长及受精结实阶段特别容易遇到高温天气，导致玉米减产[8].高温对花粉(雄配子)发育的影响要大于花丝(雌配子)[9].在雌雄穗分化至籽粒成熟阶段如果遭遇持续的高温胁迫，可能造成败育小花数量增加、花粉结构破坏、雌穗结实率降低，最终导致减产.研究发现，不同玉米材料耐高温胁迫能力存在差异[10-11]，目前在国内耐高温的玉米材料较少，仅有郑单958和中地88等较少的耐高温玉米品种被推广应用.因此，以不同耐高温胁迫能力的玉米花粉作为材料，研究在高温环境下基因的差异表达具有重要意义.

基于此，本研究以耐高温胁迫的中地88和高温敏感的先玉335两个玉米品种作为试验材料，在高温胁迫环境下取其花粉作转录组测序分析，了解基因的差异表达情况，并进一步对差异表达基因作功能注释和代谢通路分析，以期为进一步阐明玉米耐高温胁迫的分子遗传机制奠定理论基础，为克隆应用玉米耐高温胁迫基因的分子育种提供参考.

1 材料与方法

1.1 试验材料

本试验选用了两个玉米材料，一个是耐高温胁迫的中地88，另一个是对高温敏感的先玉335，两个玉米材料均来自山西省农业科学院玉米研究所.

1.2 试验方法

2019年将中地88和先玉335两个玉米品种种植于山西省农业科学院玉米研究所实验大棚，每个材料种植3行，行长3 m，每行10株.从播种到开始散粉前第5天，未作增温处理，玉米生长环境包括温度、湿度和光照等条件与田间自然生长的玉米材料一样.在开始散粉前的第4天作增温处理，但湿度和光照等条件与室外保持一致，在盛花期头一天下午对雄穗进行套袋，第2天12时30分时监测到雄穗所处位置的温度达到42 ℃，通过大棚通风手段维持该温度1 h，在13时30分，取两个材料的花粉，各取4个重复，每个重复取约400 mg花粉装进2 mL离心管，立即液氮速冻，该过程在实验大棚内完成.然后在实验室提取花粉总RNA，对符合要求的总RNA每个材料各选3份(作为3个生物学重复)，中地88作为处理，用T1、T2和T3表示；先玉335作为对照，用CK1、CK2和CK3表示；再进一步完成mRNA纯化、反转录、建库及测序等工作(该工作由上海派森诺生物科技股份有限公司完成).

1.3 数据统计分析

1.3.1 参考基因组选择 试验结果分析选取玉米B73作为参考基因组，该参考物种在数据库中注释的基因数目及注释比例为：NCBI_GI：16 607(41.94%)；KEGG：6 715(16.96%)；GO：30 113(76.06%)；EC：2 921(7.37%)；UniProtID：16 263(41.07%)；Ensembl：39 591(100%).

1.3.2 序列比对分析 对试验返回的结果，使用TopHaT2的升级版HISAT2(http://ccb.jhu.edu/software/hisaT2/index.shtml) 软件分析，将过滤后的Reads比对到参考基因组上，并对获得的基因在基因组上的区域分布进行统计.

1.3.3 基因表达量分析 使用HTSeq统计比对到每一个基因上的Read Count值作为基因的原始表达量.采用FPKM(Fragments Per Kilobase Million)对表达量进行标准化，再绘制表达基因的FPKM密度分布的小提琴图.

1.3.4 基因差异表达分析 根据基因在两种玉米材料中的表达量，以先玉335材料为对照，以中地88材料为处理，分析上调基因、下调基因及无显著差异表达基因的数量，并绘制差异表达基因的火山图和在不同染色上的基因组圈图.

对差异表达基因的生物学功能进行GO(gene ontology)富集及功能注释，对差异表达基因的GO富集分析结果，按照分子功能(molecular function，MF)、生物过程(biological process，BP) 和细胞组分(cellular component，CC) 进行GO 分类，每个GO 分类中挑选P-value最小，即富集最显著的前10个GO term条目进行展示.此外，还对差异表达基因作KEGG(kyoto encyclopedia of genes and genomes)富集.根据差异表达基因的KEGG富集分析结果，挑选P-value最小即富集最显著的前20个Pathway进行展示.

2 结果与分析

2.1 Reads与参考基因组的比对

2.1.1 基本信息比对 使用TopHaT2的升级版HISAT2软件将过滤后的 Reads 比对到参考基因组上.序列比对的基本信息统计结果见表1.对T1、T2、T3、CK1、CK2和CK3面言，比对上参考基因组的序列总数占比对的序列总数的百分比分别为92.50%、92.83%、92.24%、90.62%、91.25%和92.54%，平均值为91.96%；比对到多个位置的序列总数占比对上参考基因组的序列总数的百分比分别为5.88%、5.87%、6.29%、6.69%、6.47%和6.69%，平均值为6.34%；只比对到一个位置的序列总数占比对上参考基因组的序列总数分别为94.12%、94.13%、93.71%、93.31%、93.53%和93.31%，平均值为93.66%.

表1 Reads与参考基因组间的基本信息比对Table 1 Basic information comparison between reads and reference genome

2.1.2 比对区域分布 将比对到基因组上的Reads分布情况进行统计，结果见表2.对T1、T2、T3、CK1、CK2和CK3面言，比对到基因区域的Reads总数占比对事件发生的总次数的百分比分别为94.24%、94.07%、93.53%、93.76%、94.10%和94.07%，平均值为93.96%；比对到基因间区的Reads总数占比对事件发生的总次数的百分比分别为5.76%、5.93%、6.47%、6.24%、5.90%和5.93%，平均值为6.04%；比对到外显子区域的Reads总数占比对到基因区域的 Reads 总数的百分比分别为99.57%、99.59%、99.56%、99.45%、99.50%和99.45%，平均值为99.52%.

表2 比对到参考基因组上的Reads分布Table 2 Distribution of the Reads mapped to reference genome

2.2 基因表达的结果与分析

2.2.1 基因表达量 采用FPKM对表达量进行标准化后，将表达量分为不同的区间，对各样本在不同表达量区间内基因的数目进行统计，结果见表3.表达量值在0～0.01区间的基因数量最多，所有6个样本都超过了20 000个，平均为21 228个，占总量的60.96%；随表达量范围值的提高，对应基因数量逐渐减少，在>1000表达量区间所有6个样本的基因数量都低于160个，平均为155.5个，仅占总量的0.45%.由此表得出，低表达的基因数量较多，而高表达的基因数量较少.

表3 不同表达量区间的基因数量统计Table 3 Statistics of the number for the genes in different expression quantity intervals

2.2.2 表达基因的FPKM密度分布 基因在各个样品中的表达量利用小提琴图进行统计分析，结果如图1所示.从图1可以看出，这6个样本具有相似的基因表达模式，即中等表达的基因占绝大多数，而低表达和高表达的基因仅占一小部分.

2.3 基因差异表达结果与分析

2.3.1 差异表达基因数量 中地88(Case:Traitment)相对于先玉335(Control：CK)材料而言，差异表达基因共计有1 822个，相比于Control，上调表达基因有857个，下调表达基因有965个.另还有22 731个基因在这两个玉米材料中表达差异不显著.绘制的差异表达基因的火山图见图2.从图2可看出，左侧(蓝色)为Case相比于Control下调基因，右侧(红色)为Case相比于Control上调基因，左右差异基因分布大致对称.

2.3.2 差异表达基因的染色体分布 对在这两个玉米材料中表达的基因在染色体上的分布作基因组圈图，结果如图3所示.从上调基因、下调基因及无显著差异表达基因在整个玉米10条染色体上的分布来看，分布较均匀.

x为FPKM值；图中盒型中间的横线是中位数，盒型的上下边缘为75%，上下限为90%，外部形状为核密度估计.x means FPKM value；The horizontal line in the middle of the box is the median，the upper and lower edges of the box are 75%，the upper and lower limits are 90%，and the outer shape is the kernel density estimation.图1 FPKM密度分布Figure 1 Density distribution of FPKM

2.4 差异表达基因功能富集分析

2.4.1 GO富集分析 对差异表达基因进行GO富集，结果显示，与细胞组分(CC)相关的GO条目有491个，占总量的12.68%；与分子功能(MF)相关的GO条目有1 062个，占总量的27.42%；与生物过程(BP)相关的GO条目有2 320个，占总量的59.90%.对GO富集分析结果，每种分类挑选P-value最小，即富集最显著的前10个GO条目进行展示，结果见图4.CC分类中，P-value值最小的基因的ID号为GO:0005887，与质膜组成成分相关；MF分类中，P-value值最小的基因的ID号为GO:0004046，与氨酰基转移酶活性相关；BP分类中，P-value值最小的基因的ID号为GO:0044281，与小分子代谢过程相关.

x表示差异倍数；图中两条竖虚线为2倍表达差异阈值；横虚线为P值的阈值(0.05).红点表示该组上调基因，蓝点表示该组下调基因，灰点表示非显著差异表达基因.x means fold change；The two vertical dashed lines in the figure are 2 times of the expression difference threshold；the horizontal dashed line is the threshold value of P-value(0.05).The red，blue and gray dots represent the up-regulated，down-regulated and non-significantly differentially expressed genes，respectively.图2 差异表达基因的火山图Figure 2 Volcano map of differentially expressed genes

2.4.2 KEGG富集分析 参考KEGG通路系统的基因组信息数据库，对检测到的在这两个玉米花粉材料中呈现差异表达的基因的代谢通路作富集分析.富集到的差异基因数目为349个，共涉及到102种通路，这些小通路可进一步归为6个大通路，包括细胞过程3个、环境信息处理4个、遗传信息处理20个、人类疾病1个、新陈代谢72个和有机系统2个.上调和下调表达基因分别有143和206个(注：一个基因可能同时富集于不同通路).挑选出的P-value最小即富集最显著的前30个Pathway展示结果见图5，其中1个与环境信息处理有关，4个与遗传信息处理有关，25个与新陈代谢有关.

最外圈是染色体条带.红色和绿色分别为上调和下调基因的log2x(Fold Change) 值的柱状图，灰色为无差异表达基因的log2x(Fold Change)值的散点图.The outermost circle is the chromosomal banding.The red and green are the histogram of the log2x(Fold Change) value of up-regulated and down-regulated genes respectively，and the gray is the scatter graph of the log2x(Fold Change) value of the undifferentially expressed genes.图3 基因组圈图Figure 3 Genome circos

3 讨论与结论

温度是影响植物生长非常重要的因素[12-13]，高温胁迫是最常见的非生物胁迫之一，它能显著抑制植物生长发育[14].高温胁迫一般会降低水分含量，影响植物光合作用和呼吸作用[15].研究发现，高温胁迫已导致全球植物的产量下降和品质降低[16-17].

1：质膜组成成分;2：DNA指导的RNA聚合酶II，全酶；3：胞质；4：核小体；5：胞质大核糖体亚基；6：DNA包装复合物；7：核DNA指导RNA聚合酶复合物；8：大核糖体亚基；9：MLL1/2复合物；10：MLL1复合物；11：氨基酰化酶活性；12：蛋白异源二聚化活性；13：基本转录机制结合；14：基本RNA聚合酶II转录机制结合；15：肽基转移酶活性；16：谷胱甘肽水解酶活性；17：辅因子结合；18：核苷跨膜转运体活性；19：磷酸吡哆醛结合；20：维生素B6结合；21：小分子代谢过程；22：赤霉素生物合成过程；23：应答低氧含量；24：应答氧水平；25：糖脂分解代谢过程；26：鞘糖脂分解代谢过程；27：萜类生物合成过程；28：磷代谢过程的正向调节过程；29：磷酸盐代谢过程的正向调节过程；30：有机酸代谢过程.1：Integral component of plasma membrane；2：DNA-directed RNA polymerase II，homoenzyme；3：Cytosol；4：Nucleosome；5：Cytosolic large ribosomal subunit；6：DNA packaging complex；7：Nuclear DNA-directed RNA polymerase complex；8：Large ribosomal subunit；9：MLL1/2 complex；10：MLL1 complex；11：Aminoacylase activity；12：Protein heterodimerization activity；13：Basal transcription machinery binding；14：Basal RNA polymerase II transcription machinery binding；15：Peptidyltrasferase activity；16：Glutathione hydrolase activity；17：Cofactor binding；18：Nucleoside transmembrane transporter activity；19：Pyridoxal phosphate binding；20：Vitamin B6 binding；21：Small molecule metabolic process；22：Gibberellin biosynthetic process；23：Response to decreased oxygen levels；24：Response to oxygen levels；25：Glycolipid catabolic process；26：Glycosphingolipid catabolic process；27：Terpenoid biosynthetic process；28：Positive regulation of phosphorus metabolic process；29：Positive regulation of phosphate metabolic process；30：Organic acid metabolic process.图4 差异表达基因的GO注释和分类Figure 4 GO annotation and classification of differentially expressed genes

玉米是世界上最重要的农作物之一，高温胁迫同样会显著降低其产量和品质[18].因此，研究玉米对高温胁迫响应的分子机制及选育耐高温胁迫的玉米品种具有重要意义.

虽然已有关于玉米耐高温胁迫的基因克隆或表达方面的研究，但是以前关于基因克隆的研究主要限于单个或少数几个基因的克隆表达分析，如Wang等[19]从玉米中克隆了转录因子ZmWRKY106基因，将其转入拟兰芥，过量表达的结果显示该基因会明显提高对干旱和高温胁迫的抗性；Ma等[20]研究发现ZmbZIP4基因的过表达可以调节ABA的生物合成及根的发育，从而提高玉米包括高温胁迫等在内的多种抗性；Li等[21]也发现玉米ZmbZIP60基因的表达会应答于高温胁迫.由于玉米对高温胁迫的抗逆性属多基因控制的数量性状，因而对单个基因的克隆表达研究还不能完全满足育种需要，有必要在基于转录组测序的在全基因组水平下研究耐高温胁迫基因.

基于转录组测序对玉米耐高温胁迫的基因差异表达以前也有类似研究，且发现了大量与高温胁迫相关的差异表达基因，但以前的研究取材部位主要限于玉米叶片[1，18，22-23].然而在包括中国在内的许多国家，在玉米散粉期常遭遇高温天气，高温胁迫导致花粉活力下降，从而影响玉米的授粉、结实，使其大幅减产，因此，选择以花粉为材料，基于转录组高通量测序，研究不同玉米材料间的基因差异表达对于认识玉米耐高温胁迫的分子遗传机制具有重要意义.

本试验在研究方法、材料选择或取材部位等方面，不同于以前的研究[19-23]，本试验是以中地88(耐高温胁迫)和先玉335(高温敏感)两个玉米品种的花粉为研究对象，研究其基因表达情况，共发现差异表达基因1 822个，相对于对照(先玉335)而言，中地88上调表达基因有857个，下调表达基因965个.另外还发现有22 731个基因在两个材料中表达差异不显著.对差异表达基因进行GO功能注释，结果显示大部分基因都与细胞组分、分子功能或生物过程功能相关.KEGG功能富集结果显示，这些差异表达基因可以归为6个Level 1分类(包括102个Levle 2通路).综合分析这些差异表达基因，发现有7个基因直接应答于热激效应，其中2个上调表达，5个下调表达；有3个基因与热激蛋白结合有关，其中上调表达基因1个，下调表达基因2个.在这12个与玉米耐高温胁迫相关的基因中，多数为在玉米中新发现的与耐高温胁迫相关的基因.依据基因的ID号和玉米参考基因组序列信息，参考以前的研究[24-26]，进一步对其基因全长序列、cDNA序列及表达蛋白的生物学特性等作进一步研究，关于这些基因的分子结构及遗传功能研究正在进行之中.本实验的研究结果为进一步阐明玉米耐高温胁迫的分子遗传机制奠定了一些理论基础，为在植物中克隆应用这些耐高温胁迫基因提供了参考.

1：植物MAPK信号通路；2：核糖体；3：碱基切除修复；4：硫传递系统；5：RNA运输；6：其他多糖降解；7：糖酵解/糖异生；8：α-亚麻酸代谢；9：淀粉和蔗糖代谢；10：苯并噁唑嗪酮类化合物生物合成；11：鞘脂代谢；12：丙酮酸代谢；13：丁酸代谢；14：脂肪酸降解；15：磷酸戊糖途径；16：氰氨基酸代谢；17：抗坏血酸和醛糖酸代谢；18.苯丙氨酸代谢；19：咖啡因代谢；20：缬氨酸、亮氨酸和异亮氨酸降解；21：组氨酸代谢；22：苯丙氨酸、酪氨酸、色氨酸生物合成；23：果糖和甘露糖代谢；24：亚油酸代谢；25：牛磺酸与次牛磺酸代谢；26：烟酸和烟酰胺代谢；27：半乳糖代谢；28：硫胺代谢；29：有机含硒化合物代谢；30：泛醌和其他萜醌生物合成1：MAPK signaling pathway - plant；2：Ribosome；3：Base excision repair；4：Sulfur relay system；5：RNA transport；6：Other glycan degradation；7：Glycolysis / Gluconeogenesis；8：alpha-Linolenic acid metabolism；9：Starch and sucrose metabolism；10：Benzoxazinoid biosynthesis；11：Sphingolipid metabolism；12：Pyruvate metabolism；13：Butanoate metabolism；14：Fatty acid degradation；15：Pentose phosphate pathway；16：Cyanoamino acid metabolism；17：Ascorbate and aldarate metabolism；18：Phenylalanine metabolism；19：Caffeine metabolism；20：Valine，leucine and isoleucine degradation；21：Histidine metabolism；22：Phenylalanine，tyrosine and tryptophan biosynthesis；23：Fructose and mannose metabolism；24：Linoleic acid metabolism；25：Taurine and hypotaurine metabolism；26：Nicotinate and nicotinamide metabolism；27：Galactose metabolism；28：Thiamine metabolism；29：Selenocompound metabolism；30：Ubiquinone and other terpenoid-quinone biosynthesis图5 KEGG通路富集结果柱状图Figure 5 Histogram of enrichment results for KEGG pathway