基于转录组测序的花烟草低温胁迫响应转录因子挖掘

鲁 琳 杨尚谕 刘维东 鲁黎明*

（1. 四川农业大学风景园林学院，成都 611130；2. 四川农业大学农学院，成都 611130）

植物在生长过程中，往往会遭遇到干旱、盐碱、高温以及冷害等不利环境的胁迫，从而生长发育受阻，对农作物而言，则会导致产量的损失及品质的下降［1～2］。其中，低温冷害是影响农业生产正常进行的最为关键的逆境之一。低温会造成细胞严重脱水，当胞间溶液结冰时，还会导致细胞的机械损伤，植物则表现出生长迟缓、萎蔫、弱苗，最终影响作物产量与品质形成［3］。

在长期的进化过程中，植物形成了一套复杂的抵御冷害的机制。在冷害来临时，植物会发生一系列生理生化变化，如细胞膜组成改变、物质转运活性增强、活性氧（ROS）清除加快及细胞质保护物质合成增加等，以适应寒冷的环境，提高生存能力［4～6］。在此过程中，包含着大量的冷调节基因（Cold regulated gene，COR）表达水平的改变。在分子层面上，从植物感受到冷胁迫，到COR基因表达水平发生变化，蕴藏了复杂的信号转导途径，构成了植物耐冷机制的重要组成部分。

COR 是一大类基因，在植物的耐冷机制中扮演着十分重要的角色，参与了植物的代谢、蛋白质合成、信号转导、转录调控、激素合成等诸多的生物学过程［4，7］，其中，就包含众多的转录因子。研究表明，DREB/CBF 转录因子家族基因，是植物中普遍存在的参与耐冷调节的最为关键的因子［8～11］。拟南芥DREB1b 过表达植株表现出耐冷性增强［12～13］，而DREB1a、DREB1b 及DREB1c 三突变体植株，则显示出对冷极度敏感［14］，表明这3 个基因对拟南芥的耐冷性至关重要。除此之外，其他转录因子也在植物的耐冷调节机制中占有重要的地位，例如，NAC 类转录因子TaNAC2［15］、AtLOV1［16］、SiNAC1［17］；WRKY 类 转 录 因 子AtWRKY34［18］、Os⁃WRKY45［19］、OsWRKY71［20］及CsWRKY46［21］；MYB 类转 录 因 子OsMYB4［22］；bZIP 类 转 录 因 子Os⁃bZIP52［23］；ZFPs 类 转 录 因 子 SCOF1［24］及GmZF1［25］等。

转录水平的调控是植物生命活动调节的重要调控形式，转录组测序则是了解植物基因总体表达情况的强大工具［26～27］。目前，高通量转录组测序技术已经在植物响应生物及非生物胁迫的机制研究中得到了广泛的应用［26，28～32］，一大批涉及防御、物质运输、信号转导及次生代谢的基因得到鉴定，对于揭示植物抗逆的分子机制起到了很大的推动作用。

花烟草（Nicotiana alata）是属于茄科（Solana⁃ceae）烟草属（Nicotiana）的一年生野生草本植物，由于其花色艳丽多样、花期较长，在园林绿化中得到了较多的应用。然而，当环境温度低于15℃时，花烟草就不能够正常开花。为了探索花烟草抵御低温的机理，本研究采用转录组测序（RNA-Seq）技术，在转录组水平上分析了花烟草响应低温胁迫的基因表达特性，以期为花烟草耐冷的分子机制的研究，及耐低温花烟草品种的选育提供借鉴。

1 材料与方法

1.1 试验材料

以花烟草（Nicotiana alata）为试验材料，由四川农业大学烟草实验室提供。

1.2 材料处理及转录组测序

幼苗材料培养：挑选籽粒饱满且大小一致的种子，用次氯酸钠溶液表面消毒后，用无菌水清洗数次，随后播种于MS 培养基上，并置于恒温培养室中生长（光照：16 h光/8 h暗；温度：25℃）。

低温处理及总RNA 提取：待烟苗生长4 周后，挑选72 株长势一致的烟苗，分为6 组，在光照培养箱中进行低温（4℃）处理，每组烟苗为1个重复，每重复12 株烟苗。在处理后0 及12 h 后，分别进行整株取样，并置于液氮中保存，用于总RNA 的提取。总RNA 的提取参照鲁黎明等［33］的方法进行。提取的总RNA，经紫外分光光度计和琼脂糖凝胶电泳检测合格后，送北京诺禾致源科技公司进行转录组测序。

1.3 测序数据的分析及功能注释

植物样品的测序，由诺禾致源公司采用Illu⁃mina Hi Seq2500 系统进行。对测序所获得的原始数据，按照诺禾致源公司的数据处理方法，进行原始数据的过滤，获得高质量的clean reads，并进行后续的分析工作。

本研究采用无参转录组测序，因此，使用Trin⁃ity 软件进行测序数据的拼接。然后，根据拼接的结果，进行基因表达丰度的定量分析。

通过BLASTX 软件，与NCBI 的non-redundant（NR）、UniProt/ Swiss-Prot、Gene Ontology（GO）、KEGG（Kyoto Encyclopedia of Genes and Gnomes）等数据库进行比对，以获得其功能注释信息。

1.4 差异基因的认定

差异表达基因（DEGs，differentially expressed genes）的筛选在低温处理与对照的样本间进行。筛选的条件为P-value≤0.01，|log2Fold Changes|>1。其中，Fold Changes 代表不同处理之间基因表达量的比值。该值大于1 时，则该基因认定为上调表达，反之，则为下调表达。

1.5 差异表达基因的GO及KEGG富集分析

差异表达基因的GO 及KEGG 富集分析，根据叶兴状等［34］的方法，分别在GO 及KEGG 数据库中进行。

1.6 差异表达基因的验证

为了对转录组测序的结果进行验证，随机挑选6 个基因，并利用Primer 5.0 软件设计扩增引物（见表1），采用qPCR的方法，对其表达水平进行验证。将各个处理的样品总RNA 反转录合成cD⁃NA，以此为模板，进行qPCR 反应，内参为β-actin。qPCR 的反应体系、反应程序以及基因相对表达量的计算，均参照鲁黎明等［33］的方法进行。

2 结果与分析

2.1 测序质量分析

采用高通量测序平台，对低温处理组及对照组的样品进行了RNA 测序。结果显示，共获得了47 599 094～57 845 042 的clean reads（见表2）。就每个样品的GC 含量来看，其含量在39.37%～43.03%，偏差很小，而且基本呈随机分布。此外，测序的质量还可以用质量控制参数Q30 的数值的大小来判断。一般要求该值不低于80%。在本研究中，Q30 的数值均大于89%，说明本次测序的质量非常可靠（见表2）。样品重复间数据之间的相关性，也是衡量测序质量的重要指标之一。本研究样品重复间数据的相关性分析结果表明，重复之间的决定系数R2，最低为0.962，最高则达到了0.982，说明重复之间数据的重现性较好，可以进行下一步的分析。

表1 qPCR反应所用引物Table 1 Primers used for qPCR

表2 各样品测序质量Table 2 Evaluation of sequence quality

2.2 差异表达基因统计分析

将低温处理12 与0 h 样品的转录组数据进行对比分析，以|log2fold changes|≥1，以及P-value≤0.01 为筛选条件，进行差异表达基因的寻找。结果表明，低温处理12 h 后，差异表达基因的数量达到了8 388 个，其中，上调表达4 229 个，下调表达4 159个。

2.3 差异表达基因的GO功能注释分析

对差异基因所进行的GO 分析结果显示（见表3），上调的差异表达基因的功能，主要集中在生物过程（Biological process）以及分子功能（Molecular function）方面。在生物过程中，Phosphate-contain⁃ing compound metabolic process（761，18.0%）及Phos⁃phorus metabolic process（761，18.0%）富集上调表达基因最多，其次为Macromolecule modification（683，16.2%）、Cellular protein modification process（634，15.0%）和Protein modification process（634，15.0%）。在分子功能中，Transferase activity（1 558，36.8%）富集上调表达基因最多，其次为Transferaseactivity，transferring phosphorus-containing groups（948，22.4%）与Kinase activity（638，15.1%）。

与上调差异表达基因不同，下调差异表达基因的功能主要集中在细胞组分（Cellular compo⁃nent）及分子功能（Molecular function）方面。细胞组分中，Ubiquitin ligase complex 占比最多（125，3%）；分子功能中涉及到的两个GO 条目，均为泛素蛋白转移酶活性相关，其占比均为3.1%（见表4）。

2.4 差异表达基因的KEGG通路分析

对差异表达基因进行了KEGG 通路分析，结果如图2所示。

上调差异表达基因富集到20条主要的代谢通路。富集因子最高的为“Carbon fixation in photo⁃synthetic organisms”途径，其次为“Carbon metabo⁃lism”和“Glycolysis/ Gluconeogenesis”途径。而富集差异基因最多的途径为“Metabolic pathways”，富集到了807 个差异表达基因，其次为“Biosynthesis of secondary metabolites”途径，有458 个差异表达基因被富集到该途径。

下调差异表达基因显著富集的通路数量同样是20 条，其中，富集因子最高的通路为“alpha-Lin⁃olenic acid metabolism”。在富集基因的数量方面，这20 条通路所富集差异表达基因的数量，较上调基因少了很多。富集下调差异基因最多的通路为“Protein processing in endoplasmic reticulum”，仅仅富集了89个基因。其次，为“Spliceosome”途径，富集的差异表达基因数为74个。

从富集差异表达基因的数量角度也可以看出，低温处理后，花烟草代谢途径中，上调基因的数量远远高于下调基因。从而在某种程度上说明，花烟草在遭受低温胁迫后，其生理响应过程主要以上调相关基因的表达为主。

2.5 花烟草低温响应的转录因子分析

转录因子通过与基因上游启动子区域的作用元件结合，调控基因的表达，从而影响到植物的生理进程。为了探索低温胁迫下花烟草转录因子表达的变化，对差异表达基因中的转录因子进行分析。结果表明，在花烟草受到低温胁迫后，共有118 个转录因子的表达发生了显著改变，其中，上调表达82 个，下调表达36 个。这些转录因子属于28 个家族，其中，数量最多的为NAC（19），其次为ERF（16）、MYB（15）、WRKY（15）等（见图3）。

表3 上调差异表达基因的GO功能注释分析Table 3 Go functional annotation analysis of up regulated differentially expressed genes

表4 下调差异表达基因的GO功能注释分析Table 4 Go functional annotation analysis of down regulated differentially expressed genes

2.6 qRT-PCR 验证

为了验证RNA 测序结果的可靠性，随机选取了IAA17、NAC2、P450、WRKY6、CIPK6、ERF025 等6 个基因，进行了qPCR 扩增，计算了其表达量，然后，与其测序结果进行对比（见图4）。从图4 可以看出，在所选择的6 个基因中，除了P450 外，其余4 个基因的表达趋势，两种方法均表现一致，从而说明了本研究RNA 测序结果较为可靠，同时也表明，对RNA 测序结果所做的相关分析较为客观。

3 讨论

低温是植物在生长过程中所遭受到的主要的环境胁迫之一。在低温胁迫下，植物体内会发生一系列生理生化变化，以抵御低温对植物生长的影响。其中，转录因子能够通过调控下游基因的表达，调动相关代谢途径产生级联反应，从而提高植物对低温的适应能力［35］。在本研究中，花烟草在低温刺激下，有28 个家族的118 个转录因子的表达发生了显著变化，表明低温导致了花烟草众多的生理进程发生改变。

3.1 NAC 家族转录因子

NAC 转录因子家族是植物最大的转录因子家族之一，其名称来源于NAM、ATAF1/2 及CUC2 的首字母缩写。在结构上，NAC 家族转录因子N 端结构域较为保守，而其C 端的转录激活域，则由不同的相对保守区域构成。这些结构上的特点，使得NAC 转录因子的功能多样，广泛参与了植物的生长发育及非生物胁迫过程［36］。NAC 转录因子的启动子区域，包含有低温响应元件，表明其参与了植物低温的胁迫响应［36］。研究表明，小麦TaNAC2［15］、番茄SINAC1［17］、番茄SINAM1［37］等基因的过表达，均提高了转基因植物抵御低温的能力。此外，香蕉的MaNAC1 能够通过调控ICE1-CBF 信号途径，从而增强香蕉的抗冷性［38］；而苹果的Md⁃NAC029 则通过抑制MdCBF1 和MdCBF4 的表达，负调控苹果的耐寒性［39］。本研究中，花烟草在低温胁迫下，17 个NAC 转录因子的表达显著改变，其中，包括12 个上调表达，7 个下调表达（见图3）。由此说明，花烟草NAC 转录因子参与低温响应调控时，可能存在正向与负向两种方式，其中，具有正向调节作用转录因子占多数。

3.2 ERF 家族转录因子

AP2/ERF 转录因子家族在植物中广泛存在，分为AP2、RAV、ERF和DREB 4个亚类，其中，ERF和DREB 在AP2/ERF 转录因子家族中占有最重要的地位。研究证实，ERF 转录因子通过结合下游基因启动子上的GCC 盒（GCCGCC），从而调控其表达［35］。在植物中，DREB/CBF途径是关键的低温信号转导与调节响应途径，其中，ICE-CBF 途径又是DREB/CBF 反应途径中的关键一环［40］。当低温胁迫来临时，ICE 与CBF 启动子区域的MYC 元件（CACATG）相结合，上调CBF 基因的表达；而CBF则结合低温响应基因的顺式作用元件，增强这些基因的转录，从而改善相应的生理活动，以适应低温环境［40］。研究发现，番茄TERF2 和LeERF2［41］、草 地 禾PpCBF3［42］、麻 风 树JcCBF2［43］、小 麦TaDREB2 和TaDREB3［44］及水稻OsDREB6［45］的 过表达植株的耐寒性均得到增强。本研究的结果也显示，花烟草ERF 转录因子对低温环境有显著响应，其中，有13 个上调表达，3 个下调表达（见图3）。在以往的研究中，ERF 类转录因子下调响应低温胁迫还鲜有报道。所以，本研究的结果既说明了植物响应低温胁迫的机制的复杂性，也为其研究提供了有益的线索。

3.3 MYB 家族转录因子

在植物的转录因子家族中，MYB 是较大的一类，包含的成员数量最多。MYB 名称来源于其保守的DNA 结合域，其所结合的调控基因的顺式作用元件为TAACTG。研究表明，MYB 转录因子参与了植物的低温响应调控过程，但是，不同的MYB成员所起的作用并不相同［35］。例如，过表达Os⁃MYB4 的拟南芥［46］、过表达MdSIMYB1 的烟草对低温胁迫的耐受力均增强［47］。此外，梨PcMYB10 对低温下花青素的合成有正向作用［48］。然而，拟南芥的AtMYBC1［49］以及MYB15［50］却负向调控植物的低温耐受性。本研究也得出了类似的结果，在低温胁迫下，花烟草11 个MYB 基因的表达显著增强，有4 个则显著下调（见图3）。这些结果暗示，MYB 转录因子广泛参与了植物的低温胁迫响应，其成员间发挥作用的机制并不相同。

3.4 WRKY 家族转录因子

在高等植物中，WRKY 转录因子家族是最大的一类，其所结合的靶基因的作用元件为TGAC⁃CA/T（W-BOX）。据报道，WRKY 转录因子通常抑制下游靶基因的表达，例如，低温处理后拟南芥AtWRKY34 上调表达，并负向调控拟南芥成熟花粉的低温敏感性［18］。与此相反，黄瓜的CsWRKY46则是一个植物耐寒性提升的正向调节因子，其在拟南芥的过表达提高了转基因植株的低温耐受性［21］。与此相类似，花烟草在低温环境下有11 个WRKY上调表达，4个下调表达（见图3），说明在花烟草的低温胁迫响应中，WRKY 转录因子可能也起着十分重要的作用，其作用的机制有待于进一步的探索。

3.5 GATA 家族转录因子

植物GATA 转录因子也是一类较为重要的转录因子，在植物矿质营养、形态建成、信号转导以及花器官发育等生理生化过程中占有重要地位，其名称源自其特异结合的元件WGATAR（W 为T或A，R 为G 或A）［51］。GATA 转录因子参与植物冷胁迫的报道较少，仅有的研究表明，拟南芥的两个GATA 转录因子GNC 和GNL，能够与SOC1 相互作用，从而正向促进拟南芥低温耐受能力［52］。此外，班巴拉花生（Vigna subterranea）的GATA9［53］以及构树（Broussonetia papyrifera）的GATA 类部分转录因子也参与了植物的低温胁迫响应［54］。本研究的结果也证实了前人的研究结论，花烟草有8 个GATA类转录因子参与了低温胁迫的响应（见图3）。然后，有所不同的是，这8 个GATA 中，有7 个是下调表达，只有1个上调表达。这种现象在一定程度上说明，花烟草GATA 转录因子参与低温胁迫响应的方式，具有其独特的特性。

3.6 bHLH 家族转录因子

bHLH 转录因子家族，在高等植物转录因子家族中位列第二位，其名称来自于其bHLH 结构域。该结构域十分保守，由大概60个氨基酸所构成，包含了碱性区域和螺旋—环—螺旋结构域两个功能完全不同的区域［55］。与植物的其他转录因子类似，植物的bHLH 转录因子参与了信号转导、物质合成、腺毛和根毛发育及逆境胁迫等［56］。bHLH 转录因子的一个重要的功能，是通过植物中普遍存在的ICE1-CBF 冷响应通路，参与植物对低温胁迫的信号感应及信号传递。拟南芥ICE1 与ICE2 均能够激活CBF1 的表达，通过调控下游更多的冷响应表达，提高其对低温胁迫的适应能力［55］。水稻的OsICE1 和OsICE2［10］、白菜的BrICE1［57］、野生葡萄 的VabHLH1 以 及VvbHLH1［58］均 与 拟 南 芥 的ICE1有着相似的作用与功能。花烟草可能也存在着类似的响应途径。在冷胁迫下，有6 个bHLH 的表达产生了明显的响应，其中，有4 个上调表达，2个下调表达。一般来说，bHLH 的过量表达能够提高植物的低温耐受能力。然而，bHLH 如何以负调控的方式参与低温胁迫响应，目前还不清楚。这种现象一方面说明冷胁迫响应在植物中存在的普遍性，另一方面也表明了花烟草冷胁迫响应的特殊性与复杂性。