玉米CMS-S 型不育系花药的转录组分析

闫嘉颖，王久红，赵永锋，贾晓艳，郭晋杰，祝丽英

（河北农业大学 农学院/国家玉米改良中心河北分中心/河北省作物种质资源实验室，河北 保定 071001）

植物雄性不育是指植物在有性生殖过程中，雄性器官发育异常，不能产生正常功能的花药、花粉或雄配子的现象［1］。在棉花、小麦、大白菜等多种作物中，发现雄性不育系花药短小干瘪，成熟期花药不开裂，药囊内无花粉或花粉没有活力［2］。作物育性的形成依赖于花药正常发育，而花药发育过程受大量基因的精密调控［3］，任何发育时期的功能异常都会导致花粉败育。

转录组测序技术是对特定细胞在某一状态下转录出来的所有RNA 进行测序分析，能够反映细胞中基因的表达情况及其调控规律，在差异表达基因（Differentially expressed genes，DEGs）筛选、新基因挖掘、分子标记开发、代谢网络等方面的研究中广泛应用［4］。目前，多种作物雄性不育系的转录组测序工作已相继开展。曹晓聪等［5］在转录组水平上，分析了影响棉花CMS 的基因表达模式，探讨与CMS 相关的一些潜在基因和代谢途径。闫东等［6］以紫花苜蓿雄性不育系MS-JN1A 及保持系MS-JN1B 为材料，对5 个发育时期的花药进行转录组测序，筛选DEGs 进行富集分析和功能注释，探讨与紫花苜蓿雄性不育相关的代谢通路。Li 等［7］比较了CMS-C 型不育系C48-2 及其保持系N48-2在花粉母细胞和单核期的转录结果，检测到2 069个与小孢子败育过程有关的DEGs，讨论与玉米CMS-C 型不育相关的基因和代谢途径。

玉米是杂种优势利用最早也是最成功的作物之一，雄性不育在玉米杂交制种和杂种优势利用中发挥了重要作用［8］。CMS-S 型不育系是玉米雄性不育中资源最丰富的一类，虽然目前对其已有一些细胞形态和生理、生化机制方面的研究，但至今尚未完全阐明其不育机制，仍需从分子层面进行深入研究。郑58cms-Q1261 是河北农业大学玉米育种课题组选育的CMS-S 型雄性不育系，其不育特性稳定，通过细胞学观察发现不育系小孢子败育的关键时期是单核晚期［9］。本研究以不育系郑58cms-Q1261和保持系郑58 单核晚期的花药为研究对象，利用转录组测序和生物信息学分析方法鉴定DEGs，通过GO 和KEGG 富集分析DEGs 参与的生物学过程和代谢途径，为CMS-S 型雄性不育分子机理的深入研究提供参考。

1 材料与方法

1.1 试验材料

选用河北农业大学玉米育种课题组以玉米自交系郑58 为背景构建的CMS-S 型不育系郑58cms-Q1261 和保持系郑58。

试验材料于2021 年5 月播种于河北农业大学农业科技园，播种后50 d 起，取不育系和保持系的花药，利用涂片法镜检，确认小孢子的发育时期。待发育至单核晚期进行取样，每个样品由多株花药组成，3次重复。保持系郑58 和不育系郑58cms-Q1261 分别命名为A 和B，3 次重复分别用A1、A2、A3 和B1、B2、B3 表示，液氮速冻后-80 ℃保存备用。

1.2 扫描电镜观察

取单核晚期的花药置于2.5%的戊二醛固定液中固定。分别用50%、75%、85%、90%、95% 乙醇逐级脱水，每个梯度30 min，最后用无水酒精脱水重复2～3 次。在预冷的醋酸异戊酯置换20 min，临界点干燥。将干燥样品固定在样品台上，用离子溅射法对材料进行喷金镀膜。最后用扫描电镜（日立SU8100）进行观察。

1.3 花药RNA 提取

使用植物总RNA 提取试剂盒（华越洋）提取花药样品的总RNA，使用NanoDrop 2000 (Thermo Fisher Scientific,Wilmington,DE) 对RNA 纯度和浓度进行检测。利用1%的琼脂糖凝胶电泳（SDSPAGE）和Agilent 生物分析仪2100 系统的RNA Nano 6000 检测试剂盒对RNA 完整性进行评估。

1.4 转录组测序及DEGs 分析

使用Illumina 测序平台对样品测序并分析，原始数据经过处理后，得到clean reads 用于后续分析。利用HiSat2 软件将过滤后的数据比对到玉米参考基因组（ZmB73_RefGen_v4），利用stringtie 软件对转录本拼接并估算基因表达量。对样本FPKM 值进行统计，利用R 软件包DESeq2 对基因表达值进行DEGs 分析。以基因相对表达量的差异倍数|log2(fold change)| ＞1.5 且错误发现率（False discovery rate，FDR）＜ 0.05 为DEGs 筛选的阈值标准。

1.5 qRT-PCR 验证

选取12 个DEGs 进行实时定量PCR 验证。采用TRIZOL 法提取总RNA，合成相应的cDNA。利用Premier5.0 software 设计荧光实时定量PCR 引物，引物序列如表1 所示。以玉米肌动蛋白基因作为内参基因，计算相对表达量（2-ΔΔCt法），保持系A中各目标基因中的表达量为1，不育系B 目标基因中的表达量与之进行比较。

表1 qRT-PCR 中所用差异表达基因及其引物序列Table 1 Lists the DEGs and their primer sequences used in qRT-PCR

2 结果与分析

2.1 雄穗和花药特征

观察散粉期保持系A 和不育系B 的雄穗和花药，如图1 所示，发现保持系雄穗发达，颖壳开裂，花药正常散粉，花粉量大；不育系雄穗外形和保持系无明显差异，但颖壳不开裂，花药不外露。进一步对保持系和不育系的小花进行观察，保持系花药饱满，而不育系花药干瘪。

图1 郑58 和郑58cms-Q1261 的雄穗和花药Fig. 1 Tassel and anthers of Zheng 58 and Zheng 58cms Q1261

利用扫描电镜（SEM）对保持系A 和不育系B单核晚期的花药进行观察（图2）。可见保持系花药饱满（图A-1），花药外壁的表皮细胞排列相对疏松（图A-2），花药内壁表面纤维层存在突起，其表面存在乌氏小体（图A-3），药室内有大量圆球状饱满的花粉粒（图A-4，A-5）。不育系花药干瘪（图B-1），花药外壁表皮细胞排列相对紧密（图B-2），内壁较为平滑，其表面乌氏小体颗粒数量、大小、分布与保持系相似（图B-3），药室内花粉粒皱缩，形状不规则（图B-4，B-5）。

图2 郑58 和郑58cms-Q1261 花药的扫描电镜观察Fig. 2 Sem observation of anther of Zheng58 and Zheng58 cms-Q1261

图3 差异表达基因火山图Fig. 3 Volcano map of DEGs

2.2 转录组学分析

2.2.1 测序质量分析 利用Illumina HiSeq 2500 测序平台对保持系A 和不育系B 单核晚期的花药3 次重复样品进行转录组测序。所测数据经过严格的质控后，共得到41.67 Gb 的序列数据，各样品Clean Data 均达到6.35 Gb。如表2 所示，获得序列的GC含量在55.75%～56.56%，各样品Q20 碱基百分比均不低于98.00%，Q30 碱基百分比在94.32%以上。表明样品测序质量较高，能够用于后续的分析。

表2 转录组测序质量分析Table 2 Transcriptome sequencing quality analysis

2.2.2 差异表达基因分析 基于基因的表达量，以|log2(fold change)| ＞1.5 且FDR＜ 0.05 为标准筛选差异表达基因。通过对保持系和不育系单核晚期花药的转录本的比较分析，共鉴定到14 565 个DEGs，其中在不育系郑58 cms-Q1261 中7 551 个基因上调表达，7 014 个基因下调表达（图 3）。

2.2.3 差异表达基因GO 分析 对上述筛选到的DEGs 进行GO 富集分析，共有11 232 个DEGs 注释到生物学过程、细胞组分和分子功能3 个部分(图4)。有8 440 个DEGs 注释到生物学过程，涉及生物过程中的22 个功能条目，其中代谢进程（上调DEGs 有3 079 个，下调有DEDGs 2 634 个）；细胞进程（上调DEGs 有2 944 个，下调有DEDGs 2 643 个）和单组织过程（上调DEGs 有2 350 个，下调有DEDGs 1 930 个）等二级功能的基因富集比例较大。有8 539 个DEGs 注释到细胞组分，涉及16 个功能条目，细胞（上调DEGs 有3 845 个，下调有DEDGs 3 499 个）；细胞组成（上调DEGs 有3 850 个，下调有DEDGs 3 504 个）和细胞器官（上调DEGs 有3 073 个，下调有DEDGs 2 995 个）这3 个二级功能占比较大。有8 568 个DEGs 注释到分子功能的15 个条目，其中结合功能（上调DEGs 有3 292 个，下调有DEDGs 2 364 个）和催化活性（上调DEGs 有2 898 个，下调有DEDGs 2 272 个）所占比例较高。

2.2.4 差异表达基因KEGG 分析 为了进一步了解DEGs 的生物学功能，将鉴定到的DEGs 进行KEGG 分析。结果发现，DEGs 共富集到127 条代谢通路，表明玉米CMS-S 型不育系花粉败育是多通路协调调控的复杂过程。分析显著富集的前20 条通路（图5），可见DEGs 主要富集到淀粉和蔗糖代谢、碳代谢、氧化磷酸化、糖酵解/糖异生、氨基糖和核苷酸糖代谢和植物-病原互作，而脂肪酸生物合成和糖胺多糖的降解富集较少。

图5 差异表达基因 KEGG Pathway 富集气泡图Fig. 5 Bubble chart of KEGG Pathway enrichment of DEGs

2.3 差异表达基因qRT-PCR 验证

根据基因富集结果选取12 个DEGs 进行qRTPCR 验证，详细信息见表3，其中包括与戊糖和葡萄糖醛酸的相互转化(ko00040)相关的2 个果胶裂解酶超家族蛋白；与淀粉和蔗糖代谢(ko00500)相关的5 个果胶裂解酶类；与氨基糖和核苷酸糖代谢(ko00520)相关的2 个bHLH 转录因子；与内质网蛋白加工(ko04141)相关的3 个热休克蛋白。qRTPCR 数据通过2-ΔΔCT换算为相对表达量（Relative expression），并进行独立t检验。图6 结果显示，这些 DEGs 在不育系B 和保持系A 中表达量差异极显著，qRT-PCR 验证结果与 RNA-Seq 的差异表达趋势一致，说明本研究的 RNA-Seq 测序结果可靠，推测这些基因可能在花粉发育过程中发挥重要的调控作用。

图6 12 个差异基因的RT-qPCR 验证Fig. 6 RT-qPCR validation of 12 DEGs

3 讨论

前人关于不同作物雄性不育转录组的研究发现不育过程常常伴有大量基因参与，不同作物的DEGs富集到的主要基因功能和代谢途径不尽相同。魏霁桐等［10］对小麦生理型雄性不育系PHYMS-1376花药5 个时期的转录组分析，显示有36 058 个DEGs，主要富集在淀粉和蔗糖代谢、苯丙素生物合成和植物激素信号传导等通路。薛亚东等［11］在玉米CMS-C 型“三系”材料减数分裂前期Ⅰ、中期Ⅰ及四分体时期花药的转录组分析中，筛选出DEGs 5 676 个，主要富集在氧化磷酸化、碳代谢及糖酵解等能量代谢相关的途径。Wei 等［12］关于棉花雄性不育的转录组研究结果显示，在减数分裂和四分体阶段分别有9 595 个、10 407 个DEGs，主要富集在蔗糖和淀粉代谢、激素合成和糖酵解磷酸戊糖途径合成等通路。本研究对CMS-S 型不育系郑58cms-Q1261 和保持系郑58 单核晚期花药进行RNA-seq 测序，结果显示有14 565 个DEGs，显著富集在淀粉和蔗糖代谢、碳代谢、氧化磷酸化、糖酵解/糖异生、氨基糖和核苷酸糖代谢等途径中。Xing 等［13］在油菜雄性不育的转录组分析中，发现DEGs 主要参与淀粉和蔗糖代谢、苯丙烷生物合成。Lui 等［14］关于烟草雄性不育研究中，DEGs 主要富集在内质网蛋白加工和氧化磷酸化。高焕庭等［15］在小麦雄性不育研究中，KEGG 代谢途径主要有糖酵解途径和三羧酸循环。与本研究中富集到的代谢途径部分一致，表明这些代谢通路可能与花粉败育相关。进一步说明不同物种雄性不育均受众多基因影响，是一个多通路协同调控过程，育性调控机理复杂。

淀粉和蔗糖代谢为小孢子发育提供营养，糖类、淀粉等物质的缺失可能造成代谢紊乱，影响小孢子的形成，导致植物雄性不育［16］。Su 等［17］采用RNA-seq 技术，在玉米CMS-S 型雄性不育中鉴定出2 108 个DEGs，这些DEGs 在能量产生和转化、碳水化合物代谢和信号转导等生物过程中发挥作用，进一步研究了淀粉和蔗糖代谢的途径，找到一些与蔗糖生物合成过程有关的DEGs。果胶是花粉细胞壁的主要组分，在细胞壁结构、信号转导和植物生长发育中具有重要作用［18］，果胶降解是淀粉与蔗糖的代谢的一个重要分支。本研究选取与果胶酶相关的淀粉蔗糖代谢通路进行分析，发现编码果胶降解相关的一系列酶（包括2 个果胶裂解酶超家族蛋 白Zm00001d046389、Zm00001d053597 和5 个果胶裂解酶类Zm00001d015820、Zm00001d036820、Zm00001d015821、Zm00001d015825、Zm00001d036824）的DEGs 均显著下调，推测在不育系中由于果胶裂解酶基因的缺失导致果胶异常积累，影响花粉的正常成熟进而导致雄性不育。与Li 等［19］人在小麦花药和花粉发育调控相关基因的研究中，果胶降解相关基因的下调导致花药不开裂的结果一致。

转录因子对植物育性具有一定的调控作用，在拟南芥［20］、水稻［21］和玉米［22］等作物中能够控制花药分裂和雄性育性。其中bHLH 转录因子在所有转录因子家族中数量最多，已有研究表明，bHLH转录因子与花药绒毡层和小孢子发育相关，其功能异常时，绒毡层细胞提前降解，花粉形成受阻导致小孢子发育受阻，最终表现败育［23］。Jung 等［24］在水稻中鉴定了1 个与绒毡层发育相关的bHLH 转录因子UDT1，其逆转录转座子Tos17 插入形成的不育植株表现雄蕊周缘细胞分化异常，绒毡层液泡化，花粉粒无法形成。Dukowic 等［25］比较拟南芥和玉米在减数分裂时期的转录组发现，14 个拟南芥bHLH基因和3 个玉米bHLH基因上调表达，并且发现玉米和拟南芥有对应的同源基因。Nakata 等［26］在拟南芥中鉴定出参与植物激素信号传导的3 个bHLH 转录因子（JAM1、JAM2 和JAM3），发现它们在雄蕊发育中作为负调控因子发挥作用。本研究也鉴定到2 个bHLH 转录因子（Zm00001d021019、Zm00001d038357）在单核晚期的不育系中上调表达。

热休克蛋白是植物在逆境条件下产生的应激蛋白，能够平衡植物体内的代谢，在生长发育中具有重要的作用［27］。有研究表明，在没有逆境条件下，一些热休克蛋白与配子体的发育相关，TMS1是拟南芥热敏雄性不育基因，敲除TMS1 基因后发现在30℃条件下花粉管生长受阻，导致雄性生育能力显著降低［28］。Mei 等［29］在萝卜细胞质雄性不育研究中，鉴定出5 个热休克转录因子(HSFs)和19 个热休克蛋白(HSPs) 基因在萝卜CMS 系中下调表达，对小孢子发育产生影响。本研究中有3 个热休克蛋白（Zm00001d012420、Zm00001d028557、Zm00001d020898）在不育系中显著高表达，推测是其他代谢途径紊乱导致热休克蛋白高度表达，影响蛋白质加工或降解，在一定程度上抑制其他种类的蛋白质的合成，从而影响小孢子发育。

本研究对不育系郑58cms-Q1261 和保持系郑58 单核晚期花药进行转录组测序，筛选出14565 个DEGs，GO 和KEGG 分析DEGs 主要富集在淀粉和蔗糖代谢、碳代谢和氧化磷酸化等代谢通路，发现果胶裂解酶类、热休克蛋白、bHLH 转录因子等与雄性不育有关，为进一步深入解析玉米S 型细胞质雄性不育的分子机制提供了理论基础。