玉米自交系籽粒响应高温胁迫的转录组分析

汤 彬,耿存娟,曾 强,郭欢乐,李 涵,曹钟洋,邓力超,彭 明,周 虹,陈志辉

(1.湖南省作物研究所,湖南 长沙 410125;2.湖南智行三湘农林科技有限公司,湖南 长沙 410205)

在全球变暖的大背景下,气候系统的不稳定性加剧,全球气候变化引起的高温热害已经逐渐成为限制玉米生产的主要非生物胁迫因子[1]。解析玉米耐高温的分子机理,提高玉米品种的耐高温能力,兼具重大的理论意义和应用价值,受到国内外学者的广泛关注[2-4]。玉米生长发育适宜的日平均温度为28～32 ℃,当温度超过32 ℃会显著影响玉米正常的生长发育过程,产生光合作用受阻、蒸腾速率增加、花粉失活、籽粒灌浆期缩短等不利因素,最终造成减产甚至绝收[1]。

玉米是湖南省第一大旱粮作物,常年种植面积保持在40万hm2左右。湖南省属亚热带季风湿润气候,每年7—8月为此地区常年高温时期,正好也是春玉米籽粒灌浆和产量形成的关键时期。在气候变化影响下,高温显著降低了湖南春玉米的光温生产潜力[5]。因此,提高玉米籽粒形成期的耐高温能力是选育适应本地区气候特点的耐高温玉米新品种的重要育种目标。

转录组测序是研究植物转录水平调控网络的分子生物学技术,可以快速、有效地鉴定植物响应高温胁迫相关的基因及其代谢通路[6-7]。尽管转录组学测序技术(RNA-Seq)已经鉴定了大量玉米响应高温胁迫相关的差异表达基因,但是已有研究的取样部位主要是玉米叶片[8]和花药[9],对籽粒响应高温胁迫的研究相对较少,且多数研究都只涉及单一玉米自交系[10]或杂交种[11],缺乏对不同耐高温自交系基因型差异的比较分析,无法揭示不同耐高温优良等位基因的功能差异。

表型鉴定和多组学方法(转录组、蛋白质组和代谢组)联合分析表明,糯玉米籽粒发育早期(授粉后1～15 d)的高温胁迫会影响整个籽粒发育时期,籽粒中生长素、脱落酸和水杨酸信号通路相关基因和蛋白质的差异表达可能参与调控高温胁迫下籽粒代谢过程[12]。鉴于此,本研究在前期耐高温自交系筛选的基础上,通过对耐高温自交系XN202和敏感自交系CT110授粉后1～15 d的籽粒进行高温胁迫处理,研究不同耐高温玉米自交系籽粒响应高温胁迫的基因表达差异,在基因水平解析不同玉米种质响应高温胁迫的分子机制,并将DEGs与已报道的玉米耐高温相关性状QTL/SNP定位结果进行联合分析,以期进一步挖掘玉米籽粒响应高温胁迫的关键基因,为后续耐高温玉米新品种培育提供技术支撑。

1 材料和方法

1.1 试验材料

供试材料耐高温自交系XN202和敏感自交系CT110分别来源于湖南省农业科学院作物研究所和国家玉米产业技术体系。

1.2 高温胁迫处理

盆栽试验于2021年4—7月在湖南省作物研究所进行。收集田间土壤置于盆钵(高38 cm、直径43 cm),每盆施12 g复合肥(N-P2O5-K2O为15%-15%-15%)作为基肥,并在每个盆钵中种植3粒种子,出苗后定苗至1株,拔节期每盆追施6 g尿素(N 46%),定期浇水使土壤含水量保持在70%左右。所有植株在室外自然条件下生长至开花期,吐丝后3 d进行人工授粉,授粉后1 d选择生长一致的植株移入温室大棚进行高温处理。以室外自然温度正常生长的植株为对照,每个处理9盆。高温处理后15 d,选3株玉米果穗中部籽粒混合作为1个生物学重复,包含3个生物学重复,于液氮中速冻后将样品保存在-80 ℃的超低温冰箱。高温处理下,耐高温的自交系分别命名为THT1、THT2、THT3;不耐高温的自交系分别命名为SHT1、SHT2、SHT3。室外自然温度正常生长条件下,耐高温的自交系分别命名为TCK1、TCK2、TCK3;不耐高温的自交系分别命名为SCK1、SCK2、SCK3。由北京百迈客生物科技有限公司开展籽粒RNA样品的提取和质量检测,构建测序文库并进行质量控制,最后完成所有样品的高通量测序。

1.3 数据处理与分析

使用百迈客云平台BMKCloud(www.biocloud.net) 提供的生物信息学分析流程进行分析。将测序数据与参考玉米基因组数据库(B73_RefGen_v3版本)进行比对,采用FPKM(Fragments per kilobase of transcript per million fragments mapped)值进行归一化后作为衡量基因表达水平的指标。差异表达基因通过DESeq2_edgeR软件进行筛选与检测,筛选标准为2组(正常对照和高温胁迫处理)样品间基因表达量的差异倍数(Fold change)≥2 同时错误发现率(False discovery rate,FDR)<0.01。最后,提取各差异表达基因集的注释信息,使用GOseq R和 KOBAS软件包进行差异表达基因集的GO节点或KEGG通路富集分析(FDR<0.05)。柱形图的绘制使用GraphPad Prism v.8.0.0软件(San Diego,CA,USA)完成。

2 结果与分析

2.1 高温胁迫对玉米自交系的影响

采用分期播种的方法对511份玉米自交系进行开花期耐高温综合评价。以单穗总质量、单穗籽粒质量和单穗粒数3个性状耐高温系数(耐高温系数=高温性状值/对照测定值)的平均值,作为耐高温综合系数,将玉米自交系的耐高温特性分为5个等级(极强:耐高温系数≥0.8;较强:0.8>耐高温系数≥0.6;中等:0.6>耐高温系数≥0.4;较弱:0.4>耐高温系数≥0.2;极弱:耐高温系数<0.2),筛选出耐高温自交系XN202和敏感自交系CT110(图1)。

图1 玉米自交系耐高温田间鉴定

2.2 转录组测序质控分析

分别提取XN202和CT110授粉后15 d的籽粒RNA用Illumina平台进行转录组测序,将正常条件和高温胁迫下耐高温自交系的样品分别命名为TCK-1-2-3和THT-1-2-3,不耐高温自交系的样品分别命名为SCK-1-2-3和SHT-1-2-3。12个样品(2个品种、2个处理、3次生物学重复)共获得83 061 999 806过滤碱基数(Clean bases),各样品过滤碱基数(Clean bases)均达到5 917 021 526,Q30碱基百分比在92.59%及以上,GC含量约为53%,经过质量剪切统计后获得各样品的过滤序列(Clean reads)与指定的参考基因组进行序列比对,唯一配对序列从74.40%到81.59%不等(表1),表明本研究籽粒转录测序的质量较好,参考基因组选择合理,可以进行后续分析。

表1 测序数据统计

2.3 不同耐高温玉米自交系籽粒基因表达差异

为研究正常对照和高温胁迫下不同耐高温玉米自交系籽粒基因表达变化,使用DESeq2_edgeR软件鉴定差异表达基因(Differentially expressed gene,DEG)。SCK-1-2-3 vs SHT-1-2-3对比组共筛选出1 012个表达变化水平大于2倍的DEGs,其中上调DEGs共有386个,下调DEGs共有626个;TCK-1-2-3 vs THT-1-2-3 对比组共筛选出1 517个表达变化水平大于2倍的DEGs,其中上调DEGs共有801个,下调DEGs共有716个;TCK-1-2-3 vs THT-1-2-3对比组与SCK-1-2-3 vs SHT-1-2-3对比组共检测到142个共同的DEGs,其中上调DEGs共有74个,下调DEGs共有44个,表达模式相反的DEGs共有24个(表2)。结果表明,耐高温自交系响应高温胁迫的基因数量更多,且有更多的DEGs表现为上调表达。上调和下调DEGs数目的差异说明不同耐高温玉米自交系籽粒响应高温胁迫的基因存在较大差异,可能涉及一系列基因参与的生理生化过程。

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

参考生物学过程、细胞组分和分子功能3个大类的基因GO功能注释,对高温胁迫下玉米籽粒的DEGs 进行GO富集分析,XN202和CT110分别得到353,341个DEGs被显著富集。在富集显著的前 20个GO功能注释分类中,SCK-1-2-3 vs SHT-1-2-3 对比组和TCK-1-2-3 vs THT-1-2-3 对比组中的DEGs分别富集在11,6个生物学过程,3,9个细胞组分,6,5个分子功能(图2)。

A.CT110;B.XN202。

富集在生物学过程的亚类中,SCK-1-2-3 vs SHT-1-2-3 对比组的DEGs主要集中在蛋白质折叠(GO:0006457)、热激反应(GO:0009408),分别有39,29个基因功能得到注释;TCK-1-2-3 vs THT-1-2-3 对比组的DEGs主要集中在翻译(GO:0006412),共有67个基因功能得到注释;DNA复制起始(GO:0006270)是二者共有的显著富集代谢通路,说明在玉米籽粒响应高温胁迫的生物学过程中DNA复制起始发挥重要作用。

富集在细胞组分的亚类中,SCK-1-2-3 vs SHT-1-2-3 对比组的DEGs主要集中在核小体(GO:0000786),共有28个基因功能得到注释;TCK-1-2-3 vs THT-1-2-3 对比组的DEGs主要集中在核糖体(GO:0005840)和细胞质核糖体大亚基(GO:0022625),分别有37,35个基因功能得到注释;二者共有3个显著富集代谢通路,分别为核小体(GO:0000786)、微小染色体维持复合物(GO:0042555)、DNA引物酶-多聚酶α复合体(GO:0005658)。

富集在分子功能的亚类中,SCK-1-2-3 vs SHT-1-2-3 对比组的DEGs主要集中在未折叠蛋白结合(GO:0051082)、蛋白质异二聚活性(GO:0046982),分别有42,28个基因功能得到注释;TCK-1-2-3 vs THT-1-2-3 对比组的DEGs主要集中在核糖体的结构组成(GO:0003735),共有85个基因功能得到注释;营养库活性(GO:0045735)是二者共有的显著富集代谢通路。

2.5 差异表达基因的KEGG富集分析

为进一步解析玉米自交系籽粒响应高温胁迫的代谢调控通路,分析DEGs可能参与的生物学功能,分别对不同耐高温自交系响应高温胁迫的DEGs进行KEGG代谢通路富集分析,XN202和CT110分别得到176,138个DEGs 被显著富集(表3)。

表3 差异表达基因的KEGG注释

SCK-1-2-3 vs SHT-1-2-3 对比组的DEGs主要集中在内质网蛋白质加工(ko04141)、DNA复制(ko03030)、同源重组(ko03440)、错配修复(ko03430)、丙氨酸、天冬氨酸和谷氨酸代谢(ko00250)、嘌呤代谢(ko00230)、谷胱甘肽代谢(ko00480)等途径。有52个DEGs富集到内质网蛋白质加工是KEGG富集分析最显著的通路。

TCK-1-2-3 vs THT-1-2-3 对比组的DEGs主要集中在核糖体(ko03010)、DNA复制(ko03030)、丙氨酸、天冬氨酸和谷氨酸代谢(ko00250)、精氨酸生物合成(ko00220)、乙醛酸和二羧酸代谢(ko00630)、真核生物核糖体的生物合成(ko03008)等途径。有89个DEGs富集到核糖体是KEGG富集分析最显著的通路。

二者共有2个显著富集代谢途径,分别为DNA复制和丙氨酸、天冬氨酸和谷氨酸代谢。

2.6 差异表达基因的COG功能注释分析

对不同耐高温玉米自交系籽粒响应高温胁迫的DEGs进行COG同源分类,XN202和CT110分别得到540,365个DEGs有COG功能注释信息,分别属于24,21个功能类别(图3),其中XN202的DEGs主要分布于翻译,核糖体结构和生物起源(98个DEGs)、翻译后修饰,蛋白转运(57个DEGs)、一般功能预测(46个DEGs)、氨基酸转运和代谢(44个DEGs)、碳水化合物转运和代谢(41个DEGs)、信号传导机制(34个DEGs)、次生代谢产物的合成,转运和代谢(31个DEGs)、复制,重组和修复(30个DEGs)、防御机制(26个DEGs)、细胞壁/细胞膜/包膜的生物发生(24个DEGs)10个主要分类(富集DEGs个数≥20个);CT110的DEGs主要分布于翻译后修饰,蛋白转运(81个DEGs)、碳水化合物转运和代谢(37个DEGs)、信号传导机制(33个DEGs)、复制,重组和修复(29个DEGs)、一般功能预测(26个DEGs)、氨基酸转运和代谢(20个DEGs)6个主要分类(富集DEGs个数≥20个)。COG注释结果在宏观上解析了不同耐高温玉米自交系籽粒响应高温胁迫相关基因的功能分布特征。

图3 COG功能分类

2.7 玉米自交系籽粒响应高温胁迫的转录因子分析

在鉴定的2 387个DEGs中,共有49个转录因子家族的125个转录因子,主要有AP2/ERF(12个)、AUX/IAA(3个)、bHLH(7个)、bZIP(5个)、C2H2(5个)、GARP-G2-like(5个)、HMG(4个)、HSF(6个)、MADS(5个)、MYB(9个)、NAC(9个)、WRKY(3个)等参与植物生长发育和逆境胁迫响应相关的重要转录因子。

TCK-1-2-3 vs THT-1-2-3 对比组与SCK-1-2-3 vs SHT-1-2-3 对比组检测到142个共同的DEGs中有7个转录因子,其中共同上调的转录因子5个(GRMZM2G030768、GRMZM2G033828、GRMZM2G062657、GRMZM2G125648、GRMZM2G395749),共同下调的转录因子2个(GRMZM2G011110和GRMZM2G042756均为AP2/ERF转录因子),这些基因可能在玉米籽粒响应高温胁迫过程中发挥重要作用。

2.8 差异表达基因与耐高温QTL区间的联合分析

将DEGs与已报道的玉米耐高温相关性状QTL定位和全基因组关联分析定位结果进行联合分析,共有374个DEGs与已有研究报道的QTL区间共定位,其中第1染色体上有31个,第2染色体上有117个,第3染色体上有32个,第4染色体上有11个,第5染色体上有60个,第6染色体上有8个,第7染色体上有2个,第8染色体上有11个,第9染色体上有73个,第10染色体上有29个。

位于玉米耐高温相关性状QTL的374个DEGs中,有42个为已报道的耐高温相关候选基因(表4)。谷胱甘肽S移换酶和RWP-RK转录因子在SCK-1-2-3 vs SHT-1-2-3 对比组和TCK-1-2-3 vs THT-1-2-3 对比组都上调表达,而水通道蛋白和RLK-Pelle_DLSV蛋白激酶在SCK-1-2-3 vs SHT-1-2-3 对比组和TCK-1-2-3 vs THT-1-2-3对比组都下调表达。值得注意的是HSP70、HSP20和HSP101家族基因都在高温敏感自交系CT110中显著下调表达,但是在耐高温自交系XN202中差异不显著,可能是2个自交系耐高温能力差异形成的关键基因。

表4 玉米耐高温QTL区间内的候选基因

3 结论与讨论

由于玉米耐热性是一种受多基因控制的复杂数量性状,更容易受到基因型-环境互作的影响,传统的农艺栽培技术对缓解高温胁迫的作用较小,因此,提高玉米的耐高温能力是应对全球气候变化最经济有效的方法。玉米对高温胁迫的响应,主要涉及信号感知、信号转导、细胞结构、生理生化、基因表达和调控等过程[10]。本研究通过对授粉后15 d的玉米籽粒转录组的测序数据分析,在2个耐高温能力不同的玉米自交系中获得大量的DEGs,这些基因的表达和调控涉及高温胁迫下生理生化过程和代谢途径的变化,说明玉米耐热性与基因的差异表达有关。耐高温自交系XN202中的DEGs明显多于高温敏感自交系CT110,且XN202的DEGs中上调基因比下调基因多,而CT110的DEGs中下调基因比上调基因多,这些上调和下调DEGs数目的差异可能影响玉米籽粒发育时期的耐热性。XN202和CT110基因组中共有响应高温胁迫的DEGs只有142个,可见不同基因型的玉米对高温胁迫响应的分子机制存在明显差异[15]。

通过对DEGs的GO功能注释分析,XN202和CT110分别得到353,341个DEGs被显著富集,其中生物学过程(DNA复制起始)、细胞组分(核小体、微小染色体维持复合物、DNA引物酶-多聚酶α复合体)、分子功能(营养库活性)共同参与响应高温胁迫,增强了玉米籽粒对热胁迫的响应。这表明高温胁迫主要影响玉米的生物学过程和分子功能[8]。XN202有85个DEGs显著富集在核糖体的结构组成(分子功能),67个DEGs显著富集在翻译(生物学过程),说明这些通路在玉米籽粒遭受高温胁迫时可能发挥重要作用。谷氧还蛋白是维持活性氧平衡的重要组分,在玉米中异源表达拟南芥谷氧还蛋白家族基因AtGRXS17显著影响核糖体RNA生物合成,可以提高玉米的耐热性[16]。CT110有42个DEGs显著富集在未折叠蛋白结合(分子功能),39个DEGs显著富集在蛋白质折叠(生物学过程)。未折叠蛋白反应是一种与内质网应激相关的细胞应激反应,有助于减轻热胁迫对细胞造成的损害[17]。通过对DEGs的全基因组及代谢途径数据库(KEGG)信号通路富集分析,XN202和CT110分别得到176,138个DEGs被显著富集,只有DNA复制,丙氨酸、天冬氨酸和谷氨酸代谢等共同参与响应高温胁迫的应答过程。XN202有89个DEGs显著富集在核糖体路径,CT110有52个DEGs显著富集在内质网蛋白质加工路径,进一步说明核糖体和内质网蛋白质加工代谢路径可能是2个自交系耐高温能力差异形成的重要原因。将DEGs与COG数据库比对,XN202和CT110分别得到540,365个DEGs有COG功能注释信息,其中XN202有98个DEGs与翻译、核糖体结构和生物起源有关,CT110有81个DEGs与翻译后修饰、蛋白转运有关,XN202和CT110分别有34,33个DEGs与信号传导机制有关,这些功能基因可能与玉米耐热性有关,值得进一步研究。

转录因子可以通过调控下游多个靶基因的表达,参与植物热激响应的调控网络,是响应高温胁迫的关键基因和重要分子开关[18-19]。本研究中DEGs共有49个转录因子家族的125个转录因子,其中AP2/ERF、MYB、bHLH、NAC、HSF等家族转录因子都直接或间接参与高温胁迫响应[10]。为了减轻高温胁迫对细胞的损害,植物需要HSF来快速转录激活热激蛋白(HSP),与耐高温能力密切相关[20]。HSP作为一种重要的分子伴侣,可以参与修复受损的蛋白质,减少蛋白质的错误折叠和聚集[17]。位于玉米耐高温相关性状QTL的374个DEGs中,有42个为已报道的耐高温相关候选基因,其中15个基因的功能未知,共有7个HSP。HSP70、HSP20和HSP101家族基因都在高温敏感自交系CT110中均显著下调表达,但是在耐高温自交系XN202中差异不显著,说明在高温胁迫条件下HSPs的积累量可能与耐热性正相关。

本研究选择以授粉后15 d的玉米籽粒为试验材料,基于转录组高通量测序技术,发现高温胁迫下玉米籽粒会形成复杂的细胞保护和防御系统,不同耐高温自交系的籽粒对高温胁迫的响应存在明显差异,AP2/ERF、MYB、bHLH、NAC、HSF、HSP等耐高温相关的DEGs可能在分子调控网络中发挥重要作用。研究结果为进一步解析玉米耐高温的分子遗传机制,克隆玉米耐高温关键基因,开发分子标记,筛选或培育耐高温玉米新品种奠定了很好的理论基础。