萱草叶片响应低温胁迫的转录组分析

黄东梅 陈 颖 白 露 倪迪安 徐奕扬 张志国 秦巧平

（上海应用技术大学生态技术与工程学院，上海 201418）

萱草（）为阿福花科（Aspho⁃delaceae）萱草属（）的多年生宿根草本植物，是一种深受大众喜爱的植物，被广泛种植用于药用和观赏等方面。近年来，萱草成为研究者的重点关注对象之一，萱草花发育、衰老和遗传机理研究也成为热点。低温胁迫是植物生长过程中经常会遭遇的一种灾害，是限制农作物生长发育、生产力及地理分布的重要非生物胁迫之一，影响植株的生长和发育，导致产量的损失。但是也有部分植物能够通过调节自身生理和生化的反应来提高对低温的耐受力，从而减轻甚至消除低温胁迫所带来的伤害，包括细胞质膜的稳定、低温保护性多肽的合成及可溶性糖、可溶性蛋白质和脯氨酸的积累。但事实上，植物从接受低温信号到引起生理生化的响应，再到产生抗寒能力，这整个过程存在着复杂的信号传导系统，需要调节大量相关基因的表达。

转录组测序技术（RNA-seq）可以高通量、快速地检测基因差异表达、基因特异性表达，全面地揭示逆境胁迫下整个基因组水平的表达情况，从而探明植物在逆境胁迫下的相关应答途径，目前在许多物种上已经有了相关的报道，如拟南芥（）、狗牙根（）、玉米（）等。Niu 等通过转录组测序鉴定出了24 个基因家族与桃（）冷应答相关，其中相关程度最高的基因家族是/，其次是和基因家族。Hao 等通过模拟自然气候变化建立冷适应条件，发现茶树（）转录组在冷驯化的初始阶段就开始变化发生，而与细胞壁相关基因的变化贯穿了整个冷驯化过程。拟南芥转录组学分析结果表明CBF2 和ZAT12 转录因子在低温胁迫下冷应答途径中发挥着重要作用，且除常见的基因诱导外，基因的抑制在植物冷驯化中也发挥不可或缺的作用。Hong等通过转录组学分析了水稻（）参与冷应激反应和昼夜节律模式基因在体内差异表达所产生串扰的分子调控网络，发现参与光合作用的光捕获复杂蛋白与串扰关系紧密。

目前，萱草转录组学研究虽有报道，本课题组李昊通过转录组测序筛选出了萱草抗旱品种‘猛子花’在干旱条件下较敏感品种‘害羞老虎’高表达水平的140个基因，Li等通过将转录组测序与EST-SSR 标记结合，探讨了萱草黄花菜的遗传多样性及在群体结构中的应用。但低温条件下萱草叶片转录组学分析鲜见报道。本研究通过高通量测序技术对萱草低温处理组和对照组叶片进行转录组测序，探索萱草低温胁迫分子响应机理，以期为耐寒萱草种质开发、分子育种提供一定的理论支撑。

1 材料与方法

1.1 材料和处理

以 萱 草‘阿 斯 隆’（‘Ath⁃lone’）为材料，根据多年的观察，萱草‘阿斯隆’是长三角地区栽培萱草中耐寒性较强的品种之一，可作为研究常绿萱草潜力品种。采集长势良好且生长状况一致的萱草叶片（15 ℃，设为对照组T），用蒸馏水洗净晾干，液氮速冻并于−80 ℃下储存用于后续试验。后将萱草移至室内恒温培养箱，在10、5、0 ℃下培养24 h 后采集样品，分别命名为低温处理组T、T、T组。取样部位为完全展开的功能叶片，每组设置了3次生物学重复。

1.2 RNA的提取及cDNA文库的构建

使用改良的CTAB 法提取叶片总RNA，用1%的琼脂糖凝胶电泳分析RNA降解及污染程度，NanoDrop 检测RNA 的纯度、浓度、核酸吸收峰，安捷伦2100 RNA Nano 6000 Assay Kit 精确检测总RNA 的完整性和浓度，−80 ℃保存备用。用带有Oligo（dT）的磁珠富集萱草mRNA，加入fragmenta⁃tion buffer 将mRNA 进行随机打断，再以片段化的mRNA 为模板，采用M-MuLV 第一链cDNA 合成试剂盒合成cDNA。利用AMPure XP beads 纯化cD⁃NA，选择合适大小的片段进行PCR富集得到cDNA文库，再用Agilent 2100对cDNA文库进行质控。

1.3 转录组测序及组装

将检测合格的样品提交Illumina HiSeq 平台进行测序（上海中科新生命生物科技有限公司）。对测序得到的原始测序数据（Raw Reads）通过去低质量序列、接头污染及N 比例大于5%的reads等过程完成数据处理得到高质量的有效测序数据（Clean Reads）。再采用Trinity 对clean reads 进行组装成Unique序列。

1.4 基因功能注释

将测序获得的转录组进行了五大数据库的基因功能注释。利用NCBI BLAST2.6.0+（Evalue<10）在NCBI 的NR（non-redundant protein data⁃base）数 据 库、Swiss-port（swissprot protein se⁃quenced database）数据库对基因功能进行注释；蛋白家族（Pfam）使用HMMER3.0（Evalue<10）分类；基因本体（Gene Ontology，GO）和京都基因与基因组百科全书（kyoto encyclopedia of genes and ge⁃nome，KEGG）注释分别使用Blast2Gov2（Evalue<10）和KAAS（Evalue<10）软件结合BLASTX 对Nr数据库的注释进行。

1.5 差异表达基因的筛选及相关分析

将低温处理组T、T、T分别与对照组T进行比较，以|logfold-change|≥2，-value≤0.05 为阈值进行差异表达基因（Differentially expressed genes，DEGs）的筛选。随后在1.4 中筛选出差异表达基因的GO 和KEGG 注释，并进行分类。再用jvenn（http：//www. bioinformatics. com. cn/static/others/jvenn/example.html）绘制DEGs 韦恩图，用heatmap⁃per（http：//www.heatmapper.ca/expression/）绘制3 个低温处理组共有DEGs 在低温胁迫过程中的表达变化，用String（https：//string-db.org/）和Cytoscape进行基因调控网络分析。

1.6 基因表达的qRT-PCR验证

为了验证萱草转录组数据，从差异表达基因中随机选择12 条基因进行qRT-PCR。利用Prim⁃er5（https：//sg. idtdna. com/pages/tools/primerquest）网站根据所选12 条基因的cDNA 序列设计定量表达分析所用引物。基于本课题组前期对6 个实时荧光定量PCR 内参基因的筛选分析，UBQ 在萱草‘阿斯隆’中稳定性较好，用作本研究的内参（见表1）。提取叶片总RNA，检测合格后采用MMuLV 第一链cDNA 合成试剂盒进行第一链cDNA的合成。qRT-PCR 检测在荧光PCR 定量仪中进行，反应条件为95 ℃预变性5 min；95 ℃变性10 s，60 ℃退火30 s，循环40 次；熔解曲线采集程序为95 ℃15 s，60 ℃60 s，95 ℃15 s。采用2法计算基因相对表达量，每个样品重复3次。

表1 差异表达基因qRT-PCR引物Table 1 qRT-PCR primers for 12 DEGs

2 结果与分析

2.1 萱草转录组数据的组装与数据评估

构建了低温处理组（10、5、0 ℃）和对照组（15 ℃）萱草叶片cDNA文库，去除低质量数据后，四组分别平均获得46 376 031、42 120 296、41 392 792、44 378 448 条Clean Reads，6.44GB、5.86GB、5.74GB、6.16GB 的高质量数据，平均GC 含量分别为45.47%、45.17%、45.35%、45.88%，Q20平均值分别为98.24%、98.18%、98.29%、98.30%，Q30平均值分别为94.57%、94.45%、94.68%、94.73%（见表2），表明转录组测序数据符合后续试验要求。

表2 萱草转录组测序数据评估Table 2 Evaluation of transcriptome sequencing data in H.fulva

低温处理组和对照组萱草叶片样本转录组Clean Reads 经Trinity 软件组装后共获得≥200 bp的Unigene 143 636 条，其中最大转录本长度为9 494 bp，最小转录本长度为201 bp，平均长度为602 bp，N50 值为936 bp。图1 为组装转录本的长度分布，其中长度≥1 kB 的转录本共计22 326 条，占总转录本数的15.54%。Unigenes 注释同源基因的物种分布显示，相似序列较高的物种中占比例最高的为芦笋（），其次为油棕（）、栓皮栎（）和海枣（）。

图1 3个低温处理组和对照组萱草转录组unigenes平均长度分布Fig.1 Average unigenes length distribution of transcriptome in three low temperature treatment groups and control group

2.2 差异表达基因的筛选及分类

将3个低温处理组（T、T、T）样品中表达的基因与对照组T的基因进行比较显示，随着温度的下降，上调差异表达基因减少，下调差异表达基因增多。以FDR<0.001，|logfold-change|≥2，≤0.001为筛选标准，共筛选出差异表达基因2 457个，其中上 调 基 因1 253 个，T、T、T低 温 处 理 组 分 别 为1 038、385、112个；下调基因1 204个，T、T、T冷处理分别为697、532、315个（见图2）。

图2 DEGs韦恩图Fig.2 Venn diagrams of DEGs

对DEGs 进行GO 富集分析，将差异表达基因分为生物学过程（biological process，BP）、细胞成分（cellular component，CC）和分子功能（molecular function，MF）3 大类（见图3）。结果显示，上下调差异基因的GO Terms 基本一致。差异表达基因涉及生物学过程21 种，其中占比最大的为细胞过程（cellular process，上调29.46%、下调33.49%），其次是代谢过程（metabolic process，上调26.95%、下调30.72%）和单组织过程（single-organism process，上调22.63%、下调24.58%），不同的是光合作用（photosynthesis）仅有上调基因涉及，而内膜系统（endomembrane system）为下调基因特有。在该分类中有对刺激的反应（response to stimulus）、对非生物刺激的反应（response to abiotic stimulus）、细胞对刺激的反应（cellular response to stimulus）、对温度刺激的反应（response to temperature stimulus）以及对寒冷的反应（response to cold）等冷应激相关的生物学过程。在细胞成分方面，上调差异基因涉及11 种，下调10 种，相差的Term 为捕光复合物（light-harvesting complex）。富集中细胞部分（cell part）和 细 胞（cell）占 比 最 大，均 为 上 调45.62%、下调46.14%；其次是细胞器（organelle，上调29.58%、下调27.95%）。上调差异基因涉及分子功能17 种，下调15 种，相差的2 个Term 为过氧化物酶活性（peroxidase activity）和天冬氨酸型内肽酶活性（aspartic-type endopeptidase activity）。细胞成分占比较大的前3 项为催化活性（catalytic ac⁃tivity，上调27.71%，下调35.90%）、结合（binding，上调27.30%，下调23.86%）和转移酶活性（transfer⁃ase activity，上调15.21%，下调15.54%）。

图3 DEGS的GO富集分析Fig.3 GO term enrichment analysis of DEGs

KEGG代谢通路分析结果表明，上调基因共富集到41 条代谢途径，下调基因39 条（见图4）。这些代谢途径可分为5大类：代谢、遗传信息处理、环境信息处理、细胞过程以及有机系统。基因数较多的几条为代谢途径（metabolic pathways）、次生代谢产物的生物合成（biosynthesis of secondary me⁃tabolites）、光合作用—触角蛋白（photosynthesis-an⁃tenna proteins）、植物激素信号转导（plant hormone signal transduction）和苯丙素生物合成（phenylpro⁃panoid biosynthesis）。从整体来看，下调基因的富集程度较上调基因更高，如在代谢途径（metabolic pathways）、生物合成（biosynthesis of secondary me⁃tabolites）、内质网的蛋白质加工（protein processing in endoplasmic reticulum）和-亚麻酸代谢（alpha-Linolenic acid metabolism）等下调基因均比上调基因占比高。

图4 DEGS的KEGG代谢通路分析1.代谢；2.遗传信息处理；3.环境信息处理；4.细胞过程；5.有机体系统Fig.4 KEGG analysis of DEGs 1.Metabolism；2.Genetic information processing；3.Environmental information processing；4.Cellular processes；5.Organismal systems

2.3 与植物激素信号转导相关的DEGs

对KEGG 富集通路“植物激素信号转导”相关差异基因表达模式进一步分析。发现在生长素（Auxin，IAA）信号通路中，Aux/IAA 编码的等3 个DEGs 在T和T阶段表达量均呈现下调趋势，而后在T阶段略有回升；SAUR 编码的3 个DEGs呈现相反表达模式，在T阶段表达量上升，为T阶段的5.894 倍，表达量下调，为T阶段的0.227 倍；SAUR 编码的在T阶段表达量高至T阶段的13.623 倍（见图5）。在细胞分裂素（cytokinine，CTK）信号通路中，A-ARR 编码的2 个DEGs 中，在T阶段表达量下调，T和T阶段表达量相对T阶段上调，甚至T阶段时表达量超过T阶段；另一个基因T阶段表达量上调，T和T阶段表达量下调至T阶段与T阶段持平。在茉莉酸（Jasmonic acid，JA）信号通路中，JAZ编码的等4 个DEGs 中，4 个在T阶段基因表达均显著下调趋势，而后等3个基因在T和T阶段表达量呈上升趋势。

图5 植物激素信号转导通路及其DEGs的表达模式Fig.5 Expression patterns of DEGs in plant hormone signal transduction pathway

2.4 与可溶性糖合成相关的DEGs

萱草低温胁迫下可溶性糖代谢通路相关差异基因表达模式分析表明，基因在半乳糖代谢通路、淀粉和蔗糖代谢通路中均有参与（见图6）。参与“半乳糖代谢”通路中的等5 个DEGs 表达量随着温度的下降均呈现出先上升后下降的趋势，其中在T阶段比T阶段高28.114 倍（见图6A）。“淀粉和蔗糖代谢”通路中，基因表达量随着温度的下降呈现出先下后上的趋势，和基因表达量持续下降，而另外等6 个DEGs 表达量则先上升后下降（见图6B）。“果糖和甘露糖代谢”通路中，4 个DEGs（−3、、、）表达量均在T阶段达到顶峰（分别为T阶段的6.135、4.900、4.215、5.013 倍），而后在T和T阶段呈下降趋势；另4 个DEGs（、-2 个）则在T阶段表达量最低（分别为T阶段的0.315、0.223、0.076、0.224 倍），而在T和T阶段呈上升趋势（见图6C）。整体来说，在低温胁迫下，参与可溶性糖代谢的DEGs中，大多数基因的表达量较对照组高。

图6 可溶性糖代谢通路及其DEGs的表达模式A.半乳糖代谢通路；B.淀粉和蔗糖代谢通路；C.果糖和甘露糖代谢通路Fig.6 Expression patterns of DEGs in soluble sugar metabolic pathways A.Galactose metabolism pathway；B.Starch and sucrose metabolism pathway；C.Fructose and mannose metabolism pathway

2.5 低温诱导后共有的DEGs

韦恩图（见图2）表明，在2 457个差异基因中，T、T、T低温处理组与对照组T相比共有的差异基因有29个，其中上调基因16个，下调13个（见图7）。平均表达值热图分析结果表明，29 个基因的表达模式聚集为4 个分支。分支A 和B 具有相似的表达模式，这些分支里含有-（）、-（）、/-（）、-等基因，这些基因在低温条件下表达量受抑制。分支C 和分支D 的表达模式相反，随着温度的下降，位于分支C 的（）、（）、-（）等基因的表达量表现为先上升后下调，而分支D 内的-（）、-（）、（）等基因表现为先上升后下调再有所上调。-（）基因在受到低温胁迫时表达量下调，被聚类为A 分支，而另一个UDP 糖基转移酶基因-（）却表现为上调，聚类位于C分支。

图7 3个低温处理组共有DEGs在低温胁迫下的热图Fig.7 Heat map of the DEGs shared by three cold treatment groups

对这29 个基因进行调控网络分析，发现这些基因之间存在着许多调控关系，根据相关程度可划分为4 个相对独立的类群（见图8），/-，（）、，（）基因以及家族基因在调控网络中占据主导地位。

图8 3个低温处理组共有DEGs在低温胁迫下的基因调控网络图圆的直径与基因调控的相关程度成反比；形状大小与基因调控的相关程度成正比；颜色从灰色到黑色表示基因调控的相关程度值从小到大Fig.8 Gene regulation networks among the 29 DEGs shared by three cold treatment groups The diameter of the circle is inversely proportional to the correlation degree of gene regulation.Shape size is proportional to the degree of correlation of gene regulation；The color range from gray to black indicates the degree of correlation of gene regulation from small to large

2.6 基因表达的qPCR分析

为验证萱草转录组数据的可靠性，从DEGs中随机选择了12 个差异显著的基因进行qRT-PCR，该12 个差异表达基因中6 个基因表现为上调（见图9：A～F），另6个基因表现为下调（见图9：G～L）。其中有功能基因，如钙转运ATP 酶基因（，见图9A）；也有调节基因，如含蛋白质MRL7L的硫氧还蛋白基因（--，见图9G）。qRT-PCR结果表明，8个基因的qRT-PCR结果与在转录组内的表达模式基本一致。

图9 12个差异表达基因的qRT-PCR分析柱形表示平均值±标准误，折线表示转录组数据；不同字母表示差异显著，P<0.05Fig.9 Validation of the RNA-Seq results by qRT-PCR of 12 DEGs Column denotes mean±standard error，broken line denotes transcriptome data；Different letters denotes significant difference，P<0.05

3 讨论

低温胁迫对植物的危害众多，主要包括酶活性降低、膜系统遭到破坏、光合作用效率降低、细胞失水等，导致细胞代谢紊乱，生理生化反应受抑制。在植物遭受低温胁迫带来的危害时，植物也会通过激活与冷应答相关的基因，刺激各种信号分子、蛋白质、初级和次级代谢物的变化，从而调节自身的生理生化反应达到降低危害的效果。转录组能够揭示逆境胁迫下整个基因组水平的表达情况，是研究细胞表型和功能的一个重要手段。通过分析植物材料在不同低温处理条件下与对照相比基因表达的差异情况，明确差异基因所位于的代谢途径及生物功能，有助于探索植物低温胁迫分子响应机理。

本研究以萱草‘阿斯隆’为材料，通过转录组测序方法研究了‘阿斯隆’在10、5、0 ℃（分别为T、T、T组）条件下的转录组数据与对照（15 ℃，T组）的差异。共筛选出差异表达基因2 457 个，其中上调基因1 253个，下调基因1 204 个。GO 功能富集分析结果表明，DEGs 主要富集在代谢过程、细胞器、催化活性、转移酶活性等过程，且内膜系统、过氧化物酶活性及天冬氨酸型内肽酶活性等过程为下调基因特有。KEGG 代谢通路分析结果表明，DEGs 在代谢途径、次生代谢产物的生物合成通路的富集程度尤为明显，说明低温对萱草的代谢途径影响较大。低温胁迫下，细胞膜感测到了低温，结构发生变化，细胞膜出现损伤，电解质外渗。此时，植物会通过调整碳水化合物代谢，降低呼吸作用，减少糖分消耗来维持细胞内的稳定。

植物激素是植物防御的主要调节剂之一，复杂的激素信号网络及其控制广泛生理过程的能力使其在帮助植物适应不利环境条件方面起着至关重要的作用。在本研究中植物激素信号转导也是DEGs 富集的一个代谢通路，其中与IAA、CTK和JA 信号通路有关。生长素调节基因的转录激活是生长素在细胞核中的一项快速反应，在IAA 信号通路中，生长素信号转导的第一步是生长素受体与Aux/IAA 的相互作用，从而导致Aux/IAA 降解。在本研究中，Aux/IAA 蛋白编码的DEGs、显著下调，可能参与了萱草叶片响应低温反应。细胞分裂素信号转导途径由磷酸化机制组成，该机制由细胞分裂素与组氨酸激酶受体结合而开始，并随着细胞分裂素反应性基因在细胞核中的转录而达到顶峰。应答调节器（ARR）编码在细胞分裂素应答基因的转录激活中起主要作用。本研究中ARR 编码的2 个基因在低温胁迫下表达量一个上调，一个下调，可能与萱草叶片低温响应有关。在JA 信号通路中，JAZ 蛋白是JA 信号的阻遏物，JA-Ile 通过促进JAZ 蛋白与F-box 蛋白COI1 的结合以及随后通过泛素蛋白连接酶的降解来缓解抑制作用。在本研究中，JAZ 蛋白编码的等4 个DEGS 表达量均表现为下调，且E3 泛素蛋白连接酶（-）为3 个低温处理组共有的差异基因，表达量也显著下调（图7），故推测这些基因在萱草叶片冷应答中起重要作用。

糖在植物发育过程中既具有营养作用，又有信号传导作用。作为信号传导实体，糖可以影响植物生命周期的多种生理过程，并在转录和翻译上起调节作用。糖代谢是整个生物体代谢的中心，连通蛋白质、脂类、核酸及次生物质等代谢。在本研究中，半乳糖代谢、果糖和甘露糖代谢和淀粉和蔗糖代谢通路也是DEGs 富集的几个代谢通路。蔗糖代谢主要包括3 个代谢途径，即Suc分解、己糖磷酸化以及磷酸化己糖的生物合成反应。胞壁转化酶（CWINV4）主要参与韧皮部质外体卸载时蔗糖的分解，保持源—库之间蔗糖的浓度梯度，促使蔗糖顺利卸载。在本研究中，参与Suc 分解的、2 个基因以及参与己糖磷酸化的、-基因均显著上调，且为3 个低温处理组共有的上调差异基因，推测这些基因有极大可能性参与了萱草叶片冷应答过程。参与磷酸化己糖生物合成的基因表达量有所下调，但并不显著，可能参与了萱草叶片低温响应。

将T、T、T中的DEGS 相比较，确定了3 个低温处理组均涉及的29个关键基因。基于基因调控网络分析，在此着重研究在调控网络中相关程度较高并且在上文中未涉及的、以及基因在冷应答中的作用。ABCF5 是ABC转运蛋白家族中的一员，是依赖ATP 的膜结合转运蛋白，能够通过细胞内和细胞外膜转运各种分子，如脂质、蛋白质、化学治疗药物等，对植物响应低温胁迫反应有积极作用。卵形家族蛋白（OFPs）为植物特有，是一类在植物发育中具有重要作用的负调控蛋白，被广泛应用与植物果实形状、根形态等方面的研究，在调节植物生理并赋予其对干旱和寒冷胁迫的抵抗力方面也发挥着重要作用。SWEETs 是一类具有7 个跨膜结构域的糖转运蛋白，该蛋白有着双向转运的功能，可调节葡萄糖跨膜吸收，作为糖转运载体参与糖的运输。基因的调控对生物和非生物胁迫有着多种影响，在拟南芥、茶树以及甘蓝等多种植物上均有基因参与低温响应的相关报道。本研究中，萱草基因在低温胁迫下表达量显著上升，与在番茄（）、拟南芥等研究中的结果一致。故推测这些基因极有可能在萱草叶片响应低温胁迫中起重要的调控作用。