李 川，张盼盼，张美微，牛 军，穆蔚林，郭涵潇，乔江方

（河南省农业科学院 粮食作物研究所，河南 郑州 450002）

玉米是我国主要粮食作物，也是饲料和工业原料的重要来源。随着全球气候变暖，极端高温天气频发[1]，严重影响了我国尤其是黄淮海夏玉米主产区玉米的生长发育、产量和品质[2]。一年中的高温时期（7 月下旬—8 月上旬）恰逢河南省夏玉米的散粉吐丝期，该时期玉米植株对温度尤其敏感[3]。花期高温胁迫不仅造成玉米花粉活力降低、雄穗分枝数减少、花药瘦瘪、小花退化，还影响雌穗的分化和受精作用，最终造成籽粒败育，减少产量[4]。除此之外，高温胁迫破坏玉米叶片中叶绿体结构[5]，引起气孔关闭[6]，降低光合作用相关蛋白酶的活性，最终降低光合作用[7]；高温胁迫加速玉米植株生育过程中多种生理生化反应[8]，缩短生育期。灌浆期高温严重缩短玉米灌浆时间[9]，减少籽粒干物质积累量，降低籽粒容重[10]，造成夏玉米产量和品质下降。因此，深入研究玉米高温胁迫的响应机制意义重大，挖掘响应高温胁迫的重要基因及相关调控因子，然后富集分析其本体特征和代谢通路，从而为提高玉米耐高温性奠定理论基础。

MicroRNA（miRNA）是真核生物中一类内源性的具有调控功能的非编码小RNA，长度在18～25 nt，在植物体内参与多种调节途径，包括信号传导过程、生物和非生物胁迫响应过程等[11]。miRNA 在不同植物中具有高度保守性。植物中成熟的miRNA是由长链的初级转录物经过一系列核酸酶的剪切加工而形成的[12-13]。首先，编码miRNA 的基因在细胞核内核苷酸聚合酶的作用下形成长度为几百个核苷酸的miRNA 初始体（pri-miRNA）。然后，Dicer like 1 酶对pri-miRNA 进行切割形成长度为64～303 nt 的miRNA 前体（pre-miRNA）。pre-miRNA 在Dicer like 1酶的作用下形成双链miRNA，即miRNA/miRNA*，然后再被切割成短的双链miRNA/miRNA*。短的双链miRNA/miRNA*在甲基化酶的作用下3'端的最后一个核苷酸发生甲基化修饰。之后甲基化的短双链miRNA/miRNA*由细胞核进入细胞质，完成miRNA 的生物合成过程，形成成熟的miRNA[14-15]。研究发现，植物miRNA 的主要作用原理是通过与AGO（Argonaute）蛋白结合形成RNA 诱导沉默复合体（RISC），指导其配对的下游靶基因mRNA 的切割[16]。植物miRNA 与靶基因mRNA 互补时可以进行特异性切割，miRNA 的剪切位点精确发生在miRNA 配对碱基的第10～11 个核苷酸。这是由于miRNA 的结合位点一般位于基因的开放阅读框中，且大多数植物miRNA 与靶基因mRNA 序列是完全匹配的。当植物miRNA 与其靶基因不能足够互补配对时，miRNA 则主要抑制mRNA 的翻译作用[17-18]。

miRNA 可以调控植物生长发育的多个过程，例如：器官的形态建成、细胞分化和分裂、新陈代谢过程、信号传导、激素分泌等[19]。韦懿等[20]发现，水稻miRNA167 调控水稻植株的分蘖角度和分蘖数目。miRNA 可以通过调控植物激素载体蛋白和受体蛋白的表达来调节植物激素的合成及信号传导途径，最终影响植物的生长发育[21]。miRNA 通过调控相关靶基因的表达水平在植物逆境胁迫应答过程中起到重要调控作用[22-23]，包括生物胁迫和非生物胁迫（干旱胁迫、营养胁迫、低温胁迫、高温胁迫、盐胁迫等）。随着多组学高通量测序技术的发展，miRNA 测序技术在植物研究中被广泛应用。对不同植物种类的营养生长器官进行测序发现，甘蔗中miRNA319 和miRNA396 响 应 冷 胁 迫[24]；小 麦 中miRNA159、miRNA160、miRNA166、miRNA167 响应高 温 胁 迫[25]；miRNA159、miRNA398、miRNA408、miRNA528 在耐旱水稻品种中上调表达[26]；玉米中miRNA1661在汞胁迫条件下下调表达[27]。

目前，关于玉米miRNA 高通量测序研究主要集中于茎叶[28]、根系[29]和雌穗[30]上，主要是针对病原菌侵害[31]、营养胁迫[32]和盐分胁迫[33]。对不同玉米品种高温胁迫处理的花粉进行miRNA 高通量测序的研究还没有报道。为此，以不同耐热性玉米品种生殖器官花粉为试验材料，通过对高温胁迫条件下花粉miRNA 进行高通量测序，筛选高耐热玉米品种郑单958 及低耐热玉米品种先玉335 在高温胁迫处理后花粉中的差异表达miRNA，然后预测其靶基因，并对靶基因进行富集分析，为阐明玉米花粉响应花期高温的分子机制提供参考。

1 材料和方法

1.1 试验地概况及试验材料

试验于2023 年在河南省周口市西华县农业科学研究所试验基地（33°45'44″N、114°26'49″E）进行，土壤类型为黄褐土，肥力均匀，排灌方便。前茬作物为小麦。0～30 cm 耕层土壤基础养分：有机质含量14.10 g/kg、全氮含量0.76 g/kg、碱解氮含量94.97 mg/kg、速效磷含量30.22 mg/kg、速效钾含量169.53 mg/kg，pH值8.79。

供试材料为高耐热玉米品种郑单958（ZD958）、低耐热玉米品种先玉335（XY335）。

1.2 试验设计

ZD958 和XY335 均于2023 年6 月12 日播种，种植密度为67 500 株/hm2，小区30 行，行长5 m，行距0.6 m，重复2次，标准化大田管理。

高温处理采用自建高温棚模拟，用长20 m、宽6 m、高4 m 的高温棚框架固定于田间，周围用透光率为95%的树脂薄膜覆盖，高温棚顶部密封85%，均匀留出15%的空隙便于气体交换。高温棚四周上部、中部悬挂温、湿度计记录高温棚内温、湿度。生长箱框架下部周围树脂薄膜可以拆移，当温度过高时移开薄膜降温，控制高温处理期间高温棚内气温平均高于外界田间3～4 ℃，湿度基本与外界一致。高温处理从抽雄前（2023 年8 月3 日）开始，每天8：30—17：30通过覆盖薄膜进行增温处理。

高温处理7 d后取花粉样本，每个样本由5株盛花期玉米花粉混合组成，生物学重复3次，分别命名为HT958-1、HT958-2、HT958-3、HT335-1、HT335-2、HT335-3。对照为正常生长条件下玉米植株的花粉，分 别 命 名 为CK958-1、CK958-2、CK958-3、CK335-1、CK335-2、CK335-3。花粉样本立即放于干冰中暂时保存，在实验室中于-80 ℃超低温冰箱中冷冻保存。取样结束后拆除高温棚使玉米植株恢复正常生长。

1.3 测定方法

1.3.1 玉米花粉miRNA 高通量测序 将高温处理及正常生长条件下的玉米花粉样本送至杭州联川生物技术股份有限公司进行miRNA 的文库构建和高通量测序。利用TruSeq Small RNA Sample Prep Kit（Illumina，San Diego，USA）试剂盒构建sRNA 文库。之后用Illumina Hiseq 2500进行高通量测序。

1.3.2 玉米花粉miRNA 测序原始数据分析 测序获得的原始序列通过去除3'接头序列和垃圾序列获得干净序列。筛选保留长度在18～25 nt 的序列。将筛选到的序列与mRNA、RFam（非编码RNA 家族数据库）、Repbase（真核生物DNA 中重复序列数据库）数据库进行比对，过滤掉非相关序列和重复序列，从而获得有效序列。对有效序列的长度分布情况进行统计分析。

1.3.3 高温胁迫相关差异表达miRNA 筛选 基于miRBase 数据库和公布的玉米自交系B73 基因组序列 数 据 库（https://www.maizegdb.org/）利 用ACGT101-miR（v4.2）软件进行miRNA 鉴定，并计算每个高温处理和对照样本中miRNA 的表达水平。之后依据不同玉米品种、不同处理中鉴定的miRNA表达水平，筛选差异表达miRNA。差异表达miRNA的筛选标准为|log2FC|＞1 且P＜0.05，其中FC 为差异倍数。

1.3.4 高温胁迫相关差异表达miRNA 的靶基因预测及富集分析 利用GATAr（v1.0）软件对差异表达的miRNA 进行靶基因预测，筛选阈值：不完全匹配位点的最小自由能/完全匹配位点的最小自由能＞0.65，且miRNA 与靶基因匹配位点的罚分值＜4。然后根据GO（Gene ontology）数据库，将差异表达miRNA 候选靶基因按照参与的生物学过程、细胞组分、分子功能进行富集分析。根据KEGG（ Kyoto encyclopedia of genes and genomes）数据库富集分析差异表达miRNA候选靶基因所参与的代谢通路。

1.4 数据处理

采用Excel 2020 和SPSS 25 软件进行试验数据统计分析，用LSD法和Duncan’s 新复极差法进行方差分析和多重比较。

2 结果与分析

2.1 玉米花粉中miRNA测序质量分析

对12 个玉米花粉样本miRNA 高通量测序后获得的原始数据去除接头序列及垃圾序列，获得干净序列，筛选长度在18～25 nt的序列，然后通过数据库比对过滤掉无关和重复序列后获得可信序列。由表1 可见，12 个样本的可信序列Q20 值均在97%之上，可信序列Q30 值均在95%之上，说明测序的碱基错误率很低；GC 含量比率均衡，24 nt序列比值在13.59%～33.93%，可见miRNA 高通量测序结果质量可靠。

表1 12个不同样本中miRNA 测序数据分析Tab.1 Statistical analysis of miRNA sequencing data in the 12 different samples

2.2 不同玉米花粉对比组中差异表达miRNA的筛选

基于miRBase 数据库和玉米自交系B73 基因组序列，在所有样本中共筛选到818 个miRNA 序列。其中，144 个属于玉米中已知miRNA 序列，268 个属于近源物种中已知miRNA 序列，406 个属于全新预测miRNA 序列。在郑单958 高温胁迫花粉与对照花粉对比组（HT958 vs CK958）中，共筛选到19 个差异表达miRNA 序列，其中，15 个miRNA 序列为上调表达，4 个为下调表达。其中，3 个miRNA 序列达到了极显著水平（P＜0.01），分别是（毛果杨）ptcmiR6478_R+2_2ss10GA21GA、（ 葡 萄 ） vvi-MIR3630-p5_2ss19AT20CT、（葡 萄）vvi-MIR3630-p3_2ss19AT20CT；16 个达到显著水平（P＜0.05），分别是（玉米）zma-miR1661-3p_L+2R-2、（玉米）zmamiR159f-5p、（玉米）zma-miR164f-5p、（预测全新）PC-3p-106565_83、（玉米）zma-MIR444a-p5_2、（玉米）zma-MIR444a-p5_1、（玉米）zma-MIR11970-p5_2ss14CT17TC、（水稻）osa-miR162a_R+1、（桃）ppemiR1511-3p_L+1R-4、（亚麻荠）cas-miR5139_L-1、（药 葵）aof-miR5139b_L-3、（水 稻）osa-MIR444ep3、（水 稻）osa-MIR444e-p5、（预 测 全 新）PC-3p-171082_43、（预测全新）PC-5p-221631_29、（高粱）sbi-miR5381_L+1R+1_1ss7CT。

在先玉335 高温胁迫花粉与对照花粉对比组（HT335 vs CK335）中，共筛选到15 个差异表达miRNA 序列，其中，7 个miRNA 序列为上调表达，8个为下调表达。其中，1 个miRNA 序列达到了极显著水平（P＜0.01），为（玉米）zma-miR166j-3p；14个达到显著水平（P＜0.05），分别是（预测全新）PC-3p-30995_367、（ 茶 ） tae-MIR9774-p3_2ss5AG22GT、（玉米）zma-miR166a-3p、（预测全新）PC-3p-54389_191、（ 大 豆）gma-MIR1525-p3_2ss3TA19TG、（水稻）osa-miR535-5p、（二穗短柄草）bdi-MIR444d-p5、（玉米）zma-MIR319a-p5、（玉米）zma-MIR319c-p5、（拟南芥）ath-miR8175_L-2、（预测全新）PC-3p-208306_32、（预测全新）PC-3p-142717_56、（预测全新）PC-5p-17626_673、（预测全新）PC-3p-100896_89。

在郑单958 高温胁迫花粉与先玉335 高温胁迫花粉对比组（HT958 vs HT335）中，共筛选到85 个差异表达miRNA 序列，其中，35 个miRNA 序列为上调表达，50 个为下调表达。其中，24 个miRNA 序列达到了极显著水平（P＜0.01），分别是（预测全新）PC-3p-156933_49、（玉 米）zma-miR166l-3p_L+2R-2、（预测全新）PC-5p-50849_206、（预测全新）PC-3p-134750_61、（玉 米）zma-miR167e-5p、（葡 萄）vvimiR482_1ss5TC、（预 测 全 新）PC-3p-138818_59、（ 预 测 全 新）PC-5p-69263_143、（ 玉 米）zmamiR166j-3p、（预测全新）PC-3p-40151_272、（预测全 新）PC-5p-95787_95、（ 预 测 全 新）PC-5p-106572_83、（玉米）zma-miR167j-3p_L+1、（毛果杨）ptc-miR6478_R+2_1ss21GA、（预 测 全 新）PC-3p-166198_45、（预测全新）PC-5p-118017_73、（预测全新）PC-3p-17225_687、（预测全新）PC-3p-108744_81、（二 穗 短 柄 草）bdi-miR160f_L+2R-1、（预 测 全新）PC-3p-10069_1143、（预测全新）PC-3p-99896_90、（葡萄）vvi-MIR3630-p5_2ss19AT20CT、（葡萄）vvi-MIR3630-p3_2ss19AT20CT、（预 测 全 新）PC-5p-125138_67；4 个miRNA 序列达到极其显著水平（P＜0.001），分别是（预测全新）PC-3p-30995_367、（ 预 测 全 新）PC-5p-17626_673、（ 玉 米）zmamiR166l-3p、（预测全新）PC-5p-16098_734。

在HT958 vs CK958 与HT335 vs CK335 对比组（HT958 vs CK958 VS HT335 vs CK335）中，共筛选到94 个差异表达miRNA 序列，其中，28 个miRNA序列达到了极显著水平（P＜0.01），分别为（预测全新）PC-3p-10069_1143、（预测全新）PC-3p-18335_646、（玉 米）zma-miR164f-5p、PC-3p-16028_737、（玉米）zma-miR159f-5p、（玉米）zma-MIR11970-p3_1ss8AG、ptc-miR6478_R+2_2ss10GA21GA、（高 粱）sbi-miR171a、（玉米）zma-miR171d-3p、（玉米）zma-MIR171b-p3、（玉米）zma-MIR171i-p3、（玉米）zma-MIR171f-p3、（玉米）zma-miR396a-3p、（大豆）gma-MIR1525-p3_2ss3TA19TG、（预 测 全 新）PC-3p-47902_221、（预测全新）PC-3p-96750_94、（二穗短柄草）bdi-miR160f_L+2R-1、（芦笋）aof-miR166b_L-2_1ss20CT、（预测全新）PC-5p-86714_108、（玉米）zma-miR393b-3p_L-1、（毛果杨）ptc-miR6478_R-1_1ss4AG、（ 毛 果 杨 ）ptc-miR6478_R-1_1ss20TA、（玉米）zma-miR171h-3p_L-1R+1、（高粱）sbi-miR5381_L+1R+1_1ss7CT、（ 甜 根 草）ssp-MIR444b-p3、（预测全新）PC-3p-78777_122、（预测全新）PC-3p-30810_369、（玉米）zma-miR169p-5p；17 个miRNA 序列达到极其显著水平（P＜0.001），分别是（预测全新）PC-30995_367、（预测全新）PC-3p-17225_687、（预测全新）PC-3p-7703_1444、（预测全新）PC-5p-16098_734、（预测全新）PC-5p-50849_206、（葡萄）vvi-MIR3630-p5_2ss19AT20CT、（葡萄）vvi-MIR3630-p3_2ss19AT20CT、（二穗短柄草）bdi-MIR444d-p5、（预测全新）PC-5p-90626_102、（预 测 全 新）PC-3p-40151_272、（玉 米）zma-MIR11970-p5_2ss14CT17TC、（预 测 全 新）PC-5p-118815_72、（ 二 穗 短 柄 草）bdi-MIR444d-p5_1ss2GA、（预 测 全 新）PC-3p-100896_89、（预 测 全新）PC-5p-114734_76、（预测全新）PC-3p-53028_196、（水稻）osa-miR444a-3p.2。

2.3 玉米花粉中高温胁迫相关差异表达miRNA靶基因的富集分析

2.3.1 HT958 vs CK958 对比组 利用GSTAr（v1.0）软件对HT958 vs CK958 对比组中检测到的19 个差异表达miRNA 的靶基因进行预测，共获得了503 个基因转录本。将靶基因向GO 数据库中参与的生物学过程、细胞组分、分子功能中不同条目（Term）映射，计算富集P值获得显著富集的GO 条目。由图1可以看出，HT958 vs CK958 对比组中差异表达miRNA 的靶基因富集较多的25 个生物学过程条目分别为转录调控DNA- 模板（Regulation of transcription，DNA-templated）、微 管 生 物 过 程（Microtubule-based process）、磷 酸 化 作 用（Phosphorylation）、RNA 聚合酶Ⅱ正向调控转录过程（Positive regulation of transcription by RNA polymerase Ⅱ）、甲基化作用（Methylation）、RNA 聚合酶Ⅱ调控转录过程（Regulation of transcription by RNA polymerase Ⅱ）、蛋 白 质 转 运（Protein transport）、蛋白质磷酸化（Protein phosphorylation）、有丝分裂DNA 复制纺锤体检查点信号（Mitotic DNA replication checkpoint signaling）、细胞骨架组构（Cytoskeleton organization）、有丝分裂细胞周期（Mitotic cell cycle）、微管细胞骨架组织（Microtubule cytoskeleton organization） 、DNA 复 制（DNA replication）、细胞分裂（Cell division）、信号转导（Signal transduction）、跨 膜 运 输（Transmembrane transport）、分生组织起始（Meristem initiation）、COPⅡ-包膜囊泡出芽（COPⅡ-coated vesicle budding）、韧皮部木质部组织分化（Phloem or xylem histogenesis）、COPⅡ-包膜囊泡货物装载过程（COPⅡ-coated vesicle cargo loading）、DNA 复 制 起 始（DNA replication initiation）、花 发 育（Flower development）、有丝分裂细胞周期相变（Mitotic cell cycle phase transition）、调控细胞周期蛋白依赖丝氨酸/苏氨酸激酶活性（Regulation of cyclin-dependent protein serine/threonine kinase activity）、含核碱基化合物代谢过程（Nucleobase-containing compound metabolic process）。HT958 vs CK958 对比组中差异表达miRNA 的靶基因富集较多的15 个细胞组分条目分别是细胞核（Nucleus）、细胞膜（Membrane）、细胞质（Cytoplasm）、细胞膜的必需组成成分（Integral component of membrane）、原 生 质 膜（Plasma membrane）、微管（Microtubule）、溶质（Cytosol）、细胞骨 架（Cytoskeleton）、内 质 网（Endoplasmic reticulum）、线粒体（Mitochondrion）、高尔基体（Golgi apparatus）、复制识别复合体的核起源（Nuclear origin of replication recognition complex）、内质网出口部位（Endoplasmic reticulum exit site）、COPⅡ囊泡（COPⅡvesicle coat）、内质网向高尔基体转运囊泡膜（ER to Golgi transport vesicle membrane）。HT958 vs CK958 对比组中差异表达miRNA 的靶基因富集较多的10 个分子功能条目分别为蛋白质结合（Protein binding）、DNA 结合（DNA binding）、DNA 结合转录因子活性（DNA-binding transcription factor activity）、核苷酸结合位点（Nucleotide binding）、转移 酶 活 性（Transferase activity）、ATP 结 合（ATP binding）、金属离子结合（Metal ion binding）、RNA 结合（RNA binding）、水解酶活性（Hydrolase activity）、脂质结合（Lipid binding）。

图1 HT958 vs CK958对比组中差异表达miRNA靶基因的GO富集分析Fig.1 GO enrichment analysis of target genes of differentially expressed miRNA in HT958 vs CK958 group

依据KEGG 数据库对HT958 vs CK958 对比组中差异表达miRNA 的靶基因进行代谢通路富集分析（图2）。HT958 vs CK958 对比组中差异表达miRNA 的靶基因KEGG 富集较显著的20 个代谢通路分别是谷胱甘肽代谢（Glutathione metabolism）、碳代谢（Carbon metabolism）、花生四烯酸代谢（Arachidonic acid metabolism）、糖 酵 解/糖 异 生（Glycolysis/gluconeogenesis）、叶酸生物合成（Folate biosynthesis）、丙酸盐代谢（Propanoate metabolism）、赖氨酸降解（Lysine degradation）、β-丙氨酸代谢（Beta-alanine metabolism）、泛酸盐辅酶A 生物合成（Pantothenate and CoA biosynthesis）、色氨酸代谢（Tryptophan metabolism）、缬氨酸亮氨酸及异亮氨酸降解（Valine，leucine and isoleucine degradation）、植物昼夜节奏（Circadian rhythm-plant）、苯丙氨酸酪氨酸及色 氨 酸 生 物 合 成（Phenylalanine，tyrosine and tryptophan biosynthesis）、甘氨酸丝氨酸及苏氨酸代谢（Glycine，serine and threonine metabolism）、戊糖磷酸化通路（Pentose phosphate pathway）、辅酶因子生物合成（Biosynthesis of cofactors）、内质网内蛋白质加工过程（Protein processing in endoplasmic reticulum）、嘧啶代谢（Pyrimidine metabolism）、代谢通路（Metabolic pathways）、乙醛酸及二羧酸代谢（Glyoxylate and dicarboxylate metabolism）。

图2 HT958 vs CK958对比组中差异表达miRNA靶基因的KEGG富集分析Fig.2 KEGG enrichment analysis of target genes of differentially expressed miRNA in HT958 vs CK958 group

2.3.2 HT335 vs CK335 对比组 对HT335 vs CK335 对比组中检测到的15 个差异表达miRNA 的靶基因进行预测，共获得了454 个基因转录本。由图3 可以看出，HT335 vs CK335 对比组中差异表达miRNA 的靶基因富集较多的25 个生物学过程条目分别为转录调控DNA 模板（Regulation of transcription，DNA-templated）、磷 酸 化 作 用（Phosphorylation）、蛋白质磷酸化作用（Protein phosphorylation）、蛋白质水解（Proteolysis）、DNA 修复（DNA repair）、细胞对DNA 损伤刺激的响应（Cellular response to DNA damage stimulus）、端粒维持（Telomere maintenance）、DNA 双链解螺旋（DNA duplex unwinding）、DNA 重组（DNA recombination）、跨 膜 转 运（Transmembrane transport）、防 御 反 应（Defense response）、RNA 聚合酶ⅡC 末端结构域的磷酸化作用（Phosphorylation of RNA polymerase ⅡC-terminal domain）、RNA 聚合酶Ⅱ启动子转录延长的正向调控（Positive regulation of transcription elongation from RNA polymerase Ⅱpromoter）、调控细胞周期（Regulation of cell cycle）、碳水化合物代谢过程（Carbohydrate metabolic process）、植物器官发育（Plant organ development）、代谢过程（Metabolic process）、蛋白质泛素化（Protein ubiquitination）、茉莉酸介导的系统抗性信号通路（Induced systemic resistance， jasmonic acid mediated signaling pathway）、基 于 微 管 的 运 动（Microtubule-based movement）、植物激素响应（Response to auxin）、蛋白质 折 叠（Protein folding）、蛋 白 质 运 输（Protein transport）、蛋白质脱氨（Protein desumoylation）、参与细胞蛋白质分解代谢过程的蛋白质水解（Proteolysis involved in cellular protein catabolic process）。HT335 vs CK335 对比组中差异表达miRNA 靶基因富集较多的15 个细胞组分条目分别是细胞膜（Membrane）、细胞核（Nucleus）、细胞膜必需组成成分（Integral component of membrane）、细胞质（Cytoplasm）、原生质膜（Plasma membrane）、细胞周期蛋白依赖性蛋白激酶全酶复合物（Cyclindependent protein kinase holoenzyme complex）、细胞内膜结合细胞器（Intracellular membrane-bounded organelle）、微管（Microtubule）、胞外区（Extracellular region）、细 胞 质 基 质（Cytosol）、线 粒 体（Mitochondrion）、核糖体（Ribosome）、高尔基体（Golgi apparatus）、溶酶体（Lysosome）、细胞外空间（Extracellular space）。HT335 vs CK335 对比组中差异表达miRNA 靶基因富集较多的10 个分子功能条目分别为ATP 结合（ATP binding）、核苷酸结合（Nucleotide binding）、水 解 酶 活 性（Hydrolase activity）、蛋白质结合（Protein binding）、转移酶活性（Transferase activity）、DNA 结 合 转 录 因 子 活 性（DNA-binding transcription factor activity）、DNA 结合（DNA binding）、激酶活性（Kinase activity）、金属离子结合（Metal ion binding）、蛋白激酶活性（Protein kinase activity）。

图3 HT335 vs CK335对比组中差异表达miRNA靶基因的GO富集分析Fig.3 GO enrichment analysis of target genes of differentially expressed miRNA in HT335 vs CK335 group

由图4 可以看出，HT335 vs CK335 对比组中差异表达miRNA 的靶基因富集较显著的20 个代谢通路分别是其他多糖降解（Other glycan degradation）、亚油酸代谢（Linoleic acid metabolism）、代谢通路（Metabolic pathways）、硫 胺 素 代 谢（Thiamine metabolism）、内质网内蛋白质加工过程（Protein processing in endoplasmic reticulum）、脂肪酸延长（Fatty acid elongation）、植物激素信号传导（Plant hormone signal transduction）、类固醇生物合成（Steroid biosynthesis）、精氨酸生物合成（Arginine biosynthesis）、不同类型的N- 聚糖生物合成（Various types of N-glycan biosynthesis）、酪氨酸代谢（Tyrosine metabolism）、鞘脂类代谢（Sphingolipid metabolism）、脂肪酸降解（Fatty acid degradation）、N-聚糖生物合成（N-glycan biosynthesis）、α-亚麻酸代谢（Alpha-linolenic acid metabolism）、缬氨酸亮氨酸及异亮氨酸降解（Valine，leucine and isoleucine degradation） 、类 苯 基 丙 烷 生 物 合 成（Phenylpropanoid biosynthesis） 、核 质 转 运（Nucleocytoplasmic transport） 、 内 吞 作 用（Endocytosis）、丙氨酸天冬氨酸及谷氨酸代谢（Alanine，aspartate and glutamate metabolism）。

图4 HT335 vs CK335对比组中差异表达miRNA靶基因的KEGG富集分析Fig.4 KEGG enrichment analysis of target genes of differentially expressed miRNA in HT335 vs CK335 group

2.3.3 HT958 vs HT335 对比组 对HT958 vs HT335 对比组中检测到的85 个差异表达miRNA 的靶基因进行预测，共获得了2 286 个基因转录本。由图5 可以看出，HT958 vs HT335 对比组中差异表达miRNA 的靶基因富集较多的25 个生物学过程条目分别为转录调控DNA- 模板（Regulation of transcription，DNA-templated）、磷 酸 化 作 用（Phosphorylation）、蛋 白 质 磷 酸 化（Protein phosphorylation）、蛋白质水解（Proteolysis）、跨膜转运（Transmembrane transport）、碳水化合物代谢过程（Carbohydrate metabolic process）、代 谢 过 程（Metabolic process）、蛋 白 质 泛 素 化（Protein ubiquitination）、防御反应（Defense response）、RNA聚合酶Ⅱ的转录调控（Regulation of transcription by RNA polymerase Ⅱ）、DNA 修复（DNA repair）、甲基化作用（Methylation）、细胞对DNA 损伤刺激的响应（Cellular response to DNA damage stimulus）、初级代谢过程（Primary metabolic process）、蛋白质转运（Protein transport）、多聚糖分解代谢过程（Polysaccharide catabolic process）、DNA 重组（DNA recombination）、细胞大分子代谢过程（Cellular macromolecule metabolic process）、植物激素激活的信号通路（Auxin-activated signaling pathway）、DNA双链解螺旋（DNA duplex unwinding）、信号传导（Signal transduction）、脂 质 代 谢 过 程（Lipid metabolic process）、泛素依赖性蛋白质分解代谢过程（Ubiquitin-dependent protein catabolic process）、端粒维持（Telomere maintenance）、mRNA 加工（mRNA processing）。HT958 vs HT335 对比组中差异表达miRNA 的靶基因富集较多的15 个细胞组分条目分别是细胞膜（Membrane）、细胞膜必需组成成分（Integral component of membrane）、细 胞 核（Nucleus）、细胞质（Cytoplasm）、原生质膜（Plasma membrane）、细 胞 溶 质（Cytosol）、胞 外 区（Extracellular region）、高尔基体（Golgi apparatus）、叶绿体（Chloroplast）、色素体（plastid）、内质网膜（Endoplasmic reticulum）、线粒体（Mitochondrion）、胞间连丝（Plasmodesma）、细胞内膜结合细胞器（Intracellular membrane-bounded organelle）、高尔基体膜（Golgi membrane）。HT958 vs HT335 对比组中差异表达miRNA 靶基因富集较多的10 个分子功能条目分别为转移酶活性（Transferase activity）、ATP结合（ATP binding）、DNA 结合（DNA binding）、核苷酸结合（Nucleotide binding）、蛋白质结合（Protein binding）、水解酶活性（Hydrolase activity）、金属离子结合（Metal ion binding）、激酶活性（Kinase activity）、蛋白激酶活性（Protein kinase activity）、DNA 结合转录因子活性（DNA-binding transcription factor activity）。

图5 HT958 vs HT335对比组中差异表达miRNA靶基因的GO富集分析Fig.5 GO enrichment analysis of target genes of differentially expressed miRNA in HT958 vs HT335 group

HT958 vs HT335 对比组中差异表达miRNA 的靶基因富集较显著的20 个代谢通路分别是鞘脂类代谢（Sphingolipid metabolism）、淀粉和蔗糖代谢（Starch and sucrose metabolism）、其他多糖降解（Other glycan degradation）、代谢通路（Metabolic pathways）、半胱氨酸及甲硫氨酸代谢（Cysteine and methionine metabolism）、卟 啉 代 谢（Porphyrin metabolism）、类胡萝卜素生物合成（Carotenoid biosynthesis）、类 单 萜 生 物 合 成（Monoterpenoid biosynthesis）、精 氨 酸 生 物 合 成（Arginine biosynthesis）、氨 基 酸 生 物 合 成（Biosynthesis of amino acids）、叶 酸 一 碳 库（One carbon pool by folate）、维生素B6 代谢（Vitamin B6 metabolism）、硫依赖系统（Sulfur relay system）、植物昼夜节奏（Circadian rhythm-plant） 、 核 质 转 运（Nucleocytoplasmic transport）、次级代谢物生物合成（Biosynthesis of secondary metabolites）、植物激素信号传导（Plant hormone signal transduction）、ABC 转运体（ABC transporters）、戊糖磷酸化通路（Pentose phosphate pathway）、吞噬体（Phagosome）（图6）。

图6 HT958 vs HT335对比组中差异表达miRNA靶基因的KEGG富集分析Fig 6 KEGG enrichment analysis of target genes of differentially expressed miRNA in HT335 vs CK335 group

2.3.4 HT958 vs CK958 VS HT335 vs CK335 对比组 对HT958 vs CK958 VS HT335 vs CK335 对比组中检测到的94 个差异表达miRNA 的靶基因进行预测，共获得了4 569 个基因转录本。由图7 可以看出，HT958 vs CK958 VS HT335 vs CK335 对比组中差异表达miRNA 的靶基因富集较多的25 个生物学过程条目分别为转录调控DNA-模板（Regulation of transcription，DNA-templated）、磷 酸 化 作 用（Phosphorylation）、蛋 白 质 磷 酸 化（Protein phosphorylation）、蛋白质转运（Protein transport）、蛋白质水解（Proteolysis）、RNA 聚合酶Ⅱ的转录调控（Regulation of transcription by RNA polymerase Ⅱ）、蛋白质泛素化（Protein ubiquitination）、跨膜转运（Transmembrane transport）、初级代谢过程（Primary metabolic process）、甲基化作用（Methylation）、细胞大分子代谢过程（Cellular macromolecule metabolic process）、翻 译（Translation）、DNA 修 复（DNA repair）、细胞内蛋白质转运（Intracellular protein transport）、细胞对DNA 损伤刺激的响应（Cellular response to DNA damage stimulus）、泛素依赖的蛋白质分解代谢过程（Ubiquitin-dependent protein catabolic process）、囊 泡 介 导 的 转 运（Vesiclemediated transport）、调控催化剂活性（Regulation of catalytic activity）、信号传导（Signal transduction）、离子 传 输（Ion transport）、防 御 响 应（Defense response）、碳水化合物代谢过程（Carbohydrate metabolic process）、细胞分裂（Cell division）、DNA 重组（DNA recombination）、转录正向调控DNA-模板（Positive regulation of transcription，DNA-templated）。

图7 HT958 vs CK958 VS HT335 vs CK335对比组中差异表达miRNA靶基因的GO富集分析Fig.7 GO enrichment analysis of target genes of differentially expressed miRNA in HT958 vs CK958 VS HT335 vs CK335 group

HT958 vs CK958 VS HT335 vs CK335 对比组中差异表达miRNA 的靶基因富集较多的15 个细胞组分条目分别是细胞核（Nucleus）、细胞膜（Membrane）、细胞膜必需组成成分（Integral component of membrane）、细胞质（Cytoplasm）、原生质膜（Plasma membrane）、细 胞 溶 质（Cytosol）、内 质 网（Endoplasmic reticulum）、高 尔 基 体（Golgi apparatus）、线粒体（Mitochondrion）、细胞内膜结合细胞器（Intracellular membrane-bounded organelle）、核糖体（Ribosome）、细胞核（nucleolus）、叶绿体（Chloroplast）、色素体（Plastid）、微管（Microtubule）。HT958 vs CK958 VS HT335 vs CK335 对比组中差异表达miRNA 靶基因富集较多的10 个分子功能条目分别为蛋白质结合（Protein binding）、DNA 结合（DNA binding）、ATP 结合（ATP binding）、核苷酸结合（Nucleotide binding）、转移酶活性（Transferase activity）、金属离子结合（Metal ion binding）、DNA 结合转录因子活性（DNA-binding transcription factor activity）、水解酶活性（Hydrolase activity）、激酶活性（Kinase activity）、蛋 白 激 酶 活 性（Protein kinase activity）。

由图8 可见，HT958 vs CK958 VS HT335 vs CK335 对比组中差异表达miRNA 的靶基因富集较显著的20 个代谢通路分别是内质网内蛋白质加工过程（Protein processing in endoplasmic reticulum）、真核生物核糖体生物合成（Ribosome biogenesis in eukaryotes）、剪接体（Spliceosome）、鞘脂类代谢（Sphingolipid metabolism）、内吞作用（Endocytosis）、其他多糖降解（Other glycan degradation）、吞噬体（Phagosome）、氨酰转运RNA 生物合成（AminoacyltRNA biosynthesis）、核 质 转 运（Nucleocytoplasmic transport）、ABC转运体（ABC transporters）、谷胱甘肽代谢（Glutathione metabolism）、泛酸盐及辅酶A生物合成（Pantothenate and CoA biosynthesis）、植物激素信号传导（Plant hormone signal transduction）、光合作用天线蛋 白（Photosynthesis-antenna proteins）、DNA 修复（DNA replication）、蛋白质输出（Protein export）、赖氨酸降解（Lysine degradation）、泛醌和其他萜类醌生物合成（Ubiquinone and other terpenoidquinone biosynthesis）、叶酸一碳库（One carbon pool by folate）、单 胺 菌 素 生 物（Monobactam biosynthesis）。

图8 HT958 vs CK958 VS HT335 vs CK335对比组中差异表达miRNA靶基因的KEGG富集分析Fig.8 KEGG enrichment analysis of target genes of differentially expressed miRNA in HT958 vs CK958 VS HT335 vs CK335 group

3 结论与讨论

本研究发现，在高温胁迫下，郑单958 花粉中miRNA6478 极显著差异表达。在NaCl 胁迫下，毛竹中miRNA6478 也上调表达，调控毛竹萌发[34]。甘蔗中miRNA6478 通过调控乙烯信号转导途径基因CTR1（Constitutive triple response 1）的表达来负调控乙烯信号转导响应黑穗病菌胁迫[35]。另外，miRNA3630 在高温胁迫下郑单958 花粉中极显著上调表达。研究发现，小桐子中miRNA3630响应低温胁迫显著下调表达[36]。泡桐中miRNA3630 通过调控靶基因NLR（Nucleotide biding site-leucine rich repeat）基因家族成员在抗病中发挥作用，NLR 基因家族成员启动子区含有水杨酸、茉莉酸甲酯等激素响应元件[37]。本研究发现，在高温胁迫下，低耐热玉米品种先玉335 中miRNA166 极显著下调表达。研究发现，在低温胁迫下，miRNA166 在冷敏感玉米自交系B73 叶片中上调表达[38]。42 ℃高温胁迫处理15 d 后，miRNA166 在耐热性强蕹菜品种中下调表达，在耐热性较差蕹菜品种中显著上调表达[39]。

本研究发现，在高温胁迫下，不同耐热玉米品种花粉中miRNA1661 显著下调表达。对水稻miRNA 基因上游1 000 bp 序列进行抗性元件预测发现，miRNA1661 含有GA 元件，其对光反应敏感[40]。低磷胁迫处理玉米根系24 h 后发现，miRNA1661 下调表达[41]。miRNA167 在禾本科植物中高度保守，响应非生物胁迫。拟南芥中miRNA167 调控花器官发育，其靶基因为生长素响应因子，生长素响应因子参与植物激素信号传导，调节多种植物的生长发育过程[42]。渗透胁迫条件下，拟南芥miRNA167下调表达，促进生长素积累及侧根伸长[43]。在番茄中过量表达miRNA167 造成叶片变小、雄蕊及花柱变短，引发雄性不育[44]。本研究发现，高温胁迫下，高耐热玉米品种郑单958花粉中miRNA167、miRNA482 均显著差异表达。在4 ℃储藏条件下，马铃薯中miRNA482上调表达，参与低温糖化调控[45]；在敲除miRNA482 的转基因番茄中，抗氧化酶活性升高，脯氨酸含量增加，丙二醛含量降低，最终提高耐盐性[45]。miRNA482 的靶基因为响应胁迫相关的CBF（C-repeat responsive element binding factor）转录因子家族成员和由GAI（Gibberellic acid insensitive）、RGA（Repressor of GA1-3 mutant）和SCR（Scarecrow）所组成的GRAS超级转录因子家族成员。miRNA482 在盐胁迫大豆根尖中上调表达，推测miRNA482 参与调控大豆对盐胁迫的响应[46]。本研究发现，miRNA160 在高温胁迫处理郑单958 花粉中极显著表达。miRNA160 在高等植物中普遍高度保守，因此，调控网络复杂且多样化。miRNA160 通过调控ARF（Auxin response factors）家族成员参与植物生长发育及响应生物胁迫过程。过表达miRNA160 的拟南芥开花时间提前，miRNA160 参与了花期调控[47]。预测发现，小麦中miRNA160 的靶基因为过氧化物酶基因、ARF 家族基因、钾离子转运蛋白基因等[48]。低温胁迫条件下，西瓜中miRNA160 下调表达，负调控靶基因表达，从而提高西瓜耐冷性[49]。

本研究发现，郑单958 玉米花粉中（HT958 vs CK958对比组）19个差异表达miRNA的靶基因显著富集到GO:0007165 信号传导（Signal transduction）条目和GO:0009908 花发育（Flower development）条目。 GO:0007165 信号传导（Signal transduction）条目相关的miRNA 包括miRNA1661，其靶基因为Zm00001eb223610；miRNA11970，其 靶 基 因 为Zm00001eb206990；miRNA1511，其 靶 基 因 为Zm00001eb182010；miRNA5139，其 靶 基 因 为Zm00001eb182010。由此可以推测miRNA1661、miRNA11970、miRNA1511、miRNA5139 通过信号传导模式共同调控玉米花粉的耐热性。GO:0009908花发育（Flower development）条目相关的miRNA 为miRNA164，其靶基因为Zm00001eb264380。推测miRNA164 通过正调控靶基因Zm00001eb264380 的表达来调控玉米花的发育，从而提高郑单958 花粉的耐热性。

本研究发现，先玉335 玉米花粉（HT335 vs CK335 对比组）中15 个高温胁迫差异表达miRNA的靶基因显著富集到GO：0009864 茉莉酸介导的系统抗性信号通路（Induced systemic resistance，jasmonic acid mediated signaling pathway） 、GO：0009733 植物激素响应（Response to auxin）、GO：0004672 蛋白激酶活性（Protein kinase activity）、GO：0099402 植 物 器 官 发 育（Plant organ development）、GO：0016301 激 酶 活 性（Kinase activity）。GO:0009864 茉莉酸介导的系统抗性信号通 路（Induced systemic resistance，jasmonic acid mediated signaling pathway）条 目 相 关 的miRNA 为miRNA8175，其 靶 基 因 为 Zm00001eb363100、Zm00001eb196470、Zm00001eb093200、Zm00001eb-141160、Zm00001eb405060。GO：0009733 植物激素响 应（Response to auxin）条 目 相 关 的miRNA 为miRNA8175，其 靶 基 因 为 Zm00001eb343460、Zm00001eb353400、Zm00001eb353410、Zm00001eb-353370。GO:0016301 激酶活性（Kinase activity）条目相关的miRNA 为miRNA8175、miRNA9774。由此推测miRNA8175 通过调控茉莉酸介导的系统抗性信号通路、植物激素响应等条目相关的多个靶基因的表达正调控先玉335花粉耐热性。

本研究首先通过miRNA 高通量测序技术对高温胁迫条件下不同耐热玉米品种花粉中的miRNA进行了对比分析，筛选出高温胁迫相关的显著差异表达miRNA，并对差异表达miRNA 的靶基因进行了GO 功能注释和KEGG 通路分析。高温胁迫相关miRNA 及其靶基因通过复杂调控网络在玉米花期高温胁迫响应中发挥着重要作用，为进一步探明具体miRNA功能和阐明玉米耐热性提供理论依据。