玉米是我国重要的粮食作物之一，具有高产、适应性广和营养物质丰富等特点，因此玉米高产稳产在农业生产中占有重要地位^[1]。玉米生长发育过程中经常发生倒伏现象，造成产量大幅下降，据统计^[2-3]，倒伏可造成玉米减产15%~20%，导致我国每年玉米产量损失约100万t。玉米倒伏主要包括根倒、茎倒和茎折断3种形式。导致玉米倒伏的因素有气候因素、种植密度、土壤质地以及玉米株高、茎秆柔韧性和根系发达程度等自身遗传因素^[4⇓-6]。

茎秆力学特征是影响玉米茎秆倒伏的重要因素之一，代表玉米抗倒伏能力的茎秆力学特征有茎秆拉力、抗推力、穿刺强度和压碎强度，且与玉米抗倒伏能力呈显著正相关^[7]。茎秆显微结构是影响玉米茎秆抗倒伏性能的另一个重要因素，维管束数目越多、细胞表皮越厚的玉米茎秆强度越好，抗倒性越强。冯素伟等^[8]对茎秆显微结构进行分析发现，茎秆显微结构中大维管束数量越多，玉米的抗倒能力越强。前人对玉米倒伏的研究主要集中在茎秆力学特征、形态学特征以及相应的生理生化指标方面，在分子水平对玉米抗倒伏相关基因的挖掘的研究较少。因此，本研究以高抗倒的京农科728、中抗倒的金农738和低抗倒的先玉335为材料，在玉米抽雄期取第3节间茎秆进行转录学测序，挖掘与玉米抗倒性相关的基因，为培育高产和高抗倒伏的玉米品种提供理论依据。

试验于2019年在河北省张家口市沙岭子镇农业科学院（115°05′E，40°6′N）进行，土壤质地为壤质土。供试玉米品种为京农科728（JNK728，高抗倒性，H）、金农738（JN738，中抗倒性，M）和先玉335（XY335，低抗倒性，L），试验采取完全随机区组设计，3次重复，行距55cm，株距30cm，种植密度67 500株/hm²，每个品种12行，行长10m。在玉米抽雄期，每个小区取3株长势一致的玉米，取第3节间茎秆液氮速冻后于-80℃保存待测。

使用康为全能型植物RNA提取试剂盒（OminiPlant RNA Kit）提取9个玉米茎秆样本的总RNA。使用Qubit 2.0荧光计测定RNA的浓度，用Agilent 2100生物分析仪和1%的琼脂糖凝胶电泳检测RNA的完整性。使用Illumina公司的TruSeqTM RNA试剂盒（San Diego）构建测序所需cDNA文库，用Agilent 2100生物分析仪检测文库质量。由武汉迈特维尔生物科技公司完成建库测序。

使用Illumina HiSeq高通量测序平台对cDNA文库进行测序，对生成的原始数据进行处理，除去低质量读段、含接头或引物的读段，利用HISAT2将高质量读段与参考基因组（B73 RefGen_v4）进行序列比对，获取在参考基因组或基因上的位置信息以及测序样品特有的序列特征信息。所有下游分析均基于与参考数据完美匹配或只有一处不匹配的读段。采用FPKM（fragments per kilobase per million mapped reads）计算转录本或基因的表达水平。

为了鉴定与玉米茎折相关基因，用DESeq2分别对3个品种组间的差异表达进行分析，基因表达差异倍数在2倍以上、且FDR<0.05确定为差异表达基因（differentially expressed genes，DEGs）。为了进一步鉴定DEGs的功能，使用Blast2GO（https://www.blast2go.com/）和KOBAS 2.1.1 (http://kobas.cbi.pku.edu.cn/download.php)进行GO和代谢通路富集分析（FDR<0.05）；使用Venny2.1.0构建Venn图。

为了验证RNA测序结果的准确性，选择10个DEGs进行qRT-PCR验证。采用Primer 5对DEGs进行荧光定量引物设计。内参基因为玉米GAPDH，引物信息见表1。使用Bio-Rad iQ5荧光定量PCR仪（Bio-RAD，美国），选用2×Fast Super EvaGreen ® qPCR Mastermix（EverbrightInc.，美国）进行荧光定量分析。反应体系为cDNA模板1µL、前后引物各1µL、ddH₂O 7µL、2×AugeGreenTMMaster Mix 10µL。用2^-ΔΔCT法^[9]计算基因的相对表达量。

用HiSeq6000平台对9个玉米茎秆样本进行双端测序后共获得59.26Gb高质量读段，与玉米参考基因组比对后发现，84.84%~87.28%的高质量读段比对到参考基因的唯一位置，3.04%~3.58%的有效数据比对到参考基因的多个位置，Q30和GC含量分别高于93.90%和53.16%（表2）。其次，用主成分分析（PCA）评估样本之间相关性，结果（图1）表明，生物学重复之间高度相关，每个材料间明确分离。以上结果证实转录组测序数据是可靠的，可用于进一步分析。

为了挖掘玉米抗倒伏相关基因，对3个抗倒性不同的玉米材料进行品种间基因差异表达分析。由表3可知，先玉335与京农科728比较（L-vs-H）鉴定到差异基因最多，其中4411个DEGs上调表达，3368个DEGs下调表达。先玉335与京农738比较（L-vs-M）鉴定到5373个DEGs，包括2896个DEGs下调表达，2477个DEGs上调表达。京农738与京农科728比较（M-vs-H）鉴定到4905个DEGs，2831个DEGs下调表达，2074个DEGs上调表达。3个比较组共鉴定到10 093个DEGs，其中2095、769、606个DEGs分别在L-vs-H、L-vs-M、M-vs-H分组特异表达，1341个DEGs在3个分组中共同表达（图2）。

为了探究抗倒性不同玉米品种差异表达基因的生物学功能，对鉴定到的DEGs进行GO功能富集分析。由图3可知，生物学过程分类中，L-vs-H分组显著富集到植物型细胞壁的生物发生、光合作用、木质素合成过程，L-vs-M分组显著富集到植物型次生细胞壁的生物发生，M-vs-H分组显著富集到次生代谢过程、细胞壁发生、苯丙素的代谢过程。分子功能分类中，L-vs-H、M-vs-H分组显著富集到血红素结合、单加氧酶活性、碳水化合物结合；L-vs-M分组显著富集到叶绿素结合、四吡咯结合、糖基转移酶活性。细胞成分分类中，光合体系、光合体系Ⅰ、光合体系Ⅱ在L-vs-H、L-vs-M分组中显著富集，非原质体在M-vs-H分组中显著富集。因此，3个玉米品种抗倒性不同可能是因为富集到相同GO条目的DEGs数量不同。

为了近一步了解DEGs的功能，对3个分组鉴定到的DEGs进行KEGG代谢通路富集分析。结果（图4）表明，L-vs-H和M-vs-H分组均富集到6条代谢通路，L-vs-M分组富集到5条代谢通路。其中，苯丙素生物合成、次生代谢产物的生物合成、类黄酮生物合成在L-vs-H和M-vs-H分组显著富集。L-vs-M分组中，光合作用―天线蛋白、植物―病原体相互作用、次生代谢产物的生物合成显著富集。

3个玉米品种的显微镜可视范围内茎秆维管束数目、单个维管束面积、茎秆表皮细胞厚度如图5所示。先玉335、金农738和京农科728显微镜可视范围内茎秆维管束数目分别为32、39、44个（图5a和图6a）；先玉335单个维管束面积最小，京农科728单个维管束面积最大（图5b和图6b）；茎秆表皮细胞厚度的顺序依次为京农科728>金农738>先玉335（图5c和图6c）。由此可知，京农科728的高抗倒性可能与较少的维管束数目、较大的单个维管束面积以及较厚的茎秆表皮细胞厚度有关。

为了验证转录组测序的准确性，通过qRT-PCR检测了10个已报道与茎秆抗倒相关或注释到与茎秆抗倒相关途径的DEGs表达量，结果（图7）表明，所有10个基因的表达水平与RNA-seq测序的结果一致，且相关系数达到0.9171，因此qRT-PCR的分析结果证实了转录组测序结果的可靠性。

玉米在我国粮食安全战略中占有重要地位，随着玉米种植密度增加，茎秆增高，茎秆倒伏问题日益严重，成为限制玉米产量和品质的因素之一。作物倒伏指由各种外界因素引发的植株茎秆由自然直立状态到永久错位的现象，是农作物生产中存在的一个普遍性问题。导致作物倒伏的原因主要有遗传、栽培管理和自然环境3个方面，其中由作物基因型决定的遗传因素是内因，也是决定作物自身是否抗倒的最直接、最根本的因素^[10]。植物生长处于风、雨等自然条件的影响下，当作物的茎秆弹性不足以使弯曲的植株恢复时容易发生大面积倒伏。玉米茎秆单位面积内维管束数目和表皮细胞厚度决定了玉米茎秆的抗压强度，本研究通过比较分析不同玉米品种的维管束数目和表皮细胞厚度，结果表明维管束数目越多，表皮细胞越厚，玉米茎秆抗折强度越高。杨硕等^[11]对玉米茎秆显微结构研究表明，茎秆维管束数目、大维管束数目和小维管束数目在不同环境下均与茎秆抗倒伏能力呈正相关。同样，对水稻的研究^[12]发现，抗倒伏能力越强的水稻，其茎壁厚，大、小维管束数目多。

茎秆抗倒伏是由多个基因控制的复杂的数量性状，目前对于玉米茎秆抗倒伏相关基因与代谢通路的研究较少。因此，本研究对3个抗倒性不同玉米品种进行转录学测序，并对DEGs进行了GO和KEGG富集分析，报道了玉米抗倒伏的关键差异表达基因和代谢通路。木质素是一类重要的大分子有机物质，由类苯丙酸途径的单体衍生形成，在维持细胞壁结构完整性方面起着关键性作用，茎秆的木质素含量越高作物的抗倒性越强^[13]。木质素合成途径中，CoA连接酶（4CL）引导光合代谢产物向木质素合成途径的流动，处于代谢分支点的位置^[14]。研究^[10]表明，当4CL活性比原来提升1倍时，木质素含量也会相应提高1/4。本研究中苯丙素生物合成、木聚糖生物合成过程、细胞壁生物发生和植物型次生细胞壁的生物发生在3个分组中显著富集，对大豆茎倒伏的研究^[10]表明DEGs也与这些功能相关。已报道的与茎倒伏相关的苯丙烷类代谢途径中的关键酶苯丙氨酸解氨酶、肉桂酸4-羟化酶和连接植物苯丙烷复合途径和木质素特异生物合成途径的4CL显著差异表达。推测这些显著富集的通路和差异基因可能与玉米抗倒性相关，其功能还需进一步验证。

对3个抗倒性不同的玉米品种抽雄期茎秆进行转录组学分析，3个比较组共鉴定到10 093个DEGs，其中2095、769、606个DEGs分别在L-vs-H、L-vs-M和M-vs-H分组特异表达，1341个DEGs在3个分组中共同表达。GO富集分析表明，3个玉米品种抗倒性不同可能是因为富集到相同GO条目的DEGs数量不同。代谢通路富集分析表明，L-vs-H、M-vs-H分组显著富集到苯丙素生物合成、次生代谢产物的生物合成、类黄酮生物合成途径，而L-vs-M分组显著富集到光合作用―天线蛋白、植物―病原体相互作用途径，注释到这些途径的基因为与玉米抗倒伏相关的关键候选基因。茎秆的显微结构分析表明，京农科728的抗倒性强可能与茎秆维管束数目多、单个维管束面积大、茎秆表皮细胞厚有关。此外，qRT-PCR验证了转录组测序的结果。本研究确定的关键基因和代谢途径可作为未来定向克隆和下游分析研究的宝贵遗传资源或选择目标。