转录组测序分析甘蓝型油菜DW871矮化性状

罗京，李超*，张瑞茂，赵德刚，高志宏，王芳2,，杨元雨2,，王转转2,，王敏

（1.贵州师范大学生命科学学院，贵州 贵阳，550025；2.贵州省农业科学院油料研究所，贵州 贵阳，550006；3.贵州省农业科学院油菜研究所，贵州 贵阳，550009；4.贵州大学生命科学学院/农业生物工程研究院山地植物资源保护与种质创新省部共建教育部重点实验室，山地 生态与农业生物工程协同创新中心，贵州 贵阳，550025）

油菜（Brassica napusL.）是重要的油料经济作物，也是中国第一大油料作物[1]。目前中国油菜的主要品种包括甘蓝型油菜、芥菜型油菜和白菜型油菜。其中甘蓝型油菜占我国油菜全部种植面积的90%[2]。此外，甘蓝型油菜具有高抗病性、高含油量等特点，被誉为最成功的杂交作物之一[3]。但甘蓝型油菜普遍为高秆、易倒伏，油菜倒伏后产量减少幅度为10%～30%，且不利于机械化生产，严重限制了油菜的推广[4]。培育具有优良性状的矮化品种成为解决这一问题的主要途径。

继小麦、水稻之后，油菜成为又一主导绿色革命的作物。研究发现适当矮化的甘蓝型油菜品种具有较强的抗倒伏性，更易于进行机械化收获，并具有耐肥、抗倒伏、耐密植、经济系数高等特点，从而为大规模机械化生产创造条件[5,6]。迄今，围绕油菜矮化育种，在矮化材料选育、基因定位、遗传规律、矮化机制等方面开展了大量研究[7～13]。

植物的株高受多种因素共同控制，除内部基因调控外，也受外部环境因素、各种激素和逆境胁迫效应的影响。其中，激素在植物株型的建成过程中发挥了重要的作用，且不同激素之间存在一定的互作关系。如Müller等构建的关于生长素和细胞分裂素的模型中，生长素通过调控A 型ARR 的表达介导细胞分裂素的输出[14]。安文燕等发现陆地棉矮化突变体Ari1327因DELLA 蛋白累积，使得下游基因对GA 的响应被抑制，从而出现了矮化，当施加外援激素时，植株恢复为野生型[15]。Zhao 等发现甘蓝型油菜矮化突变体ds-4中生长素受体蛋白基因BnaC05. IAA7的突变，是导致矮化的原因[16]。牛显飞等发现甘蓝型油菜矮化突变体NDF-1的矮化与赤霉素、赤霉素和油菜素内酯有关，且细胞的伸长和细胞壁的发育都存在异常[17]。此外，张玲等发现水稻矮化突变体ddf2中，维管束中细胞的数量和大小显著小于野生型，突变体茎细胞分裂和膨胀都受到了抑制，与植株的矮化密切相关[18]。Wang 等发现多效性基因座BnDWF/DCL1的突变会导致株高的降低，叶片形状和株型也会随之改变[19]。

DW871 作为优秀的甘蓝型油菜矮秆材料，与野生型高秆相比株型紧凑，木质部较厚且木质化进程较快，收获指数达到0.38，与同源高秆相当。本文通过UMI（unique molecular identifier）-RNA-seq 对甘蓝型油菜矮化突变体DW871进行转录组分析，并以甘蓝型油菜高秆材料HW871为对照，寻找在矮化性状形成过程中起关键作用的候选基因和生物学途径，探讨DW871 矮化性状形成的机制，为后续的育种工作和矮化基因的定位和应用提供依据。

1 材料与方法

1.1 材料

甘蓝型矮秆油菜DW871 和对照组甘蓝型高秆油菜HW871 由贵州省农业科学院油菜研究所提供[20]。DW871 为隐性核不育三系杂交组合ZH 117后代中的矮秆紧凑株型突变体，HW871是从DW871杂合型中分离出的正常高秆油菜[21]。种植和田间管理保持一致。

选取处于抽薹期的甘蓝型油菜用于转录组测序，为避免因选材部位的差异导致测序结果的误差，统一选取倒数第四茎节，设置4 次生物学重复，共8 个样本，以确保实验的准确性和可靠性。将茎节分别置于液氮中速冻，-80℃保存。

1.2 方法

1.2.1 总RNA 的提取和文库构建 本实验使用TRIzol 试剂法提取总RNA，获得的总RNA 经过NanoDrop 2000 和Qubit 核酸蛋白定量仪检测总RNA 质量，并经过琼脂糖凝胶电泳检测合格。使用Plant RNA Purification Reagent 纯化后，对合格的RNA 样品进行cDNA 文库的构建，最终在Illumina Novaseq测序平台上对文库进行测序。

1.2.2 测序结果的统计和功能注释 Illumina 测序结果经trimmomatic 软件对原始测序数据进行处理，获得高质量的clean reads 后，在clean reads 上连接UMI（unique molecular identifier）标签，排除PCR扩增偏好性和测序偏好性，获得UMI reads。使用Hisat2 软件将之比对到甘蓝型油菜参考基因组上。以FPKM（Reads/Fragments per kilobase of exon model per million mapped reads）值来计算基因的表达水平，并计算TPM（transcripts per kilobase of exon model per million mapped reads）值对结果进行处理，以消除exon 长度和样本间测序总read count 不同造成的差异。应用edgeR软件来分析基因的表达差异。

1.2.3 数据的统计与分析 选择与矮化性状可能存在关联的代谢通路进行分析。使用软件Goatools进行GO 富集分析，使用KOBAS 进行KEGG PATHWAY富集分析，筛选出与矮化性状形成相关的差异表达基因（DEG），从基因的表达水平和蛋白互作等方面进行分析。在FDR≤0.05 和注释到该通路基因大于等于50个条件的基础上，筛选与油菜生长发育密切相关的代谢通路[22,23]。

1.2.4 qRT-PCR 验证 利用qPrimerDB 根据geneID 设计引物[24]。将提取RNA 按cDNA 反转录试剂盒使用说明进行反转录后，参照SYBR qPCR 荧光定量PCR 试剂盒，进行荧光定量PCR，每个基因3个生物学重复。

2 结果与分析

2.1 测序数据的统计和全基因组比对结果

共对8 个样本进行测序，对原始测序数据进行质控后，共获得了52.58 G 测序数据，352.82 M reads。其中矮秆材料组为168.30 M reads，高秆材料组为185.51 M reads，矮秆材料平均Q20% 为98.35%，最低为97.74%，平均GC%为45.04%，最低为41.64%；高秆材料组平均Q20%为98.67%，最低为98.57%，说明转录组测序结果质量较高，符合后续工作的要求。

对clean reads 进行去UMI 冗余处理和比对参考基因组，共获得了275.81 M UMI reads，样本平均比对到参考基因组的比对率为89.32%，最低为80.23%。

为了对RNA-seq 整体质量进行评估，通过基因覆盖率分析来表现测序结果的均一性。结果表明测序所得序列在基因上均匀分布，没有明显的偏向性，符合要求（图1）。

图1 基因覆盖率Fig.1 Gene coverage graph

2.2 差异表达基因分析和注释

通过计算FPKM 和TPM，共对101 041个基因进行差异水平分析。获得了75 004 个存在表达差异的基因，其中41 609 个存在显著差异，对数值进行修正后，获得了8665 个显著差异基因，其中上调基因2582个，下调基因6083个（图2）。

图2 差异基因散点图Fig.2 Differential gene scatters plot

2.2.1 差异基因GO 分析及富集分析 在GO 数据库中对甘蓝型油菜DW871 转录组中存在显著差异的基因进行功能比对，结果见图3。结果表示，在DW871 转录组中，8665 个显著差异基因被注释到BP（biological process，生物过程）、CC（cell component，细胞组分）、MF（molecular function，分子功能）3个大类共22 个亚类中。共有6875 条差异表达基因注释在BP 中，其中代谢过程（1904）、细胞过程（1574）、定位（421）注释到的差异基因最多。在CC中共注释到1442 个差异表达基因，细胞结构组件（cellular anatomical entity）、含蛋白质的复合物（protein-containing complex）、细胞内组分（intracellular）分别注释到1055、204 和183 个差异表达基因。在MF 中共注释到5376 个差异表达基因，结合（binding）、催化活性（catalytic activity）、转运活性（transporter activity）中的差异表达基因最多，分别注释到2781、2048、225、214和108个。

图3 上下调基因GO柱形注释Fig.3 GO column annotations for the up-and down-regulated genes

2.2.2 差异基因KEGG 注释及富集分析 为了对甘蓝型油菜DW871 矮化性状形成的生物学途径进行分析，以KEGG 数据库为参考，对DW871 显著差异基因进行通路（passway）统计和分类。结果显示，DW871 的差异基因共注释到5 个大类31 个亚类共295 个通路，共有6291 个差异表达基因（图4）。代谢中注释到差异表达基因最多，共3414个。其次是生物系统（1118）、环境信息处理（765）、遗传信息处理（575）、细胞过程（419）。按照代谢通路的筛选条件（1.2.3），筛选出苯丙烷生物合成、戊糖和葡萄糖醛酸相互转化、植物激素信号转导3个代谢通路，进一步分析其差异表达基因的表达模式，其富集程度分别为18.20%、14.60%、11.68%。

图4 KEGG富集Fig.4 KEGG enrichment

2.3 油菜株高性状相关基因的表达模式分析

在甘蓝型油菜矮秆和高秆对照组中，我们以野生型高秆为对照，在矮秆中筛选出了337 个与株高性状相关的差异表达基因，其中80个基因的表达上调，257个基因的表达下调。KEGG富集分析的结果显示，他们被注释到3 个代谢通路的13 个代谢途径中。根据前人研究，苯丙烷生物合成、果胶降解和生长素信号转导参与植物细胞木质素合成、细胞壁代谢和植物伸长等生理过程，对株高起十分重要的作用[25]。以此为基础，进一步筛选参与生长素信号转导、苯丙烷生物合成和果胶降解合成途径的关键基因，并对其差异表达方式进行比较分析。

2.3.1 木质素单体合成相关基因的表达模式分析 在苯丙烷代谢通路中包含两个重要分支，分别为苯丙烷生物合成途径和类黄酮代谢途径[26]。其中，苯丙烷生物合成途径的产物包括对香豆醇、针叶醇和芥子醇等木质素单体。我们分析了22 个参与木质素单体合成的显著差异基因（表1），它们分别编码苯丙氨酸解氨酶（PAL），反肉桂酸-4-单加氢酶（CYP73A），4-香豆酸酯-CoA 连接酶（4CL）、莽草酸/奎宁酸羟基肉桂酰转移酶（HCT）、肉桂酰辅酶A还原酶（CCR）、肉桂醇脱氢酶（CAT）、咖啡酸3-O-甲基转移酶（COMT）和咖啡酰辅酶A-O-甲基转移酶（CCoAOMT），矮秆转录组中除PAL、COMT和CCoAOMT中分别有一个基因表达上调以外，其他基因均表现为下调。

表1 甘蓝型油菜DW871中与木质素单体合成相关基因及其表达模式Table 1 Genes related to monolignol synthesis and their expression patterns in B.napus DW871

2.3.2 果胶降解相关基因的表达模式分析 果胶作为复杂的聚合物，其降解是次生壁木质化和细胞壁解体的起始。在矮秆转录组中共有43 个显著差异基因（表2），它们参与编码果胶酯酶（PME），多聚半乳糖醛酸酶（PG），半乳聚糖1,4-α 半乳糖醛酸酶（PGA），矮秆转录组中PME有20 个下调，12 个上调，PG均为下调，PGA有3个上调，1个下调表达。

表2 甘蓝型油菜DW871中与果胶降解相关基因及其表达模式Table 2 Genes related to pectin degradation and their expression patterns in B.napus DW871

2.3.3 生长素相关基因的表达模式分析 生长素是重要的植物激素，在植物细胞分裂、伸长和分化等方面发挥着重要的作用。在矮秆转录组中共有69 个显著差异基因（表3），它们参与编码生长素输入载体1（AUX1），生长素反应蛋白（AUX/IAA），生长素应答因子（ARF），生长素酰胺合成酶（GH3），生长素上调小RNA（SAUR），矮秆转录组中AUX1、AUX/IAA、ARF、GH3和SAUR中分别有1、2、1、1 和4个基因出现上调以外，其他均下调表达。

表3 甘蓝型油菜DW871中与生长素信号转导相关基因及其表达模式Table 3 Genes related to auxin signaling and their expression patterns in B.napus DW871

2.4 qRT-PCR

为了验证转录组测序所得基因序列结果的可靠性，本研究利用荧光实时定量RT-PCR，对测序中标记的显著差异基因进行表达量的测定与分析，包括PME、nsLTPs、DEFL、AUX1、ERTP、AEP等6 个基因，结果显示上述6 个基因的qRT-PCR 检测到的表达模式均与转录组测序所获得数据的表达趋势一致（DW871 vs HW871），验证了转录组数据的可靠性。

3 讨论

甘蓝型油菜作为世界三大油料作物之一，对保障国家食用油安全至关重要。DW871 作为甘蓝型油菜矮化品种，与以往报道的矮秆材料存在明显差异，具有重要的应用价值。本研究利用RNA-seq 高通量测序技术，构建DW871 的转录组文库，并以HW871 作为高秆对照，进行转录组数据的比较分析。从整体上分析了DW871基因的差异表达情况，对参与矮化性状的相关基因的表达进行了初步分析，挖掘了其中的关键基因，为进一步研究DW871的矮化性状提供理论基础和数据支撑。

图5 DW871显著差异基因qRT-PCR验证Fig.5 DW871 significant difference gene qRT-PCR verification

在植物中，苯丙烷生物合成途径又可分为公共途径和特异途径。在苯丙烷公共途径中，苯丙氨酸解氨酶（PALase）是催化该途径的限速酶，对木质素的合成起调控作用[27,28]。已有研究发现当拟南芥中表达PALase的基因活性降低时，会导致拟南芥的极度矮化和不育[29]。本研究中编码PALase 的基因在DW871 总体上调表达，证明PALase 在DW871 矮化性状的形成过程中起重要作用。莽草酸/奎宁酸羟基肉桂酰转移酶（HCTase）作为具有双重活性的酰基转移酶催化莽草酸/奎宁酸和对-香豆酰CoA 合成对-香豆酰莽草酸/奎宁酸。Sebastien 等发现在拟南芥中当编码HCTase的基因下调时，植株出现显著矮化，并且生长素信号转导受阻同样会导致HCT表达的下调[30]。在本研究中，编码HCTase 的基因均下调表达，且编码生长素输入载体AUX1 的基因总体下调，这说明在DW871中可能由于生长素信号转导受阻导致了HCT的下调，进而使植株出现了矮化表型。在木质素特异途径中，肉桂酰辅酶A 还原酶（CCRase）、肉桂醇脱氢酶（CADase）和咖啡酸3-O-甲基转移酶（COMTase）作为关键酶参与木质素单体的合成。其中CCRase 作为该途径的限速酶，与CADase 一同参与所有木质素单体的合成，并在木质素的合成中起主导作用[31]。Goujon 等构建拟南芥CCR反义表达载体。抑制AtCCR1的表达后，在木质素含量降低，结构和成分均发生了改变的同时，拟南芥相比于野生型更加矮小[32]。在以往的报告中，CADase 的主要影响在于对植株的机械强度和抗倒伏性状，Li 等在研究水稻突变体fc1，发现CAD表达的下调，导致了植株出现半矮表型[33]。COMTase 参与芥子醇的合成和G-木质素向S-木质素的转化，通过对植物木质素进行调控，影响植株的组织结构和形态。徐洋在研究小麦的COMT的反义表达载体时，发现COMT 被抑制后，会导致木质部细胞的变形，发育成熟的转基因植株表现矮化[34]。另外，Abbott等在研究CAD-COMT-CCR 抑制的转基因烟草，发现植株生长缓慢和植株矮化[35]。在本研究中，CCR、COMT和CAD 的表达总体下调，而在前期研究中，我们初步定位到了造成DW871 矮化的候选基因，其表达同样为下调，该基因通过编码DIR1蛋白，参与木质素单体偶联形成二聚体，进而参与木质素的生物合成[36,37]。推测因为上述基因的表达发生改变，使得木质素的组成和结构发生改变，进而影响DW871 的株型，这与DW871 的矮化表型密切相关（图6）。

图6 木质素合成途径Fig.6 Monolignol synthetic pathway

在植物果胶降解代谢途径中，果胶酯酶（PMEase）起始第一步反应，与聚半乳糖醛酸酶（PGase）一同调控果胶的降解，并在细胞壁的合成和降解中发挥着重要作用[38～40]。已报道的研究中发现，PMEase 可以提高PGase 的活性，增强细胞壁的机械强度，从而影响茎的伸长等生理生化过程[41,42]。Hong等发现在OsPMEI28过表达的水稻转基因植株中，当PMEase 受到抑制时，细胞壁的理化性质发生了改变，茎的伸长出现了异常，水稻出现了矮化表型[43]。在本研究中，PG在DW871 抽薹期发生了显著下调，而PME也发生了显著变化，这可能说明PG和PME表达的变化导致了果胶无法正常降解代谢，而积累的聚半乳糖醛酸通过一系列的反应，加强了细胞壁的结构，抑制了细胞壁的降解，阻碍了细胞的伸长，最终导致DW871矮化性状的出现。

植物的生长素信号转导途径主要包括生长素输入载体（AUX1）、生长素受体蛋白（TIR1）、生长素应答因子（ARF）和生长素反应元件（AUX/IAA、GH3、SAUR）4 个部分。AUX1 承担着生长素输入目标细胞的任务，ARF 是生长素发挥作用的媒介。张赛娜等在研究水稻突变体OsARF19时，发现OsARF19的超表达导致纤维素含量的降低、细胞缩小，叶片变窄、矮化[44]。在本研究中，ARF的表达显著上调，可能导致了细胞发育异常，细胞的伸长被抑制。AUX/IAA、GH3和SAUR是三种最主要的生长素早期响应基因。其中，AUX/IAA 和ARF 之间存在互作，AUX/IAA 和GH3 与植物的长度密切相关并介导光反应，以往有研究报道AtARF8参与子叶的伸长，白光、蓝光和红光等外源光可影响拟南芥的生长发育，ydk1-D参与生长素活性的调控，能够降低植株的顶端优势，当处于蓝光和远红光下时，它的表达下调[45,46]。在本研究中，AUX/IAA和GH3的表达总体下调，这可能是当光热条件发生改变时，DW871的株高出现较大变化的原因之一[47]。

4 结论

通过对DW871 进行转录组测序和差异基因功能注释，挖掘到了在果胶降解代谢途径、苯丙烷生物合成途径以及生长素信号转导途径等代谢通路中起关键作用的基因，并分析了它们在油菜生长过程中的表达模式。而在DW871的前期研究中，发现DW871 的矮化表型与赤霉素的生物合成无关。另外，我们初步定位到了导致DW871矮化表型的候选基因，推测该基因可能通过合成DIR1与生物氧化酶共同调控木质素的生物合成。通过上述结果，我们推测，生长素应答和木质素合成的异常，造成了细胞伸长、生长和细胞壁变化的异常，并最终导致了DW871的矮化。