玉米花药不同发育时期的转录组分析及其差异表达基因的鉴定

陈慧， 孙晓靖， 郭瑶晴， 连玉杰， 张雪海， 汤继华， 陈晓阳

(河南农业大学农学院，河南 郑州 450046)

雄性育性影响玉米种子形成和生存繁衍，雄性不育系在玉米杂交制种和杂种优势利用中发挥了重要作用[1-2]。雄性育性的形成依赖于花药正常发育，其发育过程受一系列遗传因子精密调控。依据花药发育的细胞学特征，可将玉米花药发育划分为14个时期，主要分为减数分裂前期(S1～S6)，减数分裂期(S7～S8)和减数分裂后期(S9～S14)3个阶段[3]。其中，S6期花药内形成小孢子母细胞，开始进入减数分裂阶段；S8b期花药完成减数分裂过程，形成四分体细胞；S13期花药已形成成熟花粉粒，包裹花药的内外稃张开。

花药角质层和花粉壁是雄性育性形成过程中重要屏障。花药角质层覆盖于花药表面，防止各种生物胁迫和非生物胁迫破坏花药发育[4]。花药角质层主要由角质和蜡质组成，其中角质包含C16和C18脂肪酸、羟基脂肪酸、环氧脂肪酸、二元脂肪酸等物质[5]，蜡质由脂肪酸、烃类、醇类、酯类等物质构成[6]。花粉壁位于小孢子(花粉)表面，保护小孢子正常发育和成熟花粉粒形成[7]。利用雄性不育突变体已鉴定到部分参与玉米花药角质层和花粉壁发育的关键基因。例如，脂肪酸羟化酶基因APV1和MS26[8-9]，羟基脂肪酸脱氢酶基因IPE1/ZmMs20[10-11]，脂肪酰还原酶基因MS6021/ZmMs25/ZmFAR1[12-14]，GDSL脂肪酶基因ZmMs30和IPE2[15-17]，甘油-3-磷酸酰基转移酶基因ZmMs33[18]，脂类物质转蛋白基因MS2/ZmABCG26[14,19]，转录调控基因ZmMs7，ZmbHLH51，ZmbHLH122，ZmMYB84，ZmMYB33-1/2和ZmPHD11[20-22]。

转录组测序技术(RNA-Seq)是利用高通量测序技术对生物进行转录组分析,是一种检测全基因组基因表达谱的有效方法[23]。通过比较不同发育时期基因表达谱，可以挖掘到影响植物生长发育的关键基因。WAN等[24]利用玉米自交系Oh43和W23花药转录组数据，鉴定到125个参与雄性育性形成的脂类物质代谢相关基因。JIANG等[22]分析3个自交系S5～S12时期花药转录谱，发现1 100个转录因子在3个自交系中均表现出时期差异表达，进一步通过CRISPR/Cas9技术验证了12个转录因子的功能。然而，目前尚未利用转录组测序技术挖掘影响玉米花药角质层和花粉壁发育的关键基因。本研究对玉米自交系B73的3个发育时期花药和苗期叶片进行转录组测序，通过GO和KEGG等生物学信息分析鉴定影响花药角质层和花粉壁形成的关键基因及其参与的代谢途径，为玉米雄性育性调控机制的解析奠定基础，同时为雄性不育系的创制提供潜在靶标位点。

1 材料与方法

1.1 试验材料

供试的玉米自交系为B73，2020年4月种植于河南省新乡市原阳县河南农业大学原阳试验基地(35.1°N，113.6°E)。依据花药长度和碘-碘化钾染色判断花药发育时期，选取S6期(记为B1)，S8b期(记为B2)，S13期(记为B3)花药，以及苗期叶片(记为B4)于-80 ℃保存。每组样品包含4个生物学重复。

1.2 RNA提取、文库构建和测序

利用天根试剂盒提取样品总RNA。Nanodrop检测RNA纯度，琼脂糖凝胶电泳分析RNA降解程度及是否有污染，Qubit2.0对RNA浓度进行精确定量，Agilent 2100精确检测RNA完整性。样品检测合格后，按如下流程构建文库：(1)利用带有Oligo(dT)磁珠富集mRNA。(2)加入fragmentation buffer将mRNA打断成短片段。(3)以mRNA为模板，用六碱基随机引物合成第1条cDNA链，加入缓冲液、dNTPs、RNase H和DNA polymerase I合成第2条cDNA链。(4)双链cDNA末端修复且在3’端加A。(5)QiaQuick PCR试剂盒纯化并加EB缓冲液洗脱，末端修复，连接测序接头。(6)琼脂糖凝胶电泳选择片段大小，PCR扩增富集cDNA。文库质检合格后，通过Illumina Novaseq平台进行PE150测序。

1.3 测序数据过滤、参考基因组比对和基因表达定量

测序得到原始序列raw reads，通过数据过滤得到干净、有效、高质量的clean reads。数据过滤步骤：(1)当单端测序reads中N>3时，去除此对paired reads；(2)当单端测序read中质量值低于5的碱基比例≥20%时，去除此对paired reads；(3)去除adapter序列时，adapter序列至少要匹配8 bp。采用Hisat2软件把过滤后的reads与玉米自交系B73参考基因组(V4版本)比对[25]。利用featureCounts软件对每个样品进行基因水平的定量，合并得到所有样本的基因表达量矩阵[26]。

1.4 差异表达基因鉴定和功能富集分析

利用edge R软件分析基因在各样品中的差异表达，计算差异表达的Pvalue和Padj值，Padj是校正后的Pvalue[27]。差异表达基因筛选条件为：Padj<0.05，|log2FoldChange(FC)|>1。应用topGO软件进行差异表达基因GO(gene ontoloy)富集分析。应用KOBAS软件进行差异表达基因KEGG(kyoto encyclopedia of genes and genomes)通路富集分析[28]。

1.5 qRT-PCR验证

通过NCBI(https://www.ncbi.nlm.nih.gov/)设计qRT-PCR引物，引物序列见表1。利用2×SYBR Green qPCR Premix(Universal)试剂盒在伯乐C1000 TouchTMThermal Cycler仪器进行qRT-PCR分析。以ZmACTIN为内参，2-ΔΔCt法计算相对表达量。

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

2 结果与分析

2.1 玉米花药发育不同时期转录组测序数据的质量分析

通过Illumina Novaseq平台对16个RNA样品文库进行测序，测序共产生456 100 906条clean reads，质量值≥30的碱基所占比例(Q30)均大于90%，说明测序数据质量较好，可用于进一步分析。将获得的有效测序数据比对到玉米自交系B73参考基因组(V4版本)上，比对率为92.46%～96.49%，平均比对率为95.56%(表2)。

表2 玉米花药发育不同时期RNA-Seq数据的统计分析

为验证测序数据的可靠性，分析了8个玉米雄性育性基因在测序数据中的表达模式。MS32在玉米花药发育早期表达，调控绒毡层分化过程[29]。在测序数据中，MS32基因在S6期花药中特异性表达。IPE1/ZmMs20，APV1，IPE2，MS2/ZmABCG26，ZmMs7，ZmMYB84调控花药角质层和花粉外壁发育，主要在S8b和S9期花药表达[8，10-11，14，17，19-22]。RNA-Seq结果表明，6个基因在S8b期花药中表达量最高。LBD10被报道影响花粉后期发育，通过检测发现，该基因在S13时期花药表达量最高[22]。上述结果表明，8个基因在测序数据中的表达模式与已报道基因表达模式一致。同时，qRT-PCR分析显示，8个基因RNA-Seq和qRT-PCR结果吻合，说明测序数据真实可靠(图1)。

FPKM:每千个碱基的转录每百万映射读取的碎片。FPKM: Fragments per kilobase of exon model per million mapped fragments.

2.2 玉米花药发育不同时期差异表达基因分析

以Padj<0.05和|log2FC|>1作为差异表达基因的筛选标准，在B1 vs B4，B2 vs B4，B3 vs B4，B2 vs B1，B3 vs B1，B2 vs B3 6组之间分别鉴定到15 900，14 785，17 141，8 115，17 858，16 497个差异表达基因(表3)。不同发育时期花药与苗期叶片相比，S13期花药与苗期叶片差异表达基因数最多(17 141)，其次是S6期(15 900)，S8b期最少(14 785)。不同发育时期花药相比，S6与S13期花药差异表达基因最多(17 858)，其次是S8b与S13(16 497)，S6与S8b最少(8 115)，说明S6与S8b期花药转录谱相似。此外，S6与S13期花药差异表达基因多于任一时期花药与苗期叶片差异表达基因，表明不同发育时期花药基因表达差异较大。

表3 3个时期花药和苗期叶片差异表达基因数目

2.3 玉米花药角质层发育相关基因的鉴定

为挖掘影响玉米花药角质层发育的候选基因，通过KEGG方法分析不同发育时期花药差异表达基因参与的生化代谢途径。B2 vs B1和B2 vs B3上调差异表达基因均富集到角质和蜡质合成代谢途径，而两组下调差异表达基因及B1 vs B3差异表达基因均不能富集到角质和蜡质合成代谢途径(图2)。这表明S8b时期是玉米花药角质层发育的关键时期，该时期相关基因的表达对角质层形成至关重要。

图2 B2 vs B1(A)和B2 vs B3(B)上调差异表达基因KEGG富集分析

B2 vs B1中8个差异表达基因富集到角质和蜡质合成代谢途径。Zm00001d049950和Zm00001d048337编码脂肪酰还原酶，Zm00001d012274和Zm00001d042723编码ω羟基棕榈酸阿魏酸酰基转移酶，Zm00001d025833和Zm00001d017953编码过加氧酶，Zm00001d042814和Zm00001d012326编码细胞色素单加氧酶CYP86A1(表4)。

表4 B2 vs B1上调差异表达基因中角质和蜡质合成相关基因

B2 vs B3中15个差异表达基因富集到角质和蜡质合成代谢途径，其中Zm00001d049950，Zm00001d048337，Zm00001d012274，Zm00001d042723，Zm00001d025833，Zm00001d042814和Zm00001d012326在B2 vs B1中被鉴定到。Zm00001d014055，Zm00001d051932，Zm00001d017251编码醛脱羧酶，Zm00001d032284和Zm00001d020238编码ω羟基脂肪酸脱氢酶，Zm00001d051800编码过加氧酶，Zm00001d017528和Zm00001d002687分别编码细胞色素单加氧酶CYP86A35和CYP86A7(表5)。

表5 B2 vs B3上调差异表达基因中角质和蜡质合成相关基因

2.4 玉米花粉壁发育相关基因的鉴定

花粉壁位于小孢子外表面，其形成起始于四分体时期的初生外壁。为鉴定影响花粉壁发育的关键基因，分析了四分体时期花药与其他组样品之间(B2 vs B1，B2 vs B3，B2 vs B4)的差异表达基因，共鉴定到3 555个共同差异表达基因(图3)。其中，1 348个基因在S8b时期花药中均上调表达，344个基因均下调表达。

图3 B2 vs B1，B2 vs B3，B2 vs B4差异表达基因韦恩图

GO分析表明，1 348个共同上调差异表达基因在花粉壁发育相关生物过程显著富集。11个基因参与该过程，其中7个基因为已报道的玉米雄性育性基因，分别为脂肪酰还原酶基因MS6021/ZmMs25/ZmFAR1，ABCG转运蛋白基因MS2/ZmABCG26，细胞色素单加氧酶基因APV1和MS26，PHD-finger型转录因子基因ZmMs7，α-吡喃酮还原酶基因ZmDFR1和ZmDFR2。Zm00001d002342编码多聚半乳糖醛酸酶，Zm00001d047858编码异胡豆苷合酶，Zm00001d019478编码查尔酮合成酶，Zm00001d003702编码4-香豆酰辅酶A连接酶(表6)。

表6 S8b时期花药共同上调差异表达基因中花粉壁发育相关基因

1 348个共同上调差异表达基因在转录调控等分子功能显著富集，其中191个基因被注释为转录因子和转录调节因子，占比14.17%(191/1 348)。MYB，AP2/ERF-ERF和NAC家族转录因子数最多，分别为21，19，17(图4)。4个转录因子基因(ZmMYB84，ZmMs7，MS9，ZmPHD11)为已报道的玉米雄性育性基因。

图4 S8b时期花药共同上调差异表达基因中转录因子分析

3 结论与讨论

花药角质层形成过程包含多类物质代谢途径。目前，通过雄性不育突变体代谢组学和蛋白生化功能分析，部分代谢途径被证明参与玉米花药角质层发育。例如脂肪酸ω羟基化途径，羟基脂肪酸氧化途径，脂酰ACP还原途径等。

本研究通过不同发育时期花药差异表达基因KEGG分析，鉴定到调控玉米花药角质层发育的新代谢途径。Zm00001d012274和Zm00001d042723编码一种乙酰转移酶，可能介导阿魏酰辅酶A和ω羟基棕榈酸酯化途径调控花药角质层发育。水稻乙酰转移酶DPW2催化阿魏酰辅酶A和ω羟基棕榈酸形成阿魏酰棕榈酸酯，参与花药角质形成[30]。环氧脂肪酸是花药角质层角质组成成分，然而其遗传调控机制尚不清楚。Zm00001d025833，Zm00001d017953和Zm00001d051800编码过加氧酶，可能通过脂肪酸环氧化途径影响花药角质层发育。烃类物质是花药角质层蜡质组成成分，水稻OsCER1已被报道通过调节超长链烃形成而参与角质层蜡质合成[31]。Zm00001d014055，Zm00001d051932和Zm00001d017251是水稻OsCER1同源基因，可能通过影响超长链烃合成途径而调控玉米花药角质层发育。

本研究鉴定到11个基因在玉米花粉壁发育过程中发挥重要作用，7个为已报道的雄性育性相关基因。其中MS6021/ZmMs25/ZmFAR1，APV1，MS26，ZmMs7基因突变后均表现出花粉外壁发育缺陷。四聚酮α-吡喃酮还原酶基因ZmDFR1和ZmDFR2单突变体雄性育性下降，其双突变体完全不育[32]。水稻四聚酮α-吡喃酮还原酶基因OsTKPR1突变破坏花粉外壁发育过程[33]。ABCG转运蛋白基因MS2/ZmABCG26突变体完全不育，其水稻同源基因OsABCG15参与花粉外壁发育和雄性育性形成[34]。

多聚半乳糖醛酸酶基因Zm00001d002342催化果胶分子多聚α-(1，4)-聚半乳糖醛酸的裂解。棉花多聚半乳糖醛酸酶GhNSP参与花粉外壁形成过程[35]。异胡豆苷合酶基因Zm00001d047858调控萜类吲哚生物碱的合成。玉米MS45和水稻OsSTRL2为异胡豆苷合酶基因，调控花粉壁形成[36-37]。查尔酮合成酶基因Zm00001d019478以香豆酰辅酶A和丙二酰辅酶A为底物生成查尔酮。拟南芥LAP5和LAP6具有查尔酮合成酶功能，突变体表现为花粉壁和雄性育性异常[38]。4-香豆酰辅酶A连接酶基因Zm00001d003702在香豆酰辅酶A和阿魏酰辅酶A形成中发挥重要功能，可能通过芳香族化合物合成途径调控花粉壁发育。

玉米花药角质层和花粉壁发育受一系列转录因子调控。ZmMs7，ZmbHLH51，ZmbHLH122，ZmMYB84，ZmMYB33-1/-2，ZmPHD11已被证明参与花药角质层发育和雄性育性形成。ZmMs7与核因子Y亚基形成复合体直接调节ZmMT2C表达和绒毡层细胞PCD过程，同时间接调节MS6021和ZmACOS5等基因表达而参与花药角质层和花粉壁发育[20-21]。ZmMYB84通过与ZmMs25启动子结合调控花药角质层和花粉外壁形成[13]。本研究鉴定了191个转录因子和转录调节因子,这些因子可能参与玉米花药角质层和花粉壁发育。后续对这些基因进行功能解析，将完善玉米花药角质层和花粉壁发育的转录调控机制。